function ep_plot_dp(ep, Date1, Date2)
% EP_PLOT_DP - Plot derived products of an ep mooring
%
% Use as: ep_plot_dp(ep, c)
%      or ep_plot_dp(ep, Date1, Date2)
% Inputs: ep = Eqpac mooring number(1 or 2)
%         Date1 = String representation of starting date (e.g. '01 jan 1998')
%         Date2 = String representation of ending date (e.g. '01 jan 1998')


% Brian Schlining
% 02 Mar 2000

warning off
ep_ini;                             % Load the header of each column of data (ColumnS) and the units (UnitS)
t       = find(strcmp(ColumnS,'Decimal_Year'));
x       = ep_cmpl;                    % Read the data files in epX directories on tsunami

if ep == 1
   dat = ep1;
else
   dat = ep2;
end

Position = get(0, 'ScreenSize');
Position = [Position(1:3) Position(4) - Position(4)*.1];

if dat.west_longitude < 0
   lat = [num2str(abs(dat.west_longitude)) 'W'];
elseif dat.west_longitude > 0
   lat = [num2str(abs(dat.west_longitude)) 'E'];
else
   lat = '0'
end

if dat.north_latitude < 0
   lon = [num2str(abs(dat.north_latitude)) 'S'];
elseif dat.north_latitude > 0
   lon = [num2str(abs(dat.north_latitude)) 'N'];
else
   lon = '0';
end

loc = [lon ' ' lat];

%=======================================
% Get the dates, Correct for Y2K issues
%=======================================
sdn = x(ep).dat(:,t) + 1900;
sdn = dy2date(sdn);      % Get Dates

%=========================================================
% Refine the approximate dates to make them more accurate
%=========================================================
[Year Month Day Hour Minute Second] = datevec(sdn);
i        = strcmp(ColumnS,'Time');
Hour     = x(ep).dat(:,i);
Hour(find(Hour == 1200)) = 12/24;
Hour(find(Hour == 1230)) = 12.5/24;
Day      = round(Day) + Hour;  % Correct the day to the local time
sdn      = datenum(Year, Month, Day);

if nargin < 3
   Date1 = min(sdn);
   Date2 = max(sdn);
end

goodTime = find(sdn >= Date1 & sdn <= Date2); % Find Points within the Date Limits
if isempty(goodTime)
   goodTime = length(sdn);
end


%=================
% Check arguments
%=================
if nargin == 1
   Date1 = min(sdn);
   Date2 = max(sdn);
end

%==========================
% Plot attenuation coef (K). This is acutally KEd but we apply it to Lu 
%==========================
fh = figure;
set(fh, 'Position', Position)
a    = strmatch('K',ColumnS);
b    = strmatch('KL',ColumnS);
KCol = setdiff(a,b);
plot(sdn(goodTime), x(ep).dat(goodTime, KCol));
ylabel('K [m^-^1]')
xlabel('Date')
title(['Attenuation Coeff. at ' loc])
datetick('x')
legend(ColumnS{KCol}, 0)
setfontsize(7)
set(gca,'YLim',[0 1]);
grid

%=========================
% Plot attenuation length
%=========================
fh = figure;
set(fh, 'Position', Position)
plot(sdn(goodTime), 1./x(ep).dat(goodTime, KCol));
ylabel('1/K [m]')
xlabel('Date')
title(['Attenuation length at ' loc])
datetick('x')
legend(ColumnS{KCol}, 0)
setfontsize(7)
set(gca,'YLim',[1 100], 'YScale', 'log');
grid

%=================================
% Plot Water-leaving radiance (Lw)
%=================================
fh = figure;
set(fh, 'Position', Position)
LwCol = strmatch('LLw',ColumnS);
plot(sdn(goodTime), x(ep).dat(goodTime, LwCol));
ylabel('Lw [\muW cm^-^2 nm^-^1 str^-^1]')
xlabel('Date')
title(['Water-Leaving radiance at ' loc])
datetick('x')
legend(ColumnS{LwCol}, 0)
setfontsize(7)
set(gca,'YLim',[0 1.2]);
grid

%==========
% Plot Rrs
%==========
fh = figure;
set(fh, 'Position', Position)
a = strmatch('Ed[3m+]',ColumnS);
b = strmatch('Ed[3m+]PAR',ColumnS);
EsCol = setdiff(a,b);
EsCol = EsCol(1:5);
lambda = [412 443 490 510 555];
plot(sdn(goodTime), x(ep).dat(goodTime, LwCol(1:5))./x(ep).dat(goodTime, EsCol),'o')
title(['Remote Sensing Reflectance at ' loc])
ylabel('Rrs [str^-^1]')
xlabel('Date')
datetick('x')
legend(num2str(lambda'), 0)
setfontsize(7)
grid

%==============
% 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
Lw555 = x(ep).dat(goodTime, LwCol(5));
R15 = log10(x(ep).dat(goodTime, LwCol(1))./Lw555); 
R25 = log10(x(ep).dat(goodTime, LwCol(2))./Lw555); 
R35 = log10(x(ep).dat(goodTime, LwCol(3))./Lw555); 
R45 = log10(x(ep).dat(goodTime, LwCol(4))./Lw555); 

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