function [Fdn,Fup,Fdir,par]=two_delta(F0,zen,alb,tau,omega,g,TAU,method,adjust);
% Two stream aproximation 
% Based on Fiou "Introduction to atmospheric radiation" page 306
% input:
% F0 - solar flux at specified wavelength
% zen - solar zenith angle
% alb - surface albedo - albedo for flux
% tau - optical depth of output flux
% omega - single scattering albedo
% g - asymmetry parameter
% TAU - optical depth at surface - for boundary condition
% method : "two" or "eddington" see below
% output:
% Fdn - diffuse downward flux
% Fup - diffuse upward flux
% Fdir - downward direct flux
if omega==1
    omega=0.9999999999;
end    
mi=Gauss_leg(1);
m=cos(zen*pi/180);   % cosine of solar zenith angle

if strcmp(adjust,'on');   % delta adjust method
    f=g^2;
    tau=tau*(1-f*omega);
    TAU=TAU*(1-f*omega);
    omega=(1-f)*omega/(1-f*omega);
    g=(g-f)/(1-f);
end    

if strcmp(method,'two')  % standard methods
    gamma1=(1-omega*(1+g)/2)/mi;
    gamma2=omega*(1-g)/(2*mi);
    gamma3=(1-3*g*m*mi)/2;
    if g>mi/m
        disp('Unstable method for small TAU, to large g parameter')
    end    
elseif strcmp(method,'eddington')   % delta eddington methods
    gamma1=(7-(4+3*g)*omega)/4;
    gamma2=-(1-(4-3*g)*omega)/4;
    gamma3=(2-3*g*m)/4;
    if g>0.67/m
        disp('Unstable method for small TAU, to large g parameter')
    end  
end    
Fdir=F0*m*exp(-TAU/m);
k=sqrt(gamma1^2-gamma2^2);
%K=sqrt((1-omega)*(1-omega*g))/mi
par=k*TAU;
Gup=m^2*omega*F0/(m^2*k^2-1)*(gamma3*(gamma1-1/m)+gamma2*(1-gamma3));
Gdn=m^2*omega*F0/(m^2*k^2-1)*((1-gamma3)*(1/m+gamma1)+gamma2*gamma3);

G=(gamma1-k)/gamma2;

%mian=(1-G*alb)*(exp(k*TAU)-G*exp(-k*TAU));
%K1=((alb*Gdn-Gup)*exp(-TAU/m)+Gdn*(1-G*alb)*exp(-k*TAU))/mian
%K2=-(G*(alb*Gdn-Gup)*exp(-TAU/m)+Gdn*(1-G*alb)*exp(k*TAU))/mian

mian=(1-G*alb)*exp(k*TAU)-G*exp(-k*TAU)*(G-alb);
K1=((alb*Gdn-Gup)*exp(-TAU/m)+Gdn*(G-alb)*exp(-k*TAU))/mian;
K2=-Gdn-K1*G;

Fup=K1*exp(k*tau)+G*K2*exp(-k*tau)+Gup*exp(-tau/m);
Fdn=G*K1*exp(k*tau)+K2*exp(-k*tau)+Gdn*exp(-tau/m); 

