clear all
close all

%

% s = load_spray_sat('22602901');

% load('\\atlas\projectlibrary\901805_Coastal_Biogeochemical_Sensing\Spray_Data\22602901\sat_files\22602901_sat.mat', 's');

load('22602901_QC.mat', 's');

%

% index for diel cycle expeirment #1
%idiel = 127:323; % 127:323 are 150 m dives
% 
% index for diel cycle experiment #2
idiel = 360:587; 

%

% a bit of QC; dont' need this for full res data! 
% s.PAR(s.PAR>5000) = NaN; 

%

% 100% sat concentration
s.o2satconc = reshape(calcO2sat(s.tc(:), s.psal(:)), size(s.pres));
% oxygen anomaly; [O2] - [O2]_100%
s.o2anom = s.doxy - s.o2satconc;

% Y= mean(s.o2satper(1:3,:),1);
Y = mean(s.o2satper(1:3,idiel),1);
Ypar = nanmax(s.par(1:3,idiel),[],1);

%

[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);

%

amp = nanmax(diel_O2sat,[],2) - nanmin(diel_O2sat,[],2);
% par = nanmax(diel_PAR,[],2);
par = nanmean(diel_PAR,2);
% par = sum(diel_PAR,2, 'omitnan');
par = trapz(diel_PAR,2);

% max PAR

X = par;
Y = amp;
ii = ~isnan(X) & ~isnan(Y);
c = polyfit(X(ii), Y(ii), 1);

%

% look at PI curves from hourly changes in O2

dO2anom = diff(s.o2anom(1,idiel));

dO2anom(dO2anom < -0.5 | dO2anom > 1) = NaN;

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

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

% get indicies for morning, or when light is increasing. 
imorn = hr < 15; 

%

clf
figure(1);
hold on; box on;
plot(X, Y, 'ko', 'markerfacecolor', 'r');
plot(X, polyval(c,X), 'k');
% xlabel('Daily mean PAR');
xlabel('Daily integrated PAR');
ylabel('O2sat amplitude (proxy for GPP)');
% More sun = more production

% saveas(f0, 'diel_O2amp_vs_mean_PAR.jpg');


%



%

clf
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');
% saveas(f1, 'O2cycle_and_PAR.jpg');

%

clf
figure(4)

yyaxis left
    
plot(s.sdn_(2,idiel), detrend(s.pH25atm(1,idiel)), '.-');
datetick;
ylabel('pH25C');

yyaxis right
plot(s.sdn_(2,idiel), s.o2satper(1,idiel), '.-');
ylabel('O_2Sat')

%

% Why is there a trend in pH?
clf
figure(5)
hold on; box on;
plot(PAR(imorn), dO2anom(imorn), 'k.');
plot(PAR(~imorn), dO2anom(~imorn), 'r.');
legend('morning','not morning')

%

% thanks gpt
morning_data = PAR(imorn);
not_morning_data = PAR(~imorn);

% Perform t-test
[h, p] = ttest2(morning_data, not_morning_data);

% Display result
if h == 1
    disp('There is a significant difference between morning and not-morning PAR.');
else
    disp('There is no significant difference between morning and not-morning PAR.');
end

disp(['p-value: ', num2str(p)]);

clf
histogram(morning_data,Normalization="pdf");
hold on
histogram(not_morning_data,Normalization="pdf");

%

morning_data = dO2anom(imorn);
not_morning_data = dO2anom(~imorn);

% Perform t-test
[h, p] = ttest2(morning_data, not_morning_data);

% Display result
if h == 1
    disp('There is a significant difference between morning and not-morning dO2.');
else
    disp('There is no significant difference between morning and not-morning dO2.');
end

disp(['p-value: ', num2str(p)]);
clf
histogram(morning_data,Normalization="pdf");
hold on
histogram(not_morning_data,Normalization="pdf");
%Interesting result. There is a sig diff between morn and non morn PAR but
%not true for dO2. Something else is involved.

%

% Creating bins for PAR and binning O2 to match. 
PARbin = [0:100:400 600:200:1000 1250 1500];
dO2bin = bin_data_to_X(PAR, dO2anom, PARbin);
clf
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]);

%

% Normal
clf
histogram(dO2anom,Normalization="pdf")