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'});
% daily_anom = NaN(1,length(TempC));

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

%Removing daily mean seems like a very promising approach!!
% pick a better way to do this
%%
%%
figure(1)
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',Location='northwest')
title(tl,['Glider Mission ID: ', num2str(gdata.MissionID(i)), ', diel cycle anomalies (daily mean removed)'])
subtitle(tl,['Mean Layer: ', num2str(abs(s.depth(iupper(1)))), ' m - ', num2str(abs(s.depth(iupper(end)))), ' m'])


%%
[fit] = LSF(daily_anoms.pH,daily_anoms.date, 12/24);
%%

%%
