%%
clear all;
close all;
%% Calculations 
% 1. Load in spray data.
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 = 360:587;
% 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 = trapz(diel_PAR,2);

% 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

fig1 = figure(1);
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)');

fig2 = figure(2);
nexttile;
hold on; box on;
plot(s.sdn_(2,idiel), s.o2satper(1,idiel), 'k.-');
datetick;
ylabel('O2satper');

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

fig3 = figure(3);
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');

fig4 = figure(4);
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')

fig5 = figure(5);
hold on; box on;
plot(PAR(imorn), dO2anom(imorn), 'k.');
plot(PAR(~imorn), dO2anom(~imorn), 'r.');
legend('morning','not morning')

fig6 = figure(6);
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]);
%% Create a PDF doc

saveas(fig1,'xx',)