% Sucessive order rative trensfer method for shortwave and black Earth's surface
% Solve the radiative transfer problem for minochromatic wavelength including
% molecular scattering and aerosol absorption and scattering processes
clear
close all
% Input parameters
F0=1368;       % Solar constant
lambda=500;    % Wavelength in nm
Teta0=30*pi/180;  % Solar Zenith Angle in radians
Fi0=180*pi/180;   % Solar azimuth angle in radians
p=1000;           % Surface pressure in hPa
% Definition of aerosol optical properties-------------------
AOT=0.2;      % aerosol optical depth
SSA_aer=0.95; % aerosol scattering albedo
g=0.7;        % aerosol asymmetry parameter
%------------------------------------------------------------
% computes optical properties of aerosol and molecule mixture
[TAU,SSA,P11]=MIXoptics(lambda,AOT,SSA_aer,1,g,p);

nTeta=30;
nFi=72;
nTau=25;

teta=linspace(0,pi/2,nTeta);
fi=linspace(0,2*pi,nFi);
tau=linspace(0,TAU,nTau);
[Fi,Teta]=meshgrid(fi,teta);
Idif=zeros(nTeta,nFi);
m0=airmass(90-Teta0*180/pi);
Fdir=F0*cos(Teta0)*exp(-TAU*m0);

num=0;
% main loop regarding the sucessive order 
while 1>0
 num=num+1;   
 cosm0=sin(pi/2-Teta0);
 cosm=sin(pi/2-Teta);
 m0=1./airmass(90-Teta0*180/pi);
 m=1./airmass(90-Teta*180/pi);
 if Teta0>pi/2
    m0=1./airmass(90);
 end
 ind=find(Teta>pi/2);
 m(ind)=-1;
 cosTETA=cosm.*cosm0+sqrt((1-cosm.^2)*(1-cosm0.^2)).*cos(Fi-Fi0);
 [TAU,SSA,Phase]=MIXoptics(lambda,AOT,SSA_aer,cosTETA,g,p);
 for i=1:nTeta
    for j=1:nFi
     if num==1   % The first orrder 
      I0=F0*exp(tau*m0);
      for k=1:length(tau)
       J(k)=SSA/(4*pi)*Phase(i,j)*I0(k);
      end
      else
       cosTETA=cosm.*cosm(i,j)+sqrt((1-cosm.^2)*(1-cosm(i,j).^2)).*cos(Fi-Fi(i,j));
       [TAU,SSA,Phase]=MIXoptics(lambda,AOT,SSA_aer,cosTETA,g,p);
       for k=1:length(tau)
        TERM=SSA/(4*pi)*Phase.*I; 
        J(k)=trapz2d(fi,teta,TERM); 
       end 
      end 
      term=J.*exp(-(TAU-tau).*m(i,j)).*m(i,j);
      I(i,j)=trapz(tau,term);
   end
 end
 Fdif(num)=-trapz2d(fi, cos(teta),I.*cos(Teta));
 Idif=Idif+I;
 eps=Fdif(end)/sum(Fdif);
 disp([num2str(num),' -th order. EPS : ',num2str(eps)])
 if (num>=5)|(eps<1e-3)
  break
 end
end

% Simple plot of the sky radiance
[X,Y]=pol2cart(Fi,Teta*180/pi);
h=polar([0 2*pi], [0 90],'b');
delete(h);
hold on;

pcolor(X,Y,Idif);
caxis('auto');
shading interp;
colorbar('vert');
