function ea = oms_esanom(infile, DateLim)
% OMS_ESANOM - Calculate Es anomalies from frouin's model
%
% Use as: ea = oms_esanom(infile, DateLim)
%
% Inputs: infile = a netcdf object
%         Datelim = 2-element vector containing [minDate maxDate] of inclusive data to plot.
%                   Dates must be in matlab's DATENUM format.
%
% Output: ea = matrix of Es anomalies from Frouin's model.

% Brian Schlining
% 30 Mar 2000

lambda = [412 443 490 510 555 670];

OK = oms_checkfile(infile);

t = nctime(infile);
if nargin == 2
   minIndex = min(find(t >= DateLim(1) & t <= DateLim(2)));
   maxIndex = max(find(t >= DateLim(1) & t <= DateLim(2)));
   good = minIndex:maxIndex;
else
   good = 1:length(t);
end

t  = t(good);
Es = ones(length(t), length(lambda))*NaN;
S  = oms_getinfo(infile);
v  = var(infile);



%===================
% Find Es at lambda
%===================
iEs = strmatch('es',S.symbol);               % Find Es
for i = 1:length(lambda)

   iWp = find(S.wavelength(iEs) == lambda(i));     % Find Es at the desired wavelength
   goodEs = iEs(iWp);
   if isempty(goodEs)
      warning(['Es[' num2str(lambda(i)) 'nm] (Surface irradiance) sensors were not found']);
      Es(:,i) = NaN;
   else

      vEs     = v{goodEs};                         % Get the correct Es variable object
      Es(:,i) = vEs(good);
      bad     = find(Es(:,i) == vEs.FillValue_(:));
      Es(bad) = NaN;
      bad     = find(Es(:,i) == vEs.missing_value(:));
      Es(bad) = NaN;
   end

end


E0 = frouin(t, infile{'latitude'}(:), infile{'longitude'}(:));

ea = Es - E0;
bad = find(ea < -200);
ea(bad) = NaN;

plot(t, ea(:,1))
hold on
plot(t, ea(:,2),'c')
plot(t, ea(:,3),'g')
plot(t, ea(:,4),'m')
plot(t, ea(:,5),'r')
x = oms_deploys(char(infile.mooring(:)));
plot([x x]', repmat(get(gca,'YLim'), length(x), 1)','k:')
hold off
datetick('x')
legend('412nm','443nm','490nm','510nm','555nm',0)
xlabel('Year')
ylabel('Es Anomaly [\muW cm^-^2 nm^-^1]')
title(['Es Anomalies at ' char(infile.mooring(:))])
