function [Fdn,Fup]=two_stream(F0,zen,tau,omega,g,TAU);
% Two stream aproximation - with zero surface albedo
% Based on Fiou "Introduction to atmospheric radiadion" page 306
% F0 - solar flux at specified wavelength
% zen - solar zenith angle
% tau - optical depth of output flux
% omega - single scattering albedo
% g - asymmetry parameter
% TAU - optical depth at surface - for boundary condition

mi=Gauss_leg(1);
m=cos(zen*pi/180);   % cosine of solar zenith angle
k=sqrt((1-omega)*(1-omega*g)/mi^2);
Sup=F0*omega/(4*pi)*(1+3*g*mi*m);
Sdn=F0*omega/(4*pi)*(1-3*g*mi*m);
Z1=-(1-omega*g)*(Sdn+Sup)/mi^2+(Sdn-Sup)/(mi*m);
Z2=-(1-omega)*(Sdn-Sup)/mi^2+(Sdn+Sup)/(mi*m);
alfa=Z1*m^2/(1-m^2*k^2);
beta=Z2*m^2/(1-m^2*k^2);
epsilon=(alfa+beta)/2;
gamma=(alfa-beta)/2;
a=sqrt((1-omega)/(1-omega*g));
v=(1+a)/2;
u=(1-a)/2;
K=-(epsilon*v*exp(-TAU/m)-gamma*u*exp(-k*TAU))/(v^2*exp(k*TAU)-u^2*exp(-k*TAU));
H=-(epsilon*u*exp(-TAU/m)-gamma*v*exp(-k*TAU))/(v^2*exp(k*TAU)-u^2*exp(-k*TAU));

Iup=K*v*exp(k*tau)+H*u*exp(-k*tau)+epsilon*exp(-tau/m);
Idn=K*u*exp(k*tau)+H*v*exp(-k*tau)+gamma*exp(-tau/m);

Fdn=2*pi*mi*Idn;
Fup=2*pi*mi*Iup;