function E = frouinPlus(lambda, t, lat, lon)
% FROUIN   - Incident surface irradiance model for SeaWiFS bands
%
% Use as: E = frouinPlus(lambda, t, lat, lon)
%
% Inputs: lambda = wavelengths to return (nm)
%         t   = GMT time (in matlab's date format, see DATENUM and DATEVEC)
%         lat = Latitude (+N/-S)
%         lon = longitude (+E/-W)
%
% Output: E = incoming clear-sky solar irradiance just above the surface 
%             at wavelengths specifyied by lambda. 
%             Each row of E corresponds to the time in t, the columns
%             correspond to the wavelengths( i.e Column 1 is lambda(1), 2 is lmbda(2)...).
%             E will contain NaN's for times when the sun is below the horizon.
%
% FROUIN calculates the incoming clear-sky solar irradiance
%	just above the surface at SeaWiFS wavelengths (412, 443, 490, 510,
%	555, 670 nm), using this equation:
%
%	E = E0 * ESfactor * COS(SolZen) * TransG * TransA / (1 - (Sa * As))
%
%	  where:
%	E0       = extra-terrestrial solar irradiance (mW/cm^2/um)
%	ESfactor = Earth-Sun distance factor
%	SolZen   = solar zenith angle (deg).
%	TransG   = gaseous transmittance
%	TransA   = diffuse atmospheric transmittance
%	Sa       = spherical albedo of the atmosphere
%	As       = albedo of the ocean
%	[note: Routine estimates SolZen from input parameters.
%	 Valid for 0 -- 90 deg, but for SolZen > 80 deg,
%	 the above formula and the approximations used for the various
%	 atmospheric functions are not very accurate.]
%
% Frouin is a stright port of a FORTRAN program to Java. This
% matlab function uses matlab's java interface to call the java 
% class that does the calculations. The original FORTRAN source 
% can be found at: http://genius.ucsd.edu/~john/SeaWiFS_dir/ssi.f
%
% The java code does allow inputs of TauA865(aerosol optical thickness, 865 nm),
%	Angstrom (Angstrom coefficient (convention: Angstrom > 0 )),
%  and Dobson (ozone amount in Dobson units). This functionality is not 
%  currently present in this matlab interface but may be added if needed.
%
% To Install: Copy Frouin.class to a local directory. Then add a line to 
% \[matlabroot]\toolbox\local\classpath.txt pointing to the directory.
% For example if Frouin.class is in the C:\matlab\work directory, the line
% added should be c:\matlab\work\. It will be necessary to 
% restart matlab for the changes in the classpath to take effect.
% 
%
% References:
%
% Frouin, R., D. W. Ligner, and C. Gautier, 1989: A Simple analytical formula
% to compute clear sky total and photosynthetically available solar irradiance
% at the ocean surface. J. Geophys. Res., 94, 9731-9742.
%
% Tanre, D., M. Herman, P.-Y. Deschamps, and A. De Leffe, 1979: Atmospheric
% modeling for Space measurements of ground reflectances, including
% bi-directional properties. Appl. Optics, 18, 21,3587-21,3597.

% Brian Schlining
% 30 Mar 2000

java on;

nT = length(t);
nLat = length(lat);
nLon = length(lon);

if nLat ~= nLon
   E = [];
   error('  Latitude and Longitude inputs are not the same size');
   return;
elseif nT ~= nLat
   if nT > 1 & nLat > 1
      error('  Inputs vectors are not the same size');
   elseif  nLat == 1
      lat = ones(size(t))*lat;
      lon = ones(size(t))*lon;
   else
      t = ones(size(lat))*t;
   end
   
end

n = length(t);

t2 = t;
% Java uses time as milliseconds since epoch (01 Jan 1970 00:00:00). 
t = round((t - datenum('01 Jan 1970 00:00:00'))*60*60*24*1000);
E = ones(n , length(lambda))* NaN;
fim = org.mbari.model.FrouinIrradianceModel(lambda);
%fim = java_object('org.mbari.model.FrouinIrradianceModel', lambda);
for i = 1:length(t)
   try 
      E(i,:) = fim.calcIrradiance(t(i), lat(i), lon(i))';
   catch
      E(i,:) = ones(1, length(lambda))* -999;
   end
end

E(find(E == -999.9 | E == -999)) = NaN;