%% Inputs
clear all; close all;
gdata = readtable("Diel_GPP_log.csv");
mission = 5;
exportpdf = false;
filename = char(gdata.Filename(mission));
load(filename, 's'); % .mat file
pdffilename = ['Mission_',num2str(gdata.MissionID(mission)),'_Profiles_',num2str(gdata.StartProf(mission)),'_',num2str(gdata.EndProf(mission)),'_Diel.pdf'];
%% QC
% Method question: Should I interp between removed values or just leave
% them out?
idiel = 156:255;
idiel([37, 92]) = [];
iupper = 1:5; % Defining surface layer 0 to approx 2.5 m
%% 
% if length(s.pH25atm) ~= length(s.pH25atm_QC) & height(s.pH25atm) ~= height(s.pH25atm_QC)
%     s.pH25atm(s.pH25atm_QC ~= 0) = [];
%     s.tc(s.tc_QC ~= 0) = [];
%     s.psal(s.psal_QC ~= 0) = [];
%     s.doxy(s.doxy_QC ~= 0) = [];
%     s.pH25atm(s.pH25atm_QC ~= 0) = [];
% end
% QC = char(strsplit(char(gdata.QC(mission)),','));
% QC = str2num(QC)';
% idiel(QC) = [];
%% Define Params
pH = mean(s.pH25atm(iupper,idiel),1,"omitmissing"); %pH 25C atm
TempC = mean(s.tc(iupper,idiel),1,"omitmissing");
Sal = mean(s.psal(iupper,idiel),1,"omitmissing");
s.o2satconc = reshape(calcO2sat(s.tc(:), s.psal(:)), size(s.pres)); % Is this the correct method for this calc?
s.o2anom = s.doxy - s.o2satconc;
meano2anom = mean(s.o2anom(iupper,idiel),1);
o2satper = mean(s.o2satper(iupper,idiel),1,"omitmissing");
% lon = s.lon(idiel);
% lat = s.lat(idiel);
depth = mean(s.depth(iupper,idiel),1,"omitmissing");
pres = mean(s.pres(iupper,idiel),1,"omitmissing");
% Coordinates = [lon; lat; depth]';
DIC = mean(s.dic_canb(iupper,idiel),1)';
sdn = s.sdn(idiel)-7/24; % GMT to local

%% Calculate anomalies from daily mean
dt = datetime(datevec(sdn));
dt.TimeZone = 'America/Los_Angeles';
myTimetable = table2timetable(table(dt,TempC', Sal',pH',meano2anom',DIC,VariableNames={'date','TempC', 'Sal','pH','meano2anom','DIC'}));
dailyAverages = retime(myTimetable, 'daily', 'mean');
dailyAnoms = table(dt,NaN(length(TempC),1),NaN(length(TempC),1),NaN(length(TempC),1),NaN(length(TempC),1),NaN(length(TempC),1),VariableNames={'date','TempC', 'Sal','pH','meano2anom','DIC'});

for i = 1:height(dailyAverages.date)
    ii = day(sdn) == day(dailyAverages.date(i)) & month(sdn) == month(dailyAverages.date(i));
    dailyAnoms.TempC(ii) = TempC(ii) - dailyAverages.TempC(i);
    dailyAnoms.Sal(ii) = Sal(ii) - dailyAverages.Sal(i);
    dailyAnoms.pH(ii) = pH(ii) - dailyAverages.pH(i);
    dailyAnoms.meano2anom(ii) = meano2anom(ii) - dailyAverages.meano2anom(i);
    dailyAnoms.DIC(ii) = DIC(ii) - dailyAverages.DIC(i);
end


%% 
plot(dailyAnoms.DIC,'.')
%% 
%% Plot data
if exportpdf == 1
    fig1 = figure('Visible','off');
else
    fig1 = figure('Visible','on');
end
tl = tiledlayout(5,1,"TileSpacing","compact");
title(tl,['Glider Mission ID: ', num2str(gdata.MissionID(mission))],'FontSize',10)
subtitle(tl,{[char(dailyAverages.date(1)), ' to ', char(dailyAverages.date(end))],...
    ['Mean Layer: ', num2str(abs(s.depth(iupper(1)))), ' m - ',...
    num2str(abs(s.depth(iupper(end)))), ' m']},'FontSize',8)

nexttile
plot(pH,'.')
ylabel('pH 25C')
xticklabels([])

nexttile
plot(TempC,'.')
ylabel('Temp (C)')
xticklabels([])

nexttile
plot(Sal,'.')
ylabel('Sal (PSU)')
xticklabels([])

nexttile
plot(meano2anom,'.') 
ylabel('O_2 Sat Anom')
xticklabels([])

nexttile
plot(DIC,'.') 
ylabel('DIC')

if exportpdf == 1
    exportgraphics(fig1,pdffilename, 'ContentType', 'vector', 'Append', true);
end
%% Plot anomalies from daily mean
if exportpdf == 1
    fig1 = figure('Visible','off');
else
    fig1 = figure('Visible','on');
end

tl = tiledlayout(5,1);
title(tl,'Anomalies (daily mean removed)','FontSize',8)

nexttile
plot(dailyAnoms.pH,'.')
ylabel('pH 25C')
xticklabels([])

nexttile
plot(dailyAnoms.TempC,'.')
ylabel('Temp (C)')
xticklabels([])

nexttile
plot(dailyAnoms.Sal,'.')
ylabel('Sal (PSU)')
xticklabels([])

nexttile
plot(dailyAnoms.meano2anom,'.')
ylabel('O_2 Sat Anom')
xticklabels([])

nexttile
plot(dailyAnoms.DIC,'.')
ylabel('DIC')

if exportpdf == 1
    exportgraphics(fig1,pdffilename, 'ContentType', 'vector', 'Append', true);
end
%% 
if exportpdf == 1
    fig1 = figure('Visible','off');
else
    fig1 = figure('Visible','on');
end

tl = tiledlayout(1,1);
title(tl,'profiles','FontSize',8)
plot(s.tc(iupper,idiel),s.pres(iupper,idiel))
hold on
plot(nanmedian(s.tc(iupper,idiel),2),s.pres(iupper,idiel),'-k','LineWidth',3)
legend