function [X,Y,I,x0,y0]=singlescattering(lambda)
%clear
% Single scattering aproximation for solar radiation 
% kmark Sep 2005
%------------------------------------------------------------------------

%lambda=500;
AOT=0.8;
alfa=1;
AOT=AOT*(lambda/500).^-alfa;
tau_r=rayleigh(lambda,1000);


x=4.0;  % size parameter;
n=1.45+0.08i;
Snum=100;
TetaMie=linspace(0,pi,2*Snum-1);
[Paer,ext,SSA]=MieOptics(lambda,n,0.5,2,TetaMie);

%[S1,S2,Qext,Qscat,Qback,g]=mie(x,n,Snum);
%SSA=Qscat/Qext;
%Paer=4*pi/(pi*Qscat*x^2)*0.5*((abs(S1)).^2+(abs(S2)).^2);

tau=AOT+tau_r;
Teta0=51*pi/180;
Fi0=0*pi/180;
[x0,y0]=pol2cart(Fi0,Teta0*180/pi);

nT=75;
nF=50;
if mod(nF,2)==0
    nF=nF+1;
end    
Teta=linspace(0,pi/2,nT);
Fi=linspace(0,2*pi,nF);

%[X,Y]=meshgrid(-1.5:0.07:1.5,-1.5:0.07:1.5);
%[Fi,Teta]=cart2pol(X,Y);
%[nF,nT]=size(X);

[T,F]=meshgrid(Teta,Fi);
m=cos(Teta);
m0=cos(Teta0);
for i=1:nT
 for j=1:nF
  CosTETA=m(i)*m0+sqrt((1-m(i)^2)*(1-m0.^2))*cos(Fi(j)-Fi0);  % cos of 
  TETA_s=acos(CosTETA);
  Pm=3/4*(1+CosTETA^2);
  Pa=interp1(TetaMie,Paer,TETA_s);  
  P=(Pm*1*tau_r+Pa*AOT*SSA)/tau;
  omega=(1*tau_r+SSA*AOT)/tau;
  I0=1000;
  if m(i)<0
      I(j,i)=0;
  else 
   I(j,i)=I0*m0./(m0-m(i)).*omega/(4*pi).*P*(exp(-tau/m0)-exp(-tau/m(i)));
  end   
 end
end
[X,Y]=pol2cart(F,T*180/pi);



%h=polar([0 2*pi], [0 90],'b');
%delete(h);
%hold on;

%pcolor(X,Y,log(I));
%caxis([-0.1 0.2]);
%caxis('auto');
%shading interp;
%caxis([0 1])
%colorbar('vert');
