gdata = readtable("Diel_GPP_log.csv");
delete('GliderpHAveraging.pdf')
% mndiel = NaN(height(gdata.MissionID),24);
for i = 1:height(gdata.MissionID)
    filename = char(gdata.Filename(i));
    load(filename, 's'); % .mat file
    
    idiel = gdata.StartProf(i):gdata.EndProf(i); % Chosen from ODV
    iupper = 1:5; % depth range to avg
    pH = mean(s.pH25atm(iupper,idiel),1,"omitmissing");
    sdn = s.sdn(idiel)-7/24;
    
    [diel, SDNs] = parse_to_diel(sdn, pH, 24);
    LegendNames = cellstr(datestr(SDNs, 'mm/dd/yy'));
    
    fig1 = figure('Visible','off');
    tl = tiledlayout(4,1,"TileSpacing","compact");
    title(tl,['Glider Mission ID: ', num2str(gdata.MissionID(i))])
    subtitle(tl,{['pH Diel Cycle ', char(gdata.StartDate(i)), ' to ', char(gdata.EndDate(i))],...
        ['Mean Layer: ', num2str(abs(s.depth(iupper(1)))), ' m - ', num2str(abs(s.depth(iupper(end)))), ' m']})
    nexttile
    plot(sdn,pH,'.')
    datetick('x',6,'keeplimits')
    ylabel('pH 25C')
    % xlabel('date')
    title('pH full time series')
    
    diel_detrend = detrend(diel,1,'omitmissing');
    mn = mean(diel_detrend,2,'omitmissing'); % maybe movmean
    sd = std(diel_detrend,0,2,'omitmissing');
    iweird = diel_detrend >= mn + 2*sd | diel_detrend <= mn - 2*sd; % simple outlier filter
    diel_detrend(iweird) = NaN;
    % mndiel(i,:) = mn(diel_detrend,1,"omitmissing");
    nexttile
    plot(reshape(diel_detrend.',1,[]),'.')
    title('Each Day Detrended with mean +/- 2 stds Removed')
    ylabel('Anomaly')

    nexttile
    plot(diel_detrend','.')

    nexttile
    plot(range(diel_detrend,1),'.')
    title('Range')
    ylabel('Anomaly')
    % 
    % nexttile
    % plot(mndiel(i,:),'.')
    % title('Mean')
    % ylabel('Anomaly')
    
    
    exportgraphics(fig1,'GliderpHAveraging.pdf', 'ContentType', 'vector', 'Append', true);

end
