%%
clear all;
close all;
%% Calculations 
% 1. Load in spray data.
filename = '21603401_QC.mat';
load(filename, 's'); % .mat file

% 2. Index the experimental portion of the data.
% idiel = xxx; % make clever way to index the data. maybe base this on a tbl.
idiel = 101:380;
% 3. QC if necessary.
% no need at the moment.

% 4. Calculate o2 saturation concentraion from t,s,p and reshape array
% to match the rest of the data.
s.o2satconc = reshape(calcO2sat(s.tc(:), s.psal(:)), size(s.pres));

% 5. Calculate o2 anomaly. 
	% 5a. anom = dissolved oxygen - o2 saturation concentration.
	% 5b. this normalizes o2 data because it is temp dependent.
s.o2anom = s.doxy - s.o2satconc;

% 6. Calculate mean of o2 saturation percentage in the upper 5m depth. 
% (This is indexed as 1:3 or 1:2 in barone code).

Y = mean(s.o2satper(1:3,idiel),1);

% 7. Calculate max of PAR in upper 5m.
Ypar = max(s.par(1:3,idiel),[],1,'omitmissing'); % (nanmax)

% 8. Parse_to_diel the mean of o2 saturation percentage and max of Ypar. 
% This function basically averages data into hourly sections. 
% (Can improve this function potentially).
[diel_O2sat,SDNs] = parse_to_diel(s.sdn_(2,idiel)-7/24, Y, 24);
diel_PAR = parse_to_diel(s.sdn_(2,idiel)-7/24, Ypar, 24);

% 9. Calculate the amplitude of diel o2 saturation percentage (max - min), 
% and integrate diel PAR.
amp = max(diel_O2sat,[],2,'omitmissing') - ...
    min(diel_O2sat,[],2,"omitmissing");
par = sum(diel_PAR,2,"omitmissing"); % Sum instead of trapz (lots of nans)

% 10. Fit a line to the PAR and Amplitude.
	% 10a. Barone fits 3 lines and uses all of them to estimate the diel cycle.

X = par;
Y = amp;
ii = ~isnan(X) & ~isnan(Y); % Can't calc w/ nans
c = polyfit(X(ii), Y(ii), 1);

% 11. Calculate diff O2 anom and quality control.
dO2anom = diff(s.o2anom(1,idiel));
dO2anom(dO2anom < -0.5 | dO2anom > 1) = NaN;

% 12. Calculate PAR.
PAR = s.par(1,idiel);
PAR = PAR(1:end-1)+diff(PAR)./2; % ?

% 13. Calculate hour in local time.
hr = hour(s.sdn_(2,idiel)-7/24); % hour in local time. 
hr = hr(1:end-1);

% 14. Find indices of morning hours.
imorn = hr < 15;


% 15. Define bins for PAR and bin O2 anom accordingly.
PARbin = [0:100:400 600:200:1000 1250 1500];
dO2bin = bin_data_to_X(PAR, dO2anom, PARbin);
%% Plotting
mission = char(filename);
mission = mission(1:end-4);
fig0 = figure('Visible','off');
% Add text to the center of the plot
text(.5, .5, ['Glider Mission ',mission], 'HorizontalAlignment', 'center', 'FontWeight', 'bold');
% Hide the axes
axis off;
exportgraphics(fig0, ['Glider Mission ',mission, '.pdf'], 'ContentType', 'vector');

clf
fig1 = figure('Visible','off');
hold on; box on;
plot(X, Y, 'ko', 'markerfacecolor', 'r');
plot(X, polyval(c,X), 'k');
xlabel('Daily integrated PAR');
ylabel('O2sat amplitude (proxy for GPP)');
exportgraphics(fig1, ['Glider Mission ',mission, '.pdf'], 'ContentType', 'vector', 'Append', true);

fig2 = figure('Visible','off');
nexttile;
hold on; box on;
plot(s.sdn_(2,idiel), s.o2satper(1,idiel), 'k.-');
datetick;
ylabel('O2satper');
exportgraphics(fig2, ['Glider Mission ',mission, '.pdf'], 'ContentType', 'vector', 'Append', true);

nexttile;
hold on; box on;
plot(s.sdn_(2,idiel), s.par(1,idiel), 'k.-');
datetick;
ylabel('PAR');

fig3 = figure('Visible','off');
yyaxis left
plot(s.sdn_(2,idiel), s.o2satper(1,idiel), '.-');
datetick;
ylabel('O2 %saturation');

yyaxis right
plot(s.sdn_(2,idiel), s.par(1,idiel), '.-');
datetick;
ylabel('PAR');
exportgraphics(fig3, ['Glider Mission ',mission, '.pdf'], 'ContentType', 'vector', 'Append', true);

fig4 = figure('Visible','off');
yyaxis left
plot(s.sdn_(2,idiel), s.pH25atm(1,idiel), '.-');
datetick;
ylabel('pH25C');

yyaxis right
plot(s.sdn_(2,idiel), s.o2satper(1,idiel), '.-');
ylabel('O_2Sat')
exportgraphics(fig4, ['Glider Mission ',mission, '.pdf'], 'ContentType', 'vector', 'Append', true);

fig5 = figure('Visible','off');
hold on; box on;
plot(PAR(imorn), dO2anom(imorn), 'k.');
plot(PAR(~imorn), dO2anom(~imorn), 'r.');
legend('morning','not morning')
exportgraphics(fig5, ['Glider Mission ',mission, '.pdf'], 'ContentType', 'vector', 'Append', true);

fig6 = figure('Visible','off');
hold on; box on;
plot(PAR, dO2anom, 'k.');
plot(PARbin(2:end), dO2bin, 'ko-', 'markerfacecolor', 'r');
xlabel('PAR');
ylabel('\deltaO_2');
ylim([-0.5 0.5]);
exportgraphics(fig6, ['Glider Mission ',mission, '.pdf'], 'ContentType', 'vector', 'Append', true);
%% Create a PDF doc
