%% Glider Mission 21202901

% reproduce methods from 21603401. Need a consistent method for QCing
% profiles
clear all; close all;
gdata = readtable("Diel_GPP_log.csv");
mission = 1;

filename = char(gdata.Filename(mission));
load(filename, 's'); % .mat file

idiel = gdata.StartProf(mission):gdata.StartProf(mission)+108; % sporatic after profile 227
% idiel([123,127,130,134,154,160]) = [];
idiel([13,37,47,59,74,109]) = [];

iupper = 1:5;
pH = mean(s.pH25atm(iupper,idiel),1,"omitmissing");
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);
sdn = s.sdn(idiel)-7/24; % GMT to local


[diel_pH, SDNs] = parse_to_diel(sdn, pH, 24);
[diel_TempC, ~] = parse_to_diel(sdn, TempC, 24);
[diel_Sal, ~] = parse_to_diel(sdn, Sal, 24);
[diel_meano2anom, ~] = parse_to_diel(sdn, meano2anom, 24);
LegendNames = cellstr(datestr(SDNs, 'mm/dd/yy'));

fig1 = figure('Visible','on');
tl = tiledlayout(4,1);
title(tl,['Glider Mission ID: ', num2str(gdata.MissionID(mission))])
subtitle(tl,{['pH Diel Cycle ', char(gdata.StartDate(mission)), ' to ', char(gdata.EndDate(mission))],...
    ['Mean Layer: ', num2str(abs(s.depth(iupper(1)))), ' m - ', num2str(abs(s.depth(iupper(end)))), ' m']})

nexttile
plot(pH,'.')
ylabel('pH 25C')
title('pH full time series')

nexttile
plot(TempC,'.')
ylabel('Temp (C)')
title('Temperature full time series')

nexttile
plot(Sal,'.')
ylabel('Sal (PSU)')
title('Salinity full time series')

nexttile
plot(meano2anom,'.') 
ylabel('O2')
title('mean o2 anom full time series')

%exportgraphics(fig1,['DielAnomalies_',num2str(gdata.MissionID(mission)),'.pdf'], 'ContentType', 'vector', 'Append', true);
%% Diels

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

for i = 1:height(dailyAverages.date)
    ii = day(sdn) == day(dailyAverages.date(i));
    daily_anoms.TempC(ii) = TempC(ii) - dailyAverages.TempC(i);
    daily_anoms.Sal(ii) = Sal(ii) - dailyAverages.Sal(i);
    daily_anoms.pH(ii) = pH(ii) - dailyAverages.pH(i);
    daily_anoms.meano2anom(ii) = meano2anom(ii) - dailyAverages.meano2anom(i);
end
%% 
fig1 = figure('Visible','on');
clf
tl = tiledlayout(4,1);
title(tl,['Glider Mission ID: ', num2str(gdata.MissionID(mission))])
subtitle(tl,{'Diel Anomalies (daily mean removed)',...
    ['Mean Layer: ', num2str(abs(s.depth(iupper(1)))), ' m - ', num2str(abs(s.depth(iupper(end)))), ' m']})

nexttile
plot(daily_anoms.pH,'.')
ylabel('pH 25C')
title('pH full time series')

nexttile
plot(daily_anoms.TempC,'.')
ylabel('Temp (C)')
title('Temperature full time series')

nexttile
plot(daily_anoms.Sal,'.')
ylabel('Sal (PSU)')
title('Salinity full time series')

nexttile
plot(daily_anoms.meano2anom,'.') 
ylabel('O2')
title('mean o2 anom full time series')

%exportgraphics(fig1,['DielAnomalies_',num2str(gdata.MissionID(mission)),'.pdf'], 'ContentType', 'vector', 'Append', true);
%%
fig1 = figure('Visible','on');

tl = tiledlayout(1,1);
plot(daily_anoms.date,daily_anoms.meano2anom,'.')
hold on
yyaxis right
plot(daily_anoms.date,daily_anoms.pH,'.')
legend('mean o_2 anom','pH_25 atm',Location='northwest')
title(tl,['Glider Mission ID: ', num2str(gdata.MissionID(mission)), ', diel cycle anomalies (daily mean removed)'])
subtitle(tl,['Mean Layer: ', num2str(abs(s.depth(iupper(1)))), ' m - ', num2str(abs(s.depth(iupper(end)))), ' m'])
%exportgraphics(fig1,['DielAnomalies_',num2str(gdata.MissionID(mission)),'.pdf'], 'ContentType', 'vector', 'Append', true);
%%
fig1 = figure('Visible','on');

%attempting to lsf sin and cos funcs
data = daily_anoms.pH;
time = daily_anoms.date;
P1 = 225/6; % period of cycle. 225 time steps and 6 days. If diel then 1 peak per day.
P2 = 225/12; % this is next highest frq from the power spectrum

time = (1:1:length(time))';
A2 = [ones(length(time)) sin((2*pi*time)/P1) cos((2*pi*time)/P1) sin((2*pi*time)/P2) cos((2*pi*time)/P2)];
x2 = A2\data;
B = x2(1);
A = sqrt(sum(x2(2:end).^2));
fit = A2*x2;

plot(time,daily_anoms.pH,'.')
hold on
plot(time,fit)
hold off
%exportgraphics(fig1,['DielAnomalies_',num2str(gdata.MissionID(mission)),'.pdf'], 'ContentType', 'vector', 'Append', true);
%% 
fig1 = figure('Visible','off');

geoscatter(s.lat(idiel),s.lon(idiel))
geolimits auto
exportgraphics(fig1,['DielAnomalies_',num2str(gdata.MissionID(mission)),'.pdf'], 'ContentType', 'vector', 'Append', true);
%%
%conc and dates need to match up to perform the fit
% dates = sdn;
% conc = detrend(meano2anom);
% lat = mean(s.lat(idiel));
% lon = mean(s.lon(idiel));
% timezone = -7;
% [rates,residuals,varmat,fitline] = dielFit(dates,conc,lat,lon,timezone);
%yay now I can calc GPP!