clear
% Single scattering aproximation for solar radiation 
% kmark Sep 2005
%------------------------------------------------------------------------
% Geometry
m=1;   % cosine of zenith angle
Fi=0*pi/180;
m0=0.3;
Fi0=180*pi/180;
% 
lambda=760; 
% Aerosol properties
ref=1.5+0.001i;
rm=0.5;
sigm=2;
AOT=0.02;
DU=200;

[z,p,T,ext_a,ozone]=profiles;
nlev=length(z);

CosTETA=m*m0+sqrt((1-m^2)*(1-m0))*cos(Fi-Fi0);  % cos of scattering angle
TETA=acos(CosTETA);

ext_m=molec(lambda,T,p);
Pm=3/4*(1+CosTETA^2);
aot=trapz(z,ext_a);
ext_a=ext_a*AOT/aot;


du=trapz(z,ozone);
ozone=ozone*DU/1000/du;
load spect.dat
kO3=interp1(spect(:,1),spect(:,4),lambda/1000);
ext_O3=kO3*ozone;

agas=interp1(spect(:,1),spect(:,5),lambda/1000);
kO2=1.41*agas/(1+118.3*agas)^0.45;
ext_O2=ext_m*kO2/trapz(z,ext_m);

I0=interp1(spect(:,1),spect(:,2),lambda/1000);

[Pa,ext,ssa]=MieOptics(lambda,ref,rm,sigm,TETA);
[P,ext,SSA]=MieRayleigh(Pa,Pm,ext_a,ext_m,ssa,ext_O3,ext_O2);
for i=1:nlev-1
 tau(i)=trapz(z(i:nlev),ext(i:nlev));
end 
tau(nlev)=0;
TERM=SSA.*P.*exp(tau./m).*exp(-tau/m0);
I=-I0*m0*exp(-tau(1)/m)/(4*pi*m)*trapz(tau,TERM)


return

h=polar([0 2*pi], [0 90],'b');
delete(h);
hold on;
[X,Y]=pol2cart(F,T*180/pi);
%pcolor(X,Y,Pol_l);
pcolor(X,Y,log(I));
%caxis([-0.1 0.2]);
%caxis('auto');
shading interp;
%caxis([0 1])
colorbar('vert');
