function E = ssi(sdn, lat, lon)


TauA865 = 0.1;
Angstrom = 1.0;
Dobson = 300.0;

E0    = [170.79, 189.45, 193.66, 188.35, 185.33, 153.41]; % units: mW/cm^2/um
kO3    = [0.0, 3.0e-06, 2.2e-05, 4.1e-05, 9.5e-05, 4.6e-05];
wl     = [412.0, 443.0, 490.0, 510.0, 555.0, 670.0]; %Default: wavelength in nm.  Note that microns used in TauMol calc.

doy = date2jul(sdn - gmtoffset(lon));
% Get Earth-Sun distance factor 'F':
EsFactor = varsol(doy);

% Get Solar Zenith angle:
SolZen = solarzenith(sdn - gmtoffset(lon), lat, lon);

if ( (SolZen < 0.0) | (SolZen > 90.0) )
   fprintf(1, 'Crit. Error in Frouin.posSol: SolZen = %6.2\n', SolZen);
end

% CosSolZen = COSD( SolZen )
CosSolZen = cos(SolZen * pi / 180.0);

%/////////////////////////////////////////////////
%// Estimate key parameters if missing from input:

% Get Dobson climatology if necessary:
%      if (Dobson <= 0) 
%            Dobson = estimDobson(month, lat);
% end


% Get default TauA865 if necessary:
if (TauA865 <= 0) 
   TauA865 = 0.1;
end


% Get default Angstrom coefficient (correspondinng to maritime aerosol)
% if necessary.  Also, convention is for Angstrom > 0 as input, so
% we need to switch signs in that case.

if (Angstrom <= 0) 
   Angstrom = -0.07;
else
   Angstrom = -1.0 * Angstrom;
end

% Setup for transmittance & irradiance equations:
A =  8.435e-03;
B = -1.225e-04;
C =  1.40e-04;

AsymFac = 0.6666667;
Beta    = 0.5 .* ( 1.0 + AsymFac );
Gamma   = 1.0 - AsymFac;

% Calculate albedo of the ocean:
As = 0.05 ./ ( 1.1 .* CosSolZen .^ 1.4 + 0.15 );

% Loop through wavelengths to estimate Surf.Sol.Irrad.:
for i = 1:6
   
   % Compute gaseous transmittance:
   TransG(i) = exp( -kO3(i) .* Dobson ./ CosSolZen);
   
   % Compute diffuse atmospheric transmittance:
   ww  = wl(i) / 1000.0;
   wl4 = ww .^ 4.0;
   wl5 = ww * wl4;
   wl6 = ww * wl5;
   
   TauMol(i) = (A / wl4) + (B / wl5) + (C / wl6);
   TauAer(i) = TauA865 * ((wl(i) / 865.0).^ Angstrom);
   Tau       = TauMol(i) + TauAer(i);
   TransA(i) = exp( -Tau / CosSolZen ) * exp( (0.52*TauMol(i) + Beta*TauAer(i)) / CosSolZen );
   
   % Compute spherical albedo of the atmosphere:
   Sa(i) = exp( -Tau ) * (0.92 * TauMol(i) + Gamma * TauAer(i) );
   E(i) = E0(i) * EsFactor * CosSolZen * TransG(i) * TransA(i) / ( 1.0 - (Sa(i) * As) );
end


function dsol = varsol(doy)
%public static double varSol(int jday, int month) {
om   =(0.9856 * (doy - 4.0)) * pi / 180.0;
dsol = 1 ./ (1 - 0.01673*cos(om).^2);