%% Glider Mission 21903401

%% Inputs
clear all; close all;
gdata = readtable("Diel_GPP_log.csv");
mission = 6;
exportpdf = false;
filename = char(gdata.Filename(mission));
load(filename, 's'); % .mat file

%% QC
% Method question: Should I interp between removed values or just leave
% them out?
idiel = gdata.StartProf(mission):gdata.EndProf(mission);

iupper = 1:5; % Defining surface layer 0 to approx 2.5 m

%% 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]';
Measurements = [Sal; meano2anom; TempC]';
MeasIDVec = [1, 6, 7];
TA = LIAR(Coordinates,Measurements,MeasIDVec)'; % algorithm to estimate TA
Result = CO2SYS(TA,pH,1,3,Sal,TempC,TempC,pres,pres,0, 0, 1, 10, 1);% algorithm to estimate DIC
DIC = Result(:,2); % DIC
sdn = s.sdn(idiel)-7/24; % GMT to local
vrs = mean(s.vrs(idiel,iupper),1,"omitmissing");
%% 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));
    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 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(dailyAnoms.date,pH,'.')
ylabel('pH 25C')
xticklabels([])

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

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

nexttile
plot(dailyAnoms.date,meano2anom,'.') 
ylabel('O2')
xticklabels([])

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

if exportpdf == 1
    exportgraphics(fig1,['DielAnomalies_',num2str(gdata.MissionID(mission)),'.pdf'], '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.date,dailyAnoms.pH,'.')
ylabel('pH 25C')
xticklabels([])

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

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

nexttile
plot(dailyAnoms.date,dailyAnoms.meano2anom,'.')
ylabel('O2')
xticklabels([])

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

if exportpdf == 1
    exportgraphics(fig1,['DielAnomalies_',num2str(gdata.MissionID(mission)),'.pdf'], 'ContentType', 'vector', 'Append', true);
end
%% Plot and LSF pH diel cycle
if exportpdf == 1
    fig1 = figure('Visible','off');
else
    fig1 = figure('Visible','on');
end
tl = tiledlayout(1,1);
% lsf sin and cos funcs (rough estimate)
data = dailyAnoms.pH;
time = dailyAnoms.date;
P1 = length(dailyAnoms.date)/numel(unique(day(dailyAnoms.date))); 
% If there is a diel cycle, three days should have three cycles so 
% dominant period (P1) is n_samples/#days
time = (1:1:length(time))'; % create evenly spaced time vector
A2 = [ones(length(time)) sin((2*pi*time)/P1) cos((2*pi*time)/P1)];
% 2*pi/P1 because one cycle completed every period
x2 = A2\data;
B = x2(1);
A = sqrt(sum(x2(2:end).^2));
fit_pH = A2*x2;
nexttile
plot(dailyAnoms.date,data,'.',Color='r')
hold on
plot(dailyAnoms.date,fit_pH,'k',LineWidth=1)
% plot(daily_anoms.date,movmean(daily_anoms.pH,15),Color=[.7 .7 .7],LineWidth=1,LineStyle='--') % function fit should be similar to movmean
hold off
ylim([-.2 .2])
xticks(dailyAverages.date(1):.25:dailyAverages.date(end)+1)
xtickformat('HH:mm')
xtickangle(45)
xlabel(tl,'Local time (hours:minutes)')
ylabel(tl,'pH Anomaly')
legend('pH 25C atm', 'LSF ')

if exportpdf == 1
    exportgraphics(fig1,['DielAnomalies_',num2str(gdata.MissionID(mission)),'.pdf'], 'ContentType', 'vector', 'Append', true);
end
%% Plot and LSF Mean O2 anom diel cycle
if exportpdf == 1
    fig1 = figure('Visible','off');
else
    fig1 = figure('Visible','on');
end
tl = tiledlayout(1,1);
%attempting to lsf sin and cos funcs
data = dailyAnoms.meano2anom;
time = dailyAnoms.date;
P1 = length(dailyAnoms.date)/numel(unique(day(dailyAnoms.date))); % 
% If there is a diel there should be one peak per day. This is the periodicity of the LSF.
time = (1:1:length(time))';
sf = 1;
A2 = [ones(length(time)) sf*sin((2*pi*time)/P1) sf*cos((2*pi*time)/P1)];
x2 = A2\data;
B = x2(1);
A = sqrt(sum(x2(2:end).^2));
fit_o2 = A2*x2;
nexttile
plot(dailyAnoms.date,data,'.',Color='r')
hold on
plot(dailyAnoms.date,fit_o2,'k',LineWidth=1)
% plot(dailyAnoms.date,movmean(dailyAnoms.meano2anom,15),Color=[.7 .7 .7],LineWidth=1,LineStyle='--') % function fit should be similar to movmean
hold off
ylim([-80 80])
xticks(dailyAverages.date(1):.25:dailyAverages.date(end)+1)
xtickformat('HH:mm')
xtickangle(45)
xlabel(tl,'Local time (hours:minutes)')
ylabel(tl,'O_2 Anomaly')
legend('O_2', 'LSF ')

if exportpdf == 1
    exportgraphics(fig1,['DielAnomalies_',num2str(gdata.MissionID(mission)),'.pdf'], 'ContentType', 'vector', 'Append', true);
end
%% Plot and LSF DIC diel cycle
if exportpdf == 1
    fig1 = figure('Visible','off');
else
    fig1 = figure('Visible','on');
end
tl = tiledlayout(1,1);
%attempting to lsf sin and cos funcs
data = dailyAnoms.DIC;
time = dailyAnoms.date;
P1 = length(dailyAnoms.date)/numel(unique(day(dailyAnoms.date))); % 
% If there is a diel there should be one peak per day. This is the periodicity of the LSF.
time = (1:1:length(time))';
A2 = [ones(length(time)) sin((2*pi*time)/P1) cos((2*pi*time)/P1)];
x2 = A2\data;
B = x2(1);
A = sqrt(sum(x2(2:end).^2));
fit_DIC = A2*x2;
nexttile
plot(dailyAnoms.date,data,'.',Color='r')
hold on
plot(dailyAnoms.date,fit_DIC,'k',LineWidth=1)
hold off
ylim([-120 120])
xticks(dailyAverages.date(1):.25:dailyAverages.date(end)+1)
xtickformat('HH:mm')
xtickangle(45)
xlabel(tl,'Local time (hours:minutes)')
ylabel(tl,'DIC Anomaly')
legend('DIC', 'LSF ')

if exportpdf == 1
    exportgraphics(fig1,['DielAnomalies_',num2str(gdata.MissionID(mission)),'.pdf'], 'ContentType', 'vector', 'Append', true);
end
%%
% for i = 1:height(dailyAverages.date)
%     [rs,t,d,z,a,r] = suncycle(lat(1),lon(1),dailyAverages.date);
%     rise(i) = rs(i,1) - 7;
%     set(i) = (rs(i,2) - 7) + 24;
% end
% RiseTime = dailyAverages.date;
% SetTime = dailyAverages.date;
% RiseTime.Hour = 6;
% RiseTime.Minute = 45;
% SetTime.Hour = 19;
% SetTime.Minute = 45;
% Hard coding because I can't figure this out

if exportpdf == 1
    fig1 = figure('Visible','off');
else
    fig1 = figure('Visible','on');
end

tl = tiledlayout(1,1);
nexttile
% fill([RiseTime SetTime SetTime RiseTime], [-100 -100 100 100], [.5 .5 .5], 'FaceAlpha', 0.5,'EdgeColor','none');
hold on
plot(dailyAnoms.date,dailyAnoms.meano2anom,'.',Color='r')
plot(dailyAnoms.date,fit_o2,'r',LineWidth=1)

plot(dailyAnoms.date,dailyAnoms.DIC,'.',Color='b')
plot(dailyAnoms.date,fit_DIC,'b',LineWidth=1)
plot(dailyAnoms.date,zeros(length(dailyAnoms.date)),'--',Color='k')

ylim([-100 100])
legend('Nighttime','','','', 'O_2 Fit', '', 'DIC Fit',Location='northwest')
xticks(dailyAverages.date(1):.25:dailyAverages.date(end)+1)
xtickformat('HH:mm')
xtickangle(45)
xlabel(tl,'Local time (hours:minutes)')
title('Anomaly from daily mean')
ylabel('\mu mol/kg')

% Not sure if these should be a
% seperate plot but it is really cool looking at the data without the trend
% removed
% plot(daily_anoms.date, meano2anom) 
% yyaxis right
% plot(daily_anoms.date, DIC)

if exportpdf == 1
    exportgraphics(fig1,['DielAnomalies_',num2str(gdata.MissionID(mission)),'.pdf'], 'ContentType', 'vector', 'Append', true);
end
%% Calculate GPP and CR from O2 with Barone algorithm
% dates = sdn;
% conc = dailyAnoms.meano2anom';
% lat = mean(lat);
% lon = mean(lon);
% timezone = -7;
% [rates_O2,residuals,varmat,fitline] = dielFit(dates,conc,lat,lon,timezone);
% 
% conc = dailyAnoms.DIC';
% [rates_DIC,residuals,varmat,fitline] = dielFit(dates,conc,lat,lon,timezone);
%%
% rates_O2.GPP(2)
% rates_DIC.CR(2)
% sum([rates_O2.GPP(2),rates_DIC.CR(2)])
%% Plot lat lon from data
if exportpdf == 1
    fig1 = figure('Visible','off');
else
    fig1 = figure('Visible','on');
end
geobasemap satellite
geoscatter(s.lat(idiel),s.lon(idiel))
geolimits([36.1 37.1],[-123 -121.3])
title('Location of Profiles')

if exportpdf == 1
    exportgraphics(fig1,['DielAnomalies_',num2str(gdata.MissionID(mission)),'.pdf'], 'ContentType', 'vector', 'Append', true);
end