%Code from Leathers et al., 2004, NRL report
% Optically Deep Ocean
% Monte Carlo Matlab script for computing computing upwelling and downwelling
%irradiances in a homogeneous ocean. Scattering is described by the Henyey-Greestein scattering phase
%function


c = 1; % attenuation coefficient (1/m)
w0 = 0.8; % single-scattering albedo
g = 0.9; % scattering asymmetry factor
zd = 1.0; % depth of downwelling detector
zu = 1.0; % depth of upwelling detector
n = 1.34; % relative index of refraction for water
N = 1e4; % number of photons to trace
ns = 50; % number of scatters to trace
muc = sqrt(1-1/n^2); % critical angle for internal reflection
Eu=0; Ed=0; % initiate detectors
for i=1:N
z = 0; % initial position depth
muz = sqrt(rand); % diffuse illumination
w = 1; % initial photon weight
if (muz==1),
ra= ((n-1)/(n+1))^2; % reflection normal to surface
else
ti = acos(muz); % angle of incidence
tt = asin(sqrt(1-muz^2)/n); % air-to-water trans. angle
ra = ((sin(ti-tt)/sin(ti+tt))^2+(tan(ti-tt)/tan(ti+tt))^2)/2;
muz = cos(tt); % new direction
end
w = w*(1- ra); % air-to-water transmission
for j=1:ns
zold=z; % previous position
s = -log(rand)/c; % geometric path length
z = z + muz*s; % move photon
if (z<zu)&(zold>zu), % photon passes Eu sensor
Eu=Eu+w; % update upwelling sensor
end
if (z<0) % photon reaches surface
muz = - muz; % change to positive value
z = - z; % new position of reflected photon
zold = 0; % photon now traveling down
if (muz>muc), % partial transmission
    if (muz==1), % photon vertical
ra= ((n-1)/(n+1))^2; % reflection normal to surface
else
ti = acos(muz); % angle of incidence
tt = asin(n*sqrt(1-muz^2)); % water-to-air trans. angle
rw = ((sin(ti-tt)/sin(ti+tt))^2+(tan(ti-tt)/tan(ti+tt))^2)/2;
w = w*rw; % weight of reflected photon
end
end
end
if (z>zd)&(zold<zd), % photon passes Ed sensor
Ed=Ed+w; % update downwelling sensor
end
w=w*w0; % in-water scattering
mu = (1+g^2-((1-g^2)/(1+g-2*g*rand))^2)/(2*g); % HG scattering
Phi = 2*pi*rand; % azimuthal scattering angle
muz = muz*mu - sqrt(1-muz^2)*sqrt(1-mu^2)*cos(Phi); % new direction
end
end
disp(sprintf('Eu at %g m: %g',zu,Eu/N));
disp(sprintf('Ed at %g m: %g',zd,Ed/N));