function oms_qcplot(infile, t, Filename)
% OMS_QCPLOT - Plot a single day of rediometer data against modeled data
%
% Use as: oms_qcplot(infile, t)
%         oms_qcplot(infile, t, Filename)
%
% Inputs: infile = an OMS netcdf object
%         t      = date of day to plot (in matlab's datefomrat, see DATENUM)
%         Filename = optional argument specifying the name of the gif image to create
%

% Brian Schlining
% 08 Mar 2000

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

lat = infile{'latitude'}(:);
lon = infile{'longitude'}(:);

if lat < 0
   latS = [num2str(abs(lat)) 'S'];
elseif lat > 0
   latS = [num2str(abs(lat)) 'N'];
else
   latS = '0';
end

if lon < 0
   lonS = [num2str(abs(lon)) 'W'];
elseif lon > 0
   lonS = [num2str(abs(lon)) 'E'];
else
   lonS = '0';
end

loc = [latS ' ' lonS ];

sdn  = nctime(infile);
good = near(sdn, t);
T = sdn(good);
S = oms_getinfo(infile);

%==========================
% First Figure: Ed, Lu, Es
%==========================
figure;
mA = subplot('position', [0.1300 0.11 .20 .65]);
a = strmatch('ed',S.symbol);
%b = find(S.instrumentDepth == 20);
b = find(S.instrumentDepth == 10);
EdCol = intersect(a, b);
lambdaEd = S.wavelength(EdCol);
for i = 1:length(EdCol)
   Ed(1,i) = infile{char(S.name(EdCol(i)))}(good);
end
Ed(find(Ed == -999)) = NaN;
plot(lambdaEd, Ed, 'b')
hold on
plot(lambdaEd, Ed, 'b.')
title('Downwelled Irradiance');
xlabel('Wavelength [nm]')
ylabel('E_d [\muW cm^-^2 nm^-^1]')
set(gca, 'YLim', [0.01 100], 'YScale', 'log', 'XLim', [300 700], 'YGrid', 'on')

% Es plot
subplot('position', [0.6922 0.11 .20 .65])
EsCol = strmatch('es',S.symbol);
lambdaEs = S.wavelength(EsCol);
h = [NaN NaN NaN];
for i = 1:length(EsCol)
   Es(1,i) = infile{char(S.name(EsCol(i)))}(good);
end
Es(find(Es == -999)) = NaN;

lambdaEf = [412 443 490 510 555 670];
Ef = frouin(T, lat, lon);
plot(lambdaEf, Ef, 'Color', [.7 .7 1]);
hold on
plot(lambdaEf, Ef, 'k.','Color', [.7 .7 1]);
h(3) = plot(lambdaEs, Es, 'g');
plot(lambdaEs, Es, 'g.')
set(gca, 'YLim', [10 1000], 'YScale', 'log', 'XLim', [300 700], 'YGrid', 'on')
title('Surface Irradiance');
xlabel('Wavelength [nm]')
ylabel('L_u [\muW cm^-^2 nm^-^1 sr^-^1]')


subplot('position', [0.4111 0.11 .20 .65])

a = strmatch('lu',S.symbol);
%b = find(S.instrumentDepth == 20);
b = find(S.instrumentDepth == 10);
Lu20Col = intersect(a, b);
lambda20 = S.wavelength(Lu20Col);
for i = 1:length(Lu20Col)
   Lu20(1,i) = infile{char(S.name(Lu20Col(i)))}(good);
end
Lu20(find(Lu20 == -999)) = NaN;

b = find(S.instrumentDepth == 1.5);
Lu1Col = intersect(a, b);
lambda1 = S.wavelength(Lu1Col);
for i = 1:length(Lu1Col)
   Lu1(1,i) = infile{char(S.name(Lu1Col(i)))}(good);
end
Lu1(find(Lu1 == -999)) = NaN;
lambdaX = unique(lambda1);
for i = 1:length(lambdaX)
   j = find(lambda1 == lambdaX(i));
   LuX(1,i) = nanmean(Lu1(j));
end
Lu1 = LuX;
lambda1 = lambdaX;

h = [NaN NaN];
h(1) = plot(lambda20, Lu20, 'b');
hold on
plot(lambda20, Lu20, 'b.')
h(2) = plot(lambda1, Lu1, 'r');
plot(lambda1, Lu1, 'r.')
set(gca, 'YLim', [0.001 10], 'YScale', 'log', 'XLim', [300 700], 'YGrid', 'on')
hold off
title('Upwelled Radiance');
xlabel('Wavelength [nm]')
ylabel('E_s [\muW cm^-^2 nm^-^1]')
lh = legend(h, '20m', '1.5m', 'Surface',0);
Pos = get(lh, 'Position');
set(lh, 'Position', [0.6922 0.8 Pos(3) Pos(4)])

set(gcf, 'CurrentAxes', mA)
text(300, 925, ['Position: ' loc]);
text(300, 400, ['Date: ' datestr(sdn(good)) ' [GMT]'])

setfontsize(7)
%=========================
% Make a gif if specified
%=========================
if nargin > 2
   [p fname ext] = fileparts(Filename);
   PaperPosition = [.25 2.5 4.5 3.25];
   set(gcf,'PaperPosition',PaperPosition);
   writegif([p filesep fname 'a' ext]);
   close(gcf)
end


%===============================
% Second Figure: K, Kchl, Lw
%===============================
figure

% First plot: Ke, Kl, Kw, Km
mA = subplot('position', [0.1300 0.11 .20 .65]);

h = [NaN NaN NaN NaN];



% Ke - Downwelled
lambdaKe = unique([lambdaEs lambdaEd]);
for i = 1:length(lambdaKe)
   a = find(lambdaEs == lambdaKe(i));
   b = find(lambdaEd == lambdaKe(i));
   if isempty(a) | isempty(b)
      Ke(1,i) = NaN;
   else
      Ke(1,i) = calc_k(Es(a), Ed(b), 20);
   end
end



% Kl - Upwelled
lambdaKl = unique([lambda1 lambda20]);
for i = 1:length(lambdaKl)
   a = find(lambda1 == lambdaKl(i));
   b = find(lambda20 == lambdaKl(i));
   if isempty(a) | isempty(b)
      Kl(1,i) = NaN;
   else
      Kl(1,i) = calc_k(Lu1(a), Lu20(b), 18.5);
   end
end

% Km - Modeled based on morel uisng [CHL] from morel K490
a = find(lambdaKe == 490);
Chl = morel_chl2(Ke(a), 490);
lambda = [412:683];
Km  = morel_k(Chl, lambda);

% Kw - Pure water based on Morel
Kw = morel_k(0, lambda);
h(4) = plot(lambda, Kw, 'b');
hold on
%plot(lambda, Kw, 'b.')
h(3) = plot(lambda, Km, 'g');
%plot(lambda, Km, 'g.')

h(1) = plot(lambdaKe, Ke, 'k');

plot(lambdaKe, Ke, 'k.')
h(2) = plot(lambdaKl, Kl, 'r');
plot(lambdaKl, Kl, 'r.')

title('Attenuation');
xlabel('Wavelength [nm]')
ylabel('K [m^-^1]')
set(gca, 'YLim', [-.2 0.6], 'XLim', [300 700],'YGrid', 'on')
lh = legend(h, 'Based on E_d','Based on L_u', 'Modeled', 'Pure Water');
Pos = get(lh, 'Position');
set(lh, 'Position', [ 0.6922 .80 Pos(3) Pos(4)])


% Second Plot Kw-K
subplot('position', [0.4111 0.11 .20 .65])
% Es
Ce = Ke - morel_k(0, lambdaKe);
plot(lambdaKe, Ce, 'k');
hold on
plot(lambdaKe, Ce, 'k.')
%Lu
Cl = Kl - morel_k(0, lambdaKl);
plot(lambdaKl, Cl, 'r');
plot(lambdaKl, Cl, 'r.')
% Modeled
Cm = Km - Kw;
plot(lambda, Cm, 'g');
%plot(lambda, Cm, 'g.')
set(gca, 'YLim', [-.1 .1], 'XLim', [300 700], 'YGrid', 'on')
title('K - K_{pure water} ')
xlabel('Wavelength [nm]')
ylabel('K [m^-^1]')

% Third plot: Lw
subplot('position', [0.6922 0.11 .20 .65])
lambda = [412 443 490 510 555];
for i = 1:length(lambdaKe)
   a = find(lambda20 == lambdaKe(i));
   if isempty(a)
      Lwe(1,i) = NaN;
   else
      Lwe(1,i) = calc_lw(Lu20(a), Ke(i), 20);
   end
end

for i = 1:length(lambdaKl)
   a = find(lambda20 == lambdaKl(i));
   Lwl(1,i) = calc_lw(Lu20(a), Kl(i), 20);
end

for i = 1:length(lambda20)
   Lwm(1,i) = calc_lw(Lu20(i), morel_k(Chl, lambda20(i)), 20);
end

plot(lambdaKe, Lwe, 'k')
hold on
plot(lambdaKe, Lwe, 'k.')
plot(lambdaKl, Lwl, 'r')
plot(lambdaKl, Lwl, 'r.')
plot(lambda20, Lwm, 'g')
plot(lambda20, Lwm, 'g.')
set(gca, 'YLim', [0 1.5], 'XLim', [300 700], 'YGrid', 'on')
title('Water Leaving Radiance')
xlabel('Wavelength [nm]')
ylabel('L_w [\muW cm^-^2 nm^-^1 sr^-^1]')

set(gcf, 'CurrentAxes', mA)
text(300, 1.3, ['Position: ' loc]);
text(300, 1.2, ['Date: ' datestr(sdn(good)) ' [GMT]'])
setfontsize(7)

%=========================
% Make a gif if specified
%=========================
if nargin > 2
   set(gcf,'PaperPosition',PaperPosition);
   writegif([p filesep fname 'b' ext]);
   close(gcf)
end
