function [Ftot,Fdir,Fdif]=model(doy,zen,tau500,PW,SSA,g,alfa)
% Spectral semi empirical radiative transfer model
% requires:
%           spect.dat file with molecular extincion coefficients
%           airmass.m
% input: 
%       doy - Julian Day
%       zen - solar zenith angle
%       tau500 - aerosol optical thickness at 500 nm
%       PW - total water vapor in [g/cm^3]
%       SSA - Single scatter albedo at 500 nm
%       g - asymmetry parameter
%       alfa - Angstrom exponent
% output: 
%         Ftot - total flux  at surface in [Wm^-2]
%         Fdir - direct flux  at surface in [Wm^-2]
%         Fdif - diffuse flux at surface in [Wm^-2] 

%Setting parameters
p=1000;   % surface pressure
DU=300;   % total ozone in [DU]
rg=0.05;  % Surface albedo


load spect.dat
lambda=spect(:,1);
Io=spect(:,2);
aH2O=spect(:,3);
aO3=spect(:,4);
agas=spect(:,5);

omega=SSA*ones(size(lambda));
AOT=tau500*(lambda/0.5).^-alfa;


% DIRECT FLUX
M=airmass(90-zen);
p0=1013;
Mp=M*p/p0;
ho=22;
Mo=(1+ho/6370)/sqrt((cos(zen*pi/180)).^2+2*ho/6370);
Taot=exp(-Mp*AOT);
Tray=exp(-Mp.*rayleigh(lambda*1000,p));
TH2O=exp(-0.2385.*aH2O.*PW.*M./(1+20.07.*aH2O.*M).^0.45);
TO3=exp(-DU/1000*aO3*Mo);
Tgas=exp(-1.41.*agas.*Mp./(1+118.3.*agas.*Mp).^0.45);
Itoa=Io.*cos(zen*pi/180).*sundist(doy);
Fdir=Itoa.*Tray.*TH2O.*TO3.*Tgas.*Taot;

% DIFFUSE FLUX

a=find(lambda<=0.45);
Cs=ones(size(Tray));
Cs(a)=(lambda(a)+0.55).^1.8;
ALG=log(1-g);
AFS=ALG.*(1.459+ALG.*(0.1595+ALG*0.4129));
BFS=ALG.*(0.0783+ALG.*(-0.3824-ALG*0.5874));
Fs=1-0.5*exp((AFS+BFS.*cos(zen*pi/180)).*cos(zen*pi/180));
Fs_p=1-0.5*exp((AFS+BFS/1.8)/1.8);
Taot_sca=exp(-omega.*AOT.*M);
Taot_abs=exp(-(1-omega).*AOT.*M);
Iray=Itoa.*TH2O.*TO3.*Tgas.*Taot_sca.*(1-Tray.^0.95)*0.5;

Iaot=Itoa.*TH2O.*TO3.*Tgas.*Taot_abs.*Tray.^1.5.*(1-Taot_sca).*Fs.*Cs;

%rs=TO3_p.*TH2O_p.*Taot_abs_p.*(0.5*(1-Tray_p)+(1-Fs_p).*Tray_p.*(1-Taot_sca_p))

rs=TO3.*TH2O.*Taot_abs.*(0.5*(1-Tray)+(1-Fs_p).*Tray.*(1-Taot_sca));

Igas=(Fdir+Iray+Iaot).*rs.*rg.*Cs./(1-rs.*rg);

Fdif=Iray+Iaot+Igas;
Ftot=Fdir+Fdif;

Fdir=trapz(lambda,Fdir);
Fdif=trapz(lambda,Fdif);
Ftot=trapz(lambda,Ftot);
