function Rrs = oms_plotRrs(infile, DateLim)
% oms_plotRrs - Plot Rrs data from an OMS object
%
% Use as: Rrs = oms_plotRrs(infile, DateLim)
%
% Inputs: infile = a netcdf OMS object
%         DateLim = 2 element vector specifying the time range to plot (optional)

% Brian Schlining
% 28 Jul 2000

lambda = [412 443 490 510 555];
S = oms_getinfo(infile);
t = nctime(infile);
v = var(infile);
Es = ones(length(t), length(lambda))*NaN;
Lw = Es;
useDepth = 20;

if nargin < 2
   DateLim = [min(t) max(t)];
end
goodTime = find(t >= DateLim(1) & t <= DateLim(2)); % Find Points within the Date Limits
if isempty(goodTime)
   goodTime = length(sdn);
end

loc = ' ';


for i = 1:length(lambda)
   
   %===================
   % Find Ed at lambda
   %===================
   iEs       = strmatch('es',S.symbol);					% Find Es
   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']);
   end
   
   vEs       = v{goodEs};									% Get the correct Es variable object
   Es(:,i)   = vEs(:);
   bad       = find(Es(:,i) == vEs.FillValue_(:));
   Es(bad,i) = NaN;
   
   
   %===================
   % Find Lw at lambda
   %===================
   [lw k n] = oms_lw(infile, lambda(i));
   if ~isempty(lw)
      if length(n) > 1
         for j = 1:length(n)
            z = oms_getdepth(infile, infile{n{j}});
            if isempty(z)
               Z(j) = NaN;
            else
               Z(j) = z(1);
            end
            
         end
         good = near(Z, useDepth);
         Lw(:, i) = lw(:, good);
      else
         Lw(:,i) = lw;
      end
      
   end  
   
end

Rrs = Lw./Es


%==============
% Plot OC2 CHL
%==============
%Coefficients of the OC2 MCP equations in O'Reilly et al. 1998 
c15 = [0.2457, -1.7620, 0.2830, 0.1035, -0.0388]; 
c25 = [0.1909, -1.9961, 1.3020, -0.5091, -0.0815]; 
c35 = [0.3410, -3.0010, 2.8110, -2.0410, -0.0400]; 
c45 = [0.4487, -4.3665, 2.7130, -0.2698, -0.0821]; 

%Compute ratios and log-transform
R15 = log10(Rrs(:,1)./Rrs(:,5)); 
R25 = log10(Rrs(:,2)./Rrs(:,5)); 
R35 = log10(Rrs(:,3)./Rrs(:,5)); 
R45 = log10(Rrs(:,4)./Rrs(:,5)); 

% Calculate chlorophyll concentration (ug/l) 
C_oc15 = 10.0.^(c15(1) + c15(2).*R15 +  c15(3).*R15.^2 + c15(4).*R15.^3) + c15(5); 
C_oc25 = 10.0.^(c25(1) + c25(2).*R25 +  c25(3).*R25.^2 + c25(4).*R25.^3) + c25(5); 
C_oc35 = 10.0.^(c35(1) + c35(2).*R35 +  c35(3).*R35.^2 + c35(4).*R35.^3) + c35(5);
C_oc45 = 10.0.^(c45(1) + c45(2).*R45 +  c45(3).*R45.^2 + c45(4).*R45.^3) + c45(5);

% Coefficeint of variance
cv1  = std([C_oc15 C_oc25 C_oc35 C_oc45]')'./mean([C_oc15 C_oc25 C_oc35 C_oc45]')';

fh = figure;
%set(fh, 'Position', Position)
subplot(2,1,1)
plot(t(goodTime),cv1)
title(['OC2 CHL at ' loc])

ylabel('Coeff. Variation (std/mean)')
datetick('x') %,12)
set(gca, 'YLim', [0 2]);
grid

subplot(2,1,2)
p(1) = plot(t(goodTime), C_oc15, 'ro');
hold on
p(2) = plot(t(goodTime), C_oc25, 'bo');
p(3) = plot(t(goodTime), C_oc35, 'co');
p(4) = plot(t(goodTime), C_oc45, 'ko');
%p(5) = plot(t(goodTime), morel_chl2(, 490), '-', 'Color', [.8 .42 .2], 'lineWidth', 1.5);
hold off
datetick('x') %,12);
ylabel('OC2 [CHL]');
xlabel('Date')
legend(p,'R15','R25','R35','R45',0);
set(gca, 'YLim', [0.01 100], 'YScale','log')
grid
setfontsize(7)





