% Single scattering aproximation of palarized radiative transfer eq. 
% kmark Nov 2003
clear
close all
%------------------------------------------------------------------------
% inputs:
AOT=0.05;    % aerosol optical thickness at lambda
lambda=500;  % wavlength in nm
p=1000;      % surface pressure hPa
r=0.25;      % size of aerosol particle
AOTray=rayleigh(lambda,p);

x=2*pi*r/(lambda/1000);  % aerosol size parameter;
n=1.5+0.001i;            % aerosol refractive index
Teta0=60*pi/180;         % Solar zenith angle
Fi0=0*pi/180;            % Solar azimuth angle
Snum=100;
[S1,S2,Qext,Qscat,Qback,g]=mie(x,n,Snum);
SSA=Qscat/Qext;
TetaMie=linspace(0,pi,2*Snum-1);
nT=75;
nF=100;
if mod(nF,2)==0
    nF=nF+1;
end    
Teta=linspace(0,pi/2,nT);
Fi=linspace(0,2*pi,nF);
[T,F]=meshgrid(Teta,Fi);
for i=1:nT
 for j=1:nF
  [i1,i2,CosTETA]=angles(Teta(i),Fi(j),Teta0,Fi0);
  TETA_s=acos(CosTETA);
  Mray=RayMueller(TETA_s);
  Mmie=MieMueller(x,Qscat,TetaMie,S1,S2,TETA_s);
  M=(Mray*AOTray+Mmie*AOT*SSA)/(AOTray+AOT*SSA);
  % sourse function
  II0=1000;
  I0=[II0;0;0;0];  %unpolarized sun radiatiation at the TOA
  %I0=II0/2*[1+cos(Teta0)^2;-1+cos(Teta0)^2;0;0];
  I=exp(-AOT/cos(Teta0))*SSA/(4*pi)*Ltransf(pi-i2)*M*Ltransf(-i1)*I0;
  Pol(j,i)=sqrt(I(2)^2+I(3)^2+I(4)^2)/I(1);  % degree of polarization
  Pol_l(j,i)=sqrt(I(2)^2+I(3)^2)/I(1);  % degree of linear polarization
  Pol_c(j,i)=I(4)/I(1);  % degree of circular polarization

  I1(j,i)=I(1);
  I2(j,i)=I(2);
  I3(j,i)=I(3);
  I4(j,i)=I(4);
 end 
end
% plot results
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,I1);
%caxis([-0.1 0.2]);
%caxis('auto');
shading interp;
caxis([0 1])
colorbar('vert');
%pcolor(abs(Pol)); shading interp; colorbar
figure
py=[flipud(Pol_l(ceil(nF/2),:)')',Pol_l(1,:)];
%px=[Y(fix(nT/2),:),Y(1,:)];
plot(py)