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

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


x=4.0;  % size parameter;
n=1.5+0.008i;
Snum=100;
[S1,S2,Qext,Qscat,Qback,g]=mie(x,n,Snum);
SSA=Qscat/Qext;
TetaMie=linspace(0,pi,2*Snum-1);
Paer=4*pi/(pi*Qscat*x^2)*0.5*((abs(S1)).^2+(abs(S2)).^2);
tau=AOT+tau_r;
Teta0=30*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:0.03:1,-1:0.03:1);
[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;
  I(j,i)=I0*m0./(m0-m(i)).*omega/(4*pi).*P*(exp(-tau/m0)-exp(-tau/m(i)));
 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');