%Laser shone at chosen angle into infinitely-deep, water; scattering 
%described by the Fournier-Forand scattering phase function; detect
%backscattered radiation at different radii from incidence point on 
%water surface
clear all;

i = [];

display(['We simulate the propagation of a 670nm-wavelength laser ' ...
'beam through a semi-infinite ocean with known absorption and ' ...
'scattering properties.']);
muz0 = cos(input (['Laser incidence angle, with respect to the ' ...
'normal to the water surface (rad) [0]: '])); %initial direction 
                                              %cosine z
if isempty(muz0);
    muz0 = 1;
end;
display('Light is shone in the x-z plane');
mux0 = sqrt(1-(muz0^2)); %initial direction cosine x
muy0 = 0; %initial direction cosine y (light is shone in x-z plane)

prescribed_Detect_total = ...
    input('Desired detected photon count [30000]: ');
if isempty(prescribed_Detect_total);
    prescribed_Detect_total = 30000; 
end;    
ns = input('Maximum number of scatters per photon packet [200]: ');
if isempty(ns);
    ns = 200;
end;
wmin = input('Minimum photon packet weight [1e-6]: ');
if isempty(wmin);
    wmin = 1e-6;
end;
display(['Scattering will be simulated using the Fournier-Forand ' ...
'scattering phase function']);
Bp = input(['Backscattering ratio (suggested values: 0.005 for ' ...
'phytoplankton-dominated, Case I waters; 0.02 for ' ...
'sediment-dominated, Case I waters; 0.01 for intermediate case): ']);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%establish FF parameters: the particle index of refraction %
%and the slope parameter for the particle size distribution%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (Bp==0.005);
    nparticles = 1.0664;
    mu = 3.3656;
else;
    if (Bp==0.01);
        nparticles = 1.0825;
        mu = 3.4703;
    else;
        if (Bp==0.02);
            nparticles = 1.1028;
            mu = 3.6021;
        else;
            [nu, mu, nparticles, fiterr] = Get_FF_from_Bp(Bp);
        end;
    end;
end;    

%%%%%%%%%%%%%%%%%%%%
%define values of a%
%%%%%%%%%%%%%%%%%%%%
apvec = [0.0013, 0.0026, 0.0053, 0.011, 0.021, 0.043, 0.087, 0.17, ...
         0.35, 0.7]; %particulate absorption values (m^(-1)) - equal 
                     %logarithmic intervals
asw = 0.43; %absorption coefficient for seawater at 670nm (m^(-1))
avec=apvec+asw; %total absorption = particulate + seawater absorpton

bpvec(1:10,1:10) = 0;
cvec(1:10,1:10) = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%define values of bp                                                 %
%particulate scattering coefficient (m^(-1)) - 10 possible values of %
%bp for each value of ap. The ranges correspond to an established    %
%range in ap/bp (see Babin et al., 2003 in Limnology and Oceanography%
%and in Journal of Geophysical Research). For each ap value, the     %
%values of bp are chosen to be at equal logarithmic intervals        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:10;
    bpvec(i,1) = 0.05+6.08*(apvec(i)-0.0013);
    bpvec(i,10) = 0.05+37.63*(apvec(i)-0.0013);
    for j=2:9;
        bpvec(i,j)=10^(log10(bpvec(i,j-1))+(log10(bpvec(i,10))- ...
                   log10(bpvec(i,1)))/9); %equal logarithmic intervals
    end;
end;
bsw=0.0008176; %seawater scattering coefficient at 670nm (m^(-1))
bvec=bpvec+bsw; %total scattering = particulate + seawater scattering
%%%%%%%%%%%%%%%%%%%%
%define values of c%
%%%%%%%%%%%%%%%%%%%%
for i=1:10;
    for j=1:10;
        cvec(i,j)=bvec(i,j)+avec(i);
    end;
end;

[Na, Nc] = size(cvec);

n = input ('Water refractive index [1.34]: ');
if isempty(n);
    n = 1.34;
end;
muc = sqrt(1-1/n^2); %critical angle for internal reflection
rmax = input('Maximum radius of detection (m) [1]: ');
if isempty(rmax);
    rmax = 1;
end;
Nshell = input ('Number of equally-spaced ring detectors [1000]: ');
if isempty(Nshell);
    Nshell = 1000;
end;
deltar = rmax/Nshell; %radial distance between two consecutive 
                      %detector shells (rings); detection "resolution"


%initialize detectors:
Detecttot(1:Na,1:Nc,1:Nshell) = 0;%Detecttot records total photon 
                                  %count/detector shell
Nphotvec(1:Na,1:Nc)=0;
rand('state',sum(100*clock)); %makes sure the initial seed for the 
                              %random number generation is different 
                              %for every run of the program

for j=1:Na;
    a=avec(j);
    for k=1:Nc;
        if(((j==1)&&(k==1))||(j~=1));
        c=cvec(j,k);
        bp=bpvec(j,k);
        b=c-a; %beam scattering coefficient (1/m)
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %define Fournier-Forand cumulative distribution function FF%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
          Psivec = [0.0017452778 0.001971204 0.00219713 0.00248161 ...
                    0.002766091 0.003124222 0.003482353 0.003933158 ...
                    0.004383963 0.004951528 0.005519092 0.006233609 ...
                    0.006948125 0.007847642 0.008747158 0.009879581 ...
                    0.011012005 0.012437635 0.013863265 0.015658021 ...
                    0.017452778 0.01971204 0.021971302 0.024816105 ...
                    0.027660908 0.031242218 0.034823528 0.03933158 ...
                    0.043839633 0.049515276 0.055190919 0.062336086 ...
                    0.069481254 0.078476415 0.087471577 0.098795812 ...
                    0.110120047 0.124376348 0.13863265 0.156580214 ...
                    0.174527778 0.218159722 0.261791667 0.305423611 ...
                    0.349055556 0.3926875 0.436319444 0.479951389 ...
                    0.523583333 0.567215278 0.610847222 0.654479167 ...
                    0.698111111 0.741743056 0.785375 0.829006944 ...
                    0.872638889 0.916270833 0.959902778 1.003534722 ...
                    1.047166667 1.090798611 1.134430556 1.1780625 ...
                    1.221694444 1.265326389 1.308958333 1.352590278 ...
                    1.396222222 1.439854167 1.483486111 1.527118056 ...
                    1.57075000 1.614381944 1.65801389 1.701645833 ...
                    1.745277778 1.788909722 1.832541667 1.876173611 ...
                    1.919805556 1.9634375 2.007069444 2.050701389 ...
                    2.094333333 2.137965278 2.181597222 2.225229167 ...
                    2.268861111 2.312493056 2.356125 2.399756944 ...
                    2.443388889 2.487020833 2.530652778 2.574284722 ...
                    2.617916667 2.661548611 2.705180556 2.7488125 ...
                    2.792444444 2.836076389 2.879708333 2.923340278 ...
                    2.966972222 3.010604167 3.054236111 3.097868056 ...
                    3.1415]; %polar scattering angle values

          nu=(3-mu)/2;
          delta=4*(sin(Psivec/2)).^2/(3*(nparticles-1)^2);
          delta180=4/(3*(nparticles-1)^2);
          cdfFF = (1-delta.^(nu+1)-(sin(Psivec/2)).^2.*(1-delta ...
                  .^nu))./((1-delta).*delta.^nu)+(delta180^nu-1) ...
                  *pi.*cos(Psivec).*(sin(Psivec)).^2./(2*4*pi ...
                  *delta180^nu*(1-delta180)); 
                  %this formula includes polarization effects 
                  %(modified eq.6/p.65/Fournier&Jonasz,1999, by 
                  %replacing w(pi) with b(pi)) 
          cdfTotal = bp*cdfFF/(bp+0.0008176)+2*pi*0.0000509/(bp+ ...
                     0.0008176)*(1+0.835/3-cos(Psivec)-0.835/3* ...
                     (cos(Psivec)).^3);
                     %includes scattering by water - see Alina's 
                     %notes/Sept 29, 2009&Oct 1, 2009
          %Psivec and FF contain 109 values each
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %define look-up vector for establishing the polar scattering%
        %angle                                                      %
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %musvec=cos(Psivec(1,:));
            random_no(1:109)=cdfTotal(1:109);

        display([sprintf('%g',j) ':' sprintf('%g',k) '/' ...
                 sprintf('%g',Na) ':' sprintf('%g',Nc) '---' ...
                 datestr(now) ]);
        %save current data to a file (always overwriting pre-existing 
        %file)
        %- this is a back-up system
        save('exp_name.mat');
    
        Detect_total=0; %Initialize detected photon count
        Nphot=0; %Initialize source photon count
        while (Detect_total<prescribed_Detect_total);
            %%%%%%%%%%%%%%%%%%%
            %initialize photon%
            %%%%%%%%%%%%%%%%%%%
            x = 0; %initial photon position x
            y = 0; %initial photon position y
            z = 0; %initial photon position z
            mux = mux0; %initialize photon direction cosine x
            muy = muy0;
            muz = muz0;
            w = 1; %initial photon weight
            Nphot = Nphot+w;
            nevent = 0; %counts number of events/photon traced
            %%%%%%%%%%%%%%%%%%%

            %%%%%%%%%%%%%%%%%%%%%%
            %air-water interface:%
            %%%%%%%%%%%%%%%%%%%%%%
            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 transmission 
                                            %angle
                ra = ((sin(ti-tt)/sin(ti+tt))^2+(tan(ti-tt)/tan(ti+ ...
                     tt))^2)/2;
                muz = cos(tt); %new direction cosine z
                mux = sqrt(1-(muz^2)); %new direction cosine x
            end;
            w = w*(1-ra); %air-to-water transmission
            %%%%%%%%%%%%%%%%%%%%%%

            while ((nevent<=ns) && (w>wmin)); %keep tracing photon
                %%%%%%%%%%%%
                %scattering%
                %%%%%%%%%%%%
                xold = x; %remember this value for simple calculation 
                          %of photon escape position
                yold = y; %remember this value for simple calculation 
                          %of photon escape position
                zold = z; %remember this value for simple calculation 
                          %of photon escape position
                s = -log(rand)/c; %geometric path length
                x = x+s*mux; %move photon in x direction
                y = y+s*muy; %move photon in y direction
                z = z+s*muz; %move photon in z direction
                zout = z; %remember this value for simple calculation 
                          %of photon escape position
                %%%%%%%%%%%%%%%%%%%%%
                %water-air interface%
                %%%%%%%%%%%%%%%%%%%%%
                if (z<0);
                    muz = - muz; %reflect photon back in the water
                    z = - z;
                    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 
                                                        %transmission 
                                                        %angle
                            ra = ((sin(ti-tt)/sin(ti+tt))^2+(tan(ti ...
                                 -tt)/tan(ti+tt))^2)/2;

                        end;    
                    wout = w*(1-ra); %weight of transmitted photon
                    w = w*ra; %weight of reflected photon    

                    %%%%%%%%%%%%%%%
                    %detect photon%
                    %%%%%%%%%%%%%%%   
                    xs = xold-(zold*(x-xold)/(zout-zold)); 
                         %x position of photon escape (at surface)
                    ys = yold-(zold*(y-yold)/(zout-zold)); 
                         %y position of photon escape (at surface)
                    rs = sqrt((xs^2)+(ys^2));
                    
                    
                    
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    %find out which ring detected the photon - this%
                    %can be made more efficient using smart sorting%
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    if (rs>rmax);
                    else;
                        for j1=1:Nshell;
                            if ((rs>=((j1-1)*deltar)) && (rs<(j1* ...
                                deltar)));
                                Detecttot(j,k,j1)=Detecttot(j,k,j1)+...
                                wout;
                                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                %add detected photon to the detected%
                                %photon count                       %
                                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                Detect_total=Detect_total+wout;
                                break;
                            end;    
                        end;    
                    end;  
                    end;
                end;
                %%%%%%%%%%%%%%%%%%%%%

                w = w*((c-a)/c); %update photon weight by multiplying 
                                 %with single-scattering albedo
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %establish polar scattering cosine%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                Psi=-10; %any impossible value
                die=rand; %random number that will decide the polar 
                          %scattering angle

                left=1;
                right=109;
                mid=fix((left+right)/2);
                while ((Psi==-10)&&(left<mid));
                  if (die==random_no(mid));
                     Psi=Psivec(mid);
                  else;
                     if ((die>random_no(mid-1))&&(die<random_no(mid)));
                        Psi=Psivec(mid-1)+(Psivec(mid)-Psivec(mid- ...
                            1))*(die-random_no(mid- ...
                            1))/(random_no(mid)-random_no(mid-1));
                     else;
                        if (die>random_no(mid));
                           left=mid;
                        else; %die<=random_no(mid-1)
                           right=mid-1;
                        end;
                     end;
                  end;
                  mid=fix((left+right)/2);
                end;
                if (left==mid); %right=left+1
                  if (die==random_no(right));
                    Psi=Psivec(right);
                  else; %die>random_no(left)&die<random_no(right)
                    Psi=Psivec(left)+(Psivec(right)-Psivec(left))* ...
                        (die-random_no(left))/(random_no(right)- ...
                        random_no(left));;
                  end;
                end;
              
                mus = cos(Psi); %Fournier-Forand scattering polar 
                                %scattering cosine
                                
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                phis = 2*pi*rand; %azimuthal scattering angle
                [mux,muy,muz] = chgdir(mux,muy,muz,mus,phis); 
                                %update photon direction; make sure 
                                %file chgdir.m is in the current 
                                %folder!!
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

                nevent = nevent+1;   

            end;   
        Nphotvec(j,k)=Nphot;
        
        end;
        end;
    end;
end;

%normalize Detecttot to initial photon count Nphot:
Nphotaux=repmat(Nphotvec,[1 1 Nshell]);
Detecttot=Detecttot./Nphotaux;

%define normalized photon count/detector shell
for i=1:Nshell;
    Detectnorm(:,:,i)=Detecttot(:,:,i)/((2*i-1)*pi*(deltar^2)); 
    %This value is already normalized to the initial (source) photon 
    %count
end;


%define vector with detector positions (at mid-shell) for data
%visualization and analysis:
Detectr = (deltar/2):deltar:((Nshell-0.5)*deltar);

%save all data to a file
save('exp_name.mat');

display(['Done!---' datestr(now)]);

return;