exportpdf = 0;
gdata = readtable("Diel_GPP_log.csv");
for mission = 1:height(gdata)
    if logical(gdata.PAR(mission)) == true && logical(gdata.Quality(mission)) == true
        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)),'_PIcurve.pdf'];
        idiel = gdata.StartProf(mission):gdata.EndProf(mission);
        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
        dt = datetime(datevec(sdn));
        dt.TimeZone = 'America/Los_Angeles';

        tt = timetable(dt,mean(s.o2anom(iupper,idiel),1)',max(s.par(iupper,idiel))','VariableNames',{'meano2anom','PAR'});
        tthourly = retime(tt,'hourly','mean');
        ttdaily = retime(tt,'daily','mean');
        dO2anom = diff(tthourly.meano2anom);
        dO2anom(dO2anom < -0.5 | dO2anom > 1) = NaN;
        PAR = tthourly.PAR(1:end-1)+diff(tthourly.PAR)./2;
        hr = hour(tthourly.dt(1:end-1));
        imorn = hr < 15;

        if exportpdf == 1
            fig1 = figure('Visible','off');
        else
            fig1 = figure('Visible','on');
        end
        % hist3(X,'nbins',[25 25],'CdataMode','auto')
        % view(2)
        % shading interp
        % grid off
        % colormap(flipud(hot(25)))
        % colorbar

        scatter(PAR(imorn),dO2anom(imorn))
        xlabel('PAR')
        ylabel('dO_2/dt')
        title('PI curve hourly average')
        % subtitle([char(ttdaily.dt(1)), ' to ', char(ttdaily.dt(end))])
        subtitle('June 22 to July 13')
        hold on


        if exportpdf == 1
            exportgraphics(fig1,pdffilename, 'ContentType', 'vector', 'Append', true);
        end
    end
end

