exportpdf = 0;
exportpng = 1;
gdata = readtable("Diel_GPP_log_QC_Revised.csv");

mission = 7;
filename = char(gdata.Filename(mission));
load(filename, 's'); % .mat file
pdffilename = ['Mission_',num2str(gdata.MissionID(mission)),'__PIcurve.pdf'];
pngfilename = ['Mission_',num2str(gdata.MissionID(mission)),'__PIcurve.png'];
idiel = gdata.StartProf(mission):gdata.EndProf(mission);
iupper = 1:5;
s.o2satconc = reshape(calcO2sat(s.tc(:), s.psal(:)), size(s.pres));
s.o2anom = s.doxy - s.o2satconc;
sdn = s.sdn(idiel)-7/24; % GMT to local
dt = datetime(datevec(sdn));
dt.TimeZone = 'America/Los_Angeles';

tt = timetable(dt,mean(s.o2anom(iupper,idiel),1)',mean(s.o2satper(iupper,idiel))',max(s.par(iupper,idiel))','VariableNames',{'meano2anom','meano2sat','PAR'});
tthourly = retime(tt,'hourly','mean');
ttdaily = retime(tt,'daily',@trapz);
PAR_int7 = ttdaily.PAR;
ttdaily = retime(tt,'daily',@range);
AmpGPP7 = ttdaily.meano2sat;
dO2anom = diff(tthourly.meano2anom);
dO2anom(dO2anom < -0.5 | dO2anom > 1) = NaN;
PAR = tthourly.PAR(1:end-1)+diff(tthourly.PAR)./2;
hr = hour(tthourly.dt(1:end-1));
imorn = hr < 15;
PAR(~imorn) = nan;
dO2anom(~imorn) = nan;

PAR7 = PAR;
dO2anom7 = dO2anom;
PARbin7 = [0:100:400 600:200:1000 1250 1500 2000 2500];
dO2bin7 = bin_data_to_X_gpt(PAR7, dO2anom7, PARbin7);

%
mission = 8;
filename = char(gdata.Filename(mission));
load(filename, 's'); % .mat file
pdffilename = ['Mission_',num2str(gdata.MissionID(mission)),'_Profiles_164_588_PIcurve.pdf'];
idiel = gdata.StartProf(mission):gdata.EndProf(mission);
iupper = 1:5;
s.o2satconc = reshape(calcO2sat(s.tc(:), s.psal(:)), size(s.pres));
s.o2anom = s.doxy - s.o2satconc;
sdn = s.sdn(idiel)-7/24; % GMT to local
dt = datetime(datevec(sdn));
dt.TimeZone = 'America/Los_Angeles';

tt = timetable(dt,mean(s.o2anom(iupper,idiel),1)',mean(s.o2satper(iupper,idiel))',max(s.par(iupper,idiel))','VariableNames',{'meano2anom','meano2sat','PAR'});
tthourly = retime(tt,'hourly','mean');
ttdaily = retime(tt,'daily',@trapz);
PAR_int8 = ttdaily.PAR;
ttdaily = retime(tt,'daily',@range);
AmpGPP8 = ttdaily.meano2sat;
dO2anom = diff(tthourly.meano2anom);
dO2anom(dO2anom < -0.5 | dO2anom > 1) = NaN;
PAR = tthourly.PAR(1:end-1)+diff(tthourly.PAR)./2;
hr = hour(tthourly.dt(1:end-1));
imorn = hr < 15;
PAR(~imorn) = nan;
dO2anom(~imorn) = nan;

PAR8 = PAR;
dO2anom8 = dO2anom;
PARbin8 = [0:100:400 600:200:1000 1250 1500 2000 2500];
dO2bin8 = bin_data_to_X_gpt(PAR8, dO2anom8, PARbin8);


%%

if exportpdf == 1
    fig1 = figure('Visible','off');
else
    fig1 = figure('Visible','on');
end
tl = tiledlayout(1,2,"TileSpacing","tight");
title(tl,['Mission ', num2str(gdata.MissionID(mission)), ', Profiles: ', num2str(idiel(1)), ' - ', num2str(idiel(end))])
nexttile
plot(PAR7, dO2anom7, 'b.');
hold on
plot(PAR8,dO2anom8,'r.');
plot(PARbin7(2:end), dO2bin7, 'bo-', 'markerfacecolor', 'b');
plot(PARbin8(2:end), dO2bin8, 'ro-', 'markerfacecolor', 'r');

xlabel('PAR')
ylabel('dO_2/dt')
title('PI curve')
% legend() % name stations and then add a legend

nexttile
plot(PAR_int7,AmpGPP7,'o', 'markerfacecolor', 'b')
hold on
plot(PAR_int8,AmpGPP8,'o', 'markerfacecolor', 'r')
ls = lsline;
ls(1).Color = 'b';
ls(2).Color = 'r';
xlabel('Daily Integrated PAR')
ylabel('O_2 Sat Amplitude (GPP)')
title('GPP vs PAR')
if exportpdf == 1
    exportgraphics(fig1,pdffilename, 'ContentType', 'vector', 'Append', true);
elseif exportpng == 1
    exportgraphics(fig1,pngfilename, "Resolution",600);
end

