function [Fdif,Fdir,Idif,Fi,Teta]=suc_order(F0,zen,AOT,SSA_aer,g)
% Successive order radiative transfer method
% computes radiance and flux at specific wavelengths
% input
% F0- solar constant
% zen - Solar zenith angle
% AOT - aerosol optical thickness
% SSA_aer - aerosol single scattering albedo
% g - aerosol asumetry parameter
lambda=500;
p=1000;
az=180;
[TAU,SSA,P11]=MIXoptics(lambda,AOT,SSA_aer,1,g,p);

nTeta=30;
nFi=36;
nTau=25;

Teta0=zen*pi/180;
Fi0=az*pi/180;
mi=Gauss_leg(nTeta);         
teta=fliplr(acos(mi));
fi=linspace(0,2*pi,nFi);
ind=find(teta==Teta0);
if isempty(ind)
    teta=[Teta0,teta];
    teta=sort(teta);
    nTeta=nTeta+1;
end
ind=find(fi==Fi0);
if isempty(ind)
    fi=[Fi0,fi];
    fi=sort(fi);
    nFi=nFi+1;
end
tau=linspace(0,TAU,nTau);
[Fi,Teta]=meshgrid(fi,teta);
Idif=zeros(nTeta,nFi);
m0=airmass(90-Teta0*180/pi);
Fdir=F0*cos(Teta0)*exp(-TAU*m0);

num=0;
while 1>0
 num=num+1;   
 cosm0=sin(pi/2-Teta0);
 cosm=sin(pi/2-Teta);
 m0=1./airmass(90-Teta0*180/pi);
 m=1./airmass(90-Teta*180/pi);
 if Teta0>pi/2
    m0=1./airmass(90);
 end
 ind=find(Teta>pi/2);
 m(ind)=-1;
 if num==1
   cosTETA=cosm.*cosm0+sqrt((1-cosm.^2)*(1-cosm0.^2)).*cos(Fi-Fi0);
   [TAU,SSA,Phase]=MIXoptics(lambda,AOT,SSA_aer,cosTETA,g,p);
   I0=F0*exp(tau*m0);  
   for k=1:length(tau)
      J(k,:,:)=SSA/(4*pi)*Phase*I0(k);
   end
 elseif num==2
  Phase=zeros(nTeta,nFi,nTeta,nFi);  
  for i=1:nTeta
    for j=1:nFi
       cosTETA=cosm.*cosm(i,j)+sqrt((1-cosm.^2)*(1-cosm(i,j).^2)).*cos(Fi-Fi(i,j));
       [TAU,SSA,Phase(i,j,:,:)]=MIXoptics(lambda,AOT,SSA_aer,cosTETA,g,p);
    end
  end 
end 
if num>=2 
 for k=1:length(tau)
  II=repmat(squeeze(I(:,:,k)),[1,1,nTeta,nFi]);  %powiela macierz 2D na 4D 
  TERM=SSA/(4*pi)*squeeze(Phase.*II); 
  J(k,:,:)=squeeze(trapz(fi,trapz(teta,TERM)));
 end 
end 
for k=1:length(tau)
   term(k,:,:)=squeeze(J(k,:,:)).*exp(-(TAU-tau(k)).*m).*m;
end
for k=2:length(tau)
   I(:,:,k)=squeeze(trapz(tau(1:k),squeeze(term(1:k,:,:)))); 
end 
 Isurf=I(:,:,end);
 Fdn(num)=-trapz2d(fi, cos(teta),Isurf.*cos(Teta));
 Idif=Idif+Isurf;
 eps=Fdn(end)/sum(Fdn);
 disp([num2str(num),' -st order. EPS : ',num2str(eps)])
 if (num>=5)|(eps<5e-3)
  break
 end
 if num==20
  return
 end
end
Fdif=sum(Fdn);
[X,Y]=pol2cart(Fi,Teta*180/pi);
h=polar([0 2*pi], [0 90],'b');
delete(h);
hold on;

pcolor(X,Y,Idif);
caxis('auto');
shading interp;
colorbar('vert');

%ratio=Fdir./Fdif
