gdata = readtable("Diel_GPP_log.csv");
%%
exportpdf = 1;
for i = 1:height(gdata)
    if logical(gdata.PAR(i)) == true && logical(gdata.Quality(i)) == true
        filename = char(gdata.Filename(i));
        load(filename, 's'); % .mat file
        
        idiel = gdata.StartProf(i):gdata.EndProf(i);
        iupper = 1:5;
        s.o2satconc = reshape(calcO2sat(s.tc(:), s.psal(:)), size(s.pres));
        s.o2anom = s.doxy - s.o2satconc;
        sdn = s.sdn(idiel)-7/24; % GMT to local

        Y = mean(s.o2satper(iupper,idiel),1,'omitmissing');
        Ypar = nanmax(s.par(iupper,idiel),[],1);
        
        [diel_O2sat,SDNs] = parse_to_diel(s.sdn_(2,idiel)-7/24, Y, 24);
        diel_PAR = parse_to_diel(s.sdn_(2,idiel)-7/24, Ypar, 24);
        
        amp = nanmax(diel_O2sat,[],2) - nanmin(diel_O2sat,[],2);
        % par = nanmax(diel_PAR,[],2);
        par = nanmean(diel_PAR,2);
        % par = sum(diel_PAR,2, 'omitnan');
        par = trapz(diel_PAR,2);
        
        % max PAR
        
        X = par;
        Y = amp;
        ii = ~isnan(X) & ~isnan(Y);
        c = polyfit(X(ii), Y(ii), 1);
        
        % look at PI curves from hourly changes in O2
        
        dO2anom = diff(s.o2anom(1,idiel));
        
        dO2anom(dO2anom < -0.5 | dO2anom > 1) = NaN;
        
        PAR = s.par(1,idiel);
            PAR = PAR(1:end-1)+diff(PAR)./2;
        
        hr = hour(s.sdn_(2,idiel)-7/24); % hour in local time. 
            hr = hr(1:end-1);
        
        % get indicies for morning, or when light is increasing. 
        imorn = hr < 15; 
        
        if exportpdf == 1
            fig1 = figure('Visible','off');
        else
            fig1 = figure('Visible','on');
        end
        tl = tiledlayout(2,1);
        title(tl,['Glider Mission ID: ', num2str(gdata.MissionID(i))])
        subtitle(tl,[char(gdata.StartDate(i)), ' to ', char(gdata.EndDate(i))])
        nexttile
        hold on; box on;
        plot(X, Y, 'ko', 'markerfacecolor', 'r');
        plot(X, polyval(c,X), 'k');
        % xlabel('Daily mean PAR');
        xlabel('Daily integrated PAR');
        ylabel('O2sat amplitude (proxy for GPP)');
        title('PI Curve')
        % More sun = more production
        
        nexttile
        yyaxis left
        plot(s.sdn_(2,idiel), s.o2satper(1,idiel), '.-');
        datetick;
        ylabel('O2 %saturation');
        
        yyaxis right
        plot(s.sdn_(2,idiel), s.par(1,idiel), '.-');
        datetick;
        ylabel('PAR');
        
        if exportpdf == 1
            exportgraphics(fig1,'GliderPICurves.pdf', 'ContentType', 'vector', 'Append', true);
        end
    end
end
