%% Inputs
clear all; close all;
gdata = readtable("Diel_GPP_log.csv");
%% 
for mission = 5

    FitModel = false;
    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 
    if isnan(gdata.StartProf(mission))
        continue
    end
    
    idiel = gdata.StartProf(mission):gdata.EndProf(mission); % QC from spreadsheet
    iupper = 1:5; % Defining surface layer 0 to approx 2.5 m
    % [Make script here for defining a mixed layer depth based on temp]
    
    % Auto QC, maybe include 4?
    if length(s.pH25atm) ~= length(s.pH25atm_QC) & height(s.pH25atm) ~= height(s.pH25atm_QC) & m == 12
        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) = [];
    else
    % Manual QC
        QC = char(strsplit(char(gdata.QC(mission)),','));
        QC = str2num(QC)';
        idiel(QC) = [];
    end

    %% 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,"omitmissing")';
    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'});
    % this can be done way easier with matrix operations.
    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
    
    %% LSF sin and cos funcs (rough estimate)
    fit_pH = LSF(dailyAnoms.date, dailyAnoms.pH);
    fit_o2 = LSF(dailyAnoms.date, dailyAnoms.meano2anom);
    fit_DIC = LSF(dailyAnoms.date, dailyAnoms.DIC);
    %% Calculate hourly diel (full time series)
    dailyAnomshourly = retime(myTimetable, 'hourly', 'mean'); % this could all be one function instead of a for loop
    ii = ~isnan(dailyAnomshourly.DIC);
    dailyAnomshourly = dailyAnomshourly(ii,:);
    for i = 1:height(dailyAverages.date)
        ii = day(sdn) == day(dailyAverages.date(i));
        dailyAnomshourly.TempC(ii) = dailyAnomshourly.TempC(ii) - dailyAverages.TempC(i);
        dailyAnomshourly.Sal(ii) = dailyAnomshourly.Sal(ii) - dailyAverages.Sal(i);
        dailyAnomshourly.pH(ii) = dailyAnomshourly.pH(ii) - dailyAverages.pH(i);
        dailyAnomshourly.meano2anom(ii) = dailyAnomshourly.meano2anom(ii) - dailyAverages.meano2anom(i);
        dailyAnomshourly.DIC(ii) = dailyAnomshourly.DIC(ii) - dailyAverages.DIC(i);
    end
    %% Calculate composite diel cycle
    dailyAnoms.hourOfDay = hour(dailyAnoms.date);
    hourlyAverage = varfun(@median, dailyAnoms, 'GroupingVariables', 'hourOfDay');
end