% Run Simple Monte Carlo model for radiative transfer in the atmosphere
% Single column and single wavlength case, black Earth's surface
% This version of model inludes only aerosol particle in the atmosphere
clear
close all

% Inputs:
Teta=0;          % Solar Zenith Angle
Fi=0;            % Solar Azimuth Angle
nfoton=5000;    % number of photons
nlev=100;        % number of layers

AOT=0.1;         % aerosol optical depth
g=0.7;           % asymetry parameter
SSA=0.9;         % single scattering albedo

Teta=Teta*pi/180;
Fi=Fi*pi/180;
cosTeta=cos(Teta);
sinTeta=sin(Teta);
cosFi=cos(Fi);
sinFi=sin(Fi);

z=0;
dtau=AOT/nlev;     % optical depth of each layer
zz=[0:dtau:AOT];
T=zeros(nfoton,1);
Z=zeros(nfoton,1);
K=zeros(nfoton,1);
F=zeros(nfoton,1);
Edtau=zeros(nlev,1);
Euz=zeros(nlev,1);
j=0;
for i=1:nfoton     % main loop
% disp(i/nfoton*100)
 z0=0;
 nk=0;
 cosTeta_new=cosTeta;
 sinTeta_new=sinTeta;
 cosFi_new=cosFi;
 sinFi_new=sinFi;
 while z<10000
  nk=nk+1;
  s=-log(rand);
  z=z0+s*cosTeta_new;  
  if z<=0
   j=j+1;   
   Y(i)=1;
   T(i)=cosTeta_new;
   F(i)=cosFi_new;
   if z0>0
    az=find(zz>z0);
    nz=az(1);
    Euz(1:nz-1)=Euz(1:nz-1)+1;
   end
   break   % photon escape the atmosphere
  elseif z>=AOT
   Z(i)=AOT;
   K(i)=nk;
   T(i)=cosTeta_new;
   F(i)=cosFi_new;
   az0=find(zz>z0);
   nz0=az0(1);
   Edtau(nz0-1:nlev)=Edtau(nz0-1:nlev)+1;
   break   % photon hit the surface
  else
  r=rand;
   if r>SSA
    Z(i)=z;
    az0=find(zz>z0);
    nz0=az0(1);
    az=find(zz>z);
    nz=az(1);
    if cosTeta>0
      Edtau(nz0-1:nz-1)=Edtau(nz0-1:nz-1)+1;
    else
      Euz(nz-1:nz0-1)=Euz(nz-1:nz0-1)+1;
    end
   break   % photon absorbed
   elseif r<=SSA
    cosTeta0=gen_HG(1,g); %Teta0 changes from 0 to pi
    %cosTeta0=gen_Ray(1);
    %cosTeta0=(2*rand-1);   % uniform distrubiton 
    sinTeta0=sqrt(1-cosTeta0^2);
    Fi0=rand*2*pi;
    cosFi0=cos(Fi0);
    sinFi0=sin(Fi0);
    [cosTeta_new,sinTeta_new,cosFi_new,sinFi_new]=rot_coordinate(cosTeta_new,sinTeta_new,cosFi_new,sinFi_new,cosTeta0,sinTeta0,cosFi0,sinFi0);                             
    az0=find(zz>z0);
    nz0=az0(1);
    az=find(zz>z);
    nz=az(1);
    if cosTeta>0
      Edtau(nz0-1:nz-1)=Edtau(nz0-1:nz-1)+1;
    else
      Euz(nz-1:nz0-1)=Euz(nz-1:nz0-1)+1;
    end
    z0=z;
   end  
  end %photon scattered
 end % while loop
end % nfotons

%save carlo1 Z T F AOT

Kd=-diff(log(Edtau))/dtau;
Ku=-diff(log(Euz))/dtau;

tdir=length(find(K==1));
ttot=length(find(Z==AOT));
tu=length(find(Z==0));
ta=nfoton-(ttot+tu);

Idir=tdir/nfoton;
Itot=ttot/nfoton;
Idif=Itot-Idir;
Itoa=tu/nfoton;
Iabs=ta/nfoton;

% uncertainty
dIdir=sqrt((nfoton-tdir)/(nfoton*tdir));
dItot=sqrt((nfoton-ttot)/(nfoton*ttot));
tdif=ttot-tdir;
dIdif=sqrt((nfoton-tdif)/(nfoton*tdif));
dIabs=sqrt((nfoton-ta)/(nfoton*ta));
dItoa=sqrt((nfoton-tu)/(nfoton*tu));

ad=find(Z==AOT);
T1=acos(T);
for i=1:18;
 Ad=find((T1(ad)>=5*(i-1))&(T1(ad)<5*i));
end

disp('------------------------------------------------')
disp ('*******************Results*********************')
disp('------------------------------------------------')
disp(['Direct normalized radiation = ',sprintf('%5.5f',Idir),' +- ',sprintf('%5.5f',dIdir)])
disp(['Diffuse normalized radiation = ',sprintf('%5.5f',Idif),' +- ',sprintf('%5.5f',dIdif)])
disp(['Total normalized radiation = ',sprintf('%5.5f',Itot),' +- ',sprintf('%5.5f',dItot)])
disp(['TOA reflect normalized radiation = ',sprintf('%5.5f',Itoa),' +- ',sprintf('%5.5f',dItoa)])
disp(['Normalized absorbing radiation = ',sprintf('%5.5f',Iabs),' +- ',sprintf('%5.5f',dIabs)])

