function [X,Y,I,Idir,x0,y0,F,T,I0,BRF]=singlescattering(lambda,zen,az,AOT)
% Single scattering aproximation sover 
% kmark Sep 2005

alfa=1;
Snum=100;
TetaMie=linspace(0,pi,2*Snum-1);

if 3>2  % MIE camputation
 % x=4.0;  % size parameter;
  n=1.45+0.002i;    % aerosol refractive index
  rm=0.05;          % mode radius for log-normal size distribution
 [Paer,ext,SSA]=MieOptics(lambda,n,rm,2,TetaMie);
 [Paer1,ext1,SSA1]=MieOptics(500,n,rm,2,TetaMie);
 if lambda==500
   alfa=1;
 else
  alfa=abs(log(ext/ext1)/(log(lambda/500)))
 end
else  % definition of aerosol optical properties without MIE computation
    SSA=0.95;
    g=0.7;
    Paer=HGfun(TetaMie*180/pi,g);
end    

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

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

nT=300;
nF=500;
Teta=linspace(0,pi/2,nT);
Teta=sort([Teta0,Teta]);
nT=nT+1;
phi=linspace(0,2*pi,nF);
phi=sort([phi,Fi0]);
nF=nF+1;

m0=cos(Teta0);
M=cos(Teta);
[m,Fi]=meshgrid(M,phi);
I=zeros(nF,nT);
load solar
dl=5;
index=find((sol.lambda>=lambda-dl)&(sol.lambda>=lambda+dl));
I0=mean(sol.I0(index));
%I0=1368;

CosTETA=m.*m0+sqrt((1-m.^2)*(1-m0.^2)).*cos(Fi-Fi0);
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;

ind=find(m>0);
TRANS=(exp(-tau/m0)-exp(-tau./m));
I(ind)=I0*m0./(m0-m-eps).*omega/(4*pi).*P.*TRANS;
BRF=1./(m0+m+eps).*omega/4.*P.*(1-exp(-tau./m-tau./m0));
ind=find(TRANS==0);
if isempty(ind)==0
   I(ind)=I0./(m0+eps).*omega/(4*pi).*P(ind).*tau*exp(-tau/m0); 
end    

[T,F]=meshgrid(Teta,phi);
[X,Y]=pol2cart(F,T*180/pi);

Idir=I0*exp(-tau/m0);
Fdif=-trapz2d(phi,cos(Teta),(I.*cos(T)));
ratio=Fdif/(Idir*cos(pi/180*Teta0));

return





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);
