%% Inputs
clear all; close all;
gdata = readtable("Diel_GPP_log_QC_Revised.csv");
%% 
% for mission = 1:height(gdata.MissionID)
for mission = 7

    FitModel = false;
    exportpdf = false;
    exportpng = 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'];
    %%
    % idiel = gdata.StartProf(mission):gdata.EndProf(mission); % QC from spreadsheet
    idiel = gdata.StartProf(mission):gdata.EndProf(mission); % before DIC starts decreasing
    iupper = 1:5; % Defining surface layer 0 to approx 2.5 m
    % iupper = 1:30;
    QC = char(strsplit(char(gdata.QC(mission)),','));
    QC = str2num(QC)';
    idiel(QC) = [];
    % %% 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]';
    sdn = s.sdn(idiel)-7/24; % GMT to local
    % DIC = mean(s.dic_canb(iupper,idiel),1,"omitmissing")'; % recalc with pHcorr
    %% Recalc DIC
    % need everything in vectors
    ta = s.ta_canb(:);
    pHin = s.pHin(:);
    tc = s.tc(:);
    psal = s.psal(:);
    pres = s.pres(:);
    sil = s.sil_canb(:);
    po4 = s.po4_canb(:); 
    
    pHcorr = s.pHcorr(:);


    % a bit of QC
    pHcorr(pHcorr > 8.8 | pHcorr < 7) = NaN;
    
    % 17 Oct 2021 - updated to CO2SYSv3 (J. Sharp) - should have Yui double
    % check the input options are correct -JSL
    % CO2SYS(PAR1,PAR2,PAR1TYPE,PAR2TYPE,SAL,TEMPIN,TEMPOUT,PRESIN,PRESOUT,SI,PO4,NH4,H2S,pHSCALEIN,K1K2CONSTANTS,KSO4CONSTANT,KFCONSTANT,BORON)
    trex = CO2SYSv3(ta, pHcorr, 1, 3, psal, tc, 20, pres, 0, sil, po4, 0, 0, 1, 10, 1, 2, 2);
    
    
    % extract variables
    DIC = reshape(trex(:,2),size(s.pres)); % CanyonB DIC
    pH20atm = reshape(trex(:,20), size(s.pres)); % changed from 18 to 20 following CO2SYSv3 update

    % s.pco2in = reshape(trex(:,4), size(s.pres));
    % s.fco2in = reshape(trex(:,5), size(s.pres));
    % s.satarin = reshape(trex(:,18), size(s.pres)); % changed from 16 to 18 following CO2SYSv3 update
    % s.satcain = reshape(trex(:,17), size(s.pres)); % changed from 15 to 17 following CO2SYSv3 update
    
    DIC = mean(DIC(iupper,idiel),1,"omitmissing")';
    pH = mean(pH20atm(iupper,idiel),1,"omitmissing"); % pH 20C atm 
    %% 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
    dailyAnoms = table2timetable(dailyAnoms);
    %% 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)
    hourlyAnoms = retime(dailyAnoms,'hourly','mean');
    %% Calculate composite diel cycle
    hourlyAnoms.hourOfDay = hour(hourlyAnoms.date);
    hourlyAverage = varfun(@nanmedian, hourlyAnoms, 'GroupingVariables', 'hourOfDay');
    hourlyAverage.date.Day = hourlyAverage.date.Day(1);
    %% Polyfit to hourly composite
    % try
    %     fDIC = fit(hourlyAverage.hourOfDay,hourlyAverage.mean_DIC,'poly3');
    %     fO2 = fit(hourlyAverage.hourOfDay,hourlyAverage.mean_meano2anom,'poly3');
    % catch
    %     sprintf('An error occured with polyfit')
    %     % continue
    % end
    %% Plot data
    if exportpdf == 1
        fig1 = figure('Visible','off');
    else
        fig1 = figure('Visible','on');
    end
    tl = tiledlayout(5,1,"TileSpacing","compact");
    title(tl,['Glider Mission ID: ', num2str(gdata.MissionID(mission))],'FontSize',10)
    subtitle(tl,{[char(dailyAverages.date(1)), ' to ', char(dailyAverages.date(end))],...
        ['Mean Layer: ', num2str(abs(s.depth(iupper(1)))), ' m - ',...
        num2str(abs(s.depth(iupper(end)))), ' m']},'FontSize',8)
    
    nexttile
    plot(dailyAnoms.date,pH,'.')
    ylabel('pH 20C')
    xticklabels([])
    
    nexttile
    plot(dailyAnoms.date,TempC,'.')
    ylabel('Temp (C)')
    xticklabels([])
    
    nexttile
    plot(dailyAnoms.date,Sal,'.')
    ylabel('Sal (PSU)')
    xticklabels([])
    
    nexttile
    plot(dailyAnoms.date,meano2anom,'.') 
    ylabel('O_2 Sat Anom')
    xticklabels([])
    
    nexttile
    plot(dailyAnoms.date,DIC,'.') 
    ylabel('DIC')
    
    if exportpdf == 1
        exportgraphics(fig1,pdffilename, 'ContentType', 'vector', 'Append', true);
    end
    
    %% Plot anomalies from daily mean
    if exportpdf == 1
        fig1 = figure('Visible','off');
    else
        fig1 = figure('Visible','on');
    end
    
    tl = tiledlayout(5,1);
    title(tl,'Anomalies (daily mean removed)','FontSize',8)
    
    nexttile
    plot(dailyAnoms.date,dailyAnoms.pH,'.')
    ylabel('pH 25C')
    xticklabels([])
    
    nexttile
    plot(dailyAnoms.date,dailyAnoms.TempC,'.')
    ylabel('Temp (C)')
    xticklabels([])
    
    nexttile
    plot(dailyAnoms.date,dailyAnoms.Sal,'.')
    ylabel('Sal (PSU)')
    xticklabels([])
    
    nexttile
    plot(dailyAnoms.date,dailyAnoms.meano2anom,'.')
    ylabel('O_2 Sat Anom')
    xticklabels([])
    
    nexttile
    plot(dailyAnoms.date,dailyAnoms.DIC,'.')
    ylabel('DIC')
    
    if exportpdf == 1
        exportgraphics(fig1,pdffilename, 'ContentType', 'vector', 'Append', true);
    end
    
    %% Plot and LSF pH diel cycle
    if exportpdf == 1
        fig1 = figure('Visible','off');
    else
        fig1 = figure('Visible','on');
    end
    tl = tiledlayout(1,1);
    
    nexttile
    plot(dailyAnoms.date,dailyAnoms.pH,'.',Color='r')
    if FitModel == 1
        hold on
        plot(dailyAnoms.date,fit_pH,'k',LineWidth=1)
        hold off
    end
    ylim auto
    if numel(dailyAverages.date(1):.25:dailyAverages.date(end)+1) < 15
        xticks(dailyAverages.date(1):.25:dailyAverages.date(end)+1)
    else
        xticks(dailyAverages.date(1):.5:dailyAverages.date(end)+1)
    end
    xtickformat('HH:mm')
    xtickangle(45)
    xlabel(tl,'Local time (hours:minutes)')
    ylabel(tl,'pH Anomaly')
    legend('pH 25C atm', 'LSF ')
    title('Anomaly from daily mean')
    if exportpdf == 1
        exportgraphics(fig1,pdffilename, 'ContentType', 'vector', 'Append', true);
    end
    %% Plot and LSF Mean O2 anom diel cycle
    if exportpdf == 1
        fig1 = figure('Visible','off');
    else
        fig1 = figure('Visible','on');
    end
    tl = tiledlayout(1,1);
    
    nexttile
    plot(dailyAnoms.date,dailyAnoms.meano2anom,'.',Color='r')
    if FitModel == 1
        hold on
        plot(dailyAnoms.date,fit_o2,'k',LineWidth=1)
        hold off
    end
    ylim auto
    if numel(dailyAverages.date(1):.25:dailyAverages.date(end)+1) < 15
        xticks(dailyAverages.date(1):.25:dailyAverages.date(end)+1)
    else
        xticks(dailyAverages.date(1):.5:dailyAverages.date(end)+1)
    end
    xtickformat('HH:mm')
    xtickangle(45)
    xlabel(tl,'Local time (hours:minutes)')
    ylabel(tl,'O_2 Sat Anom')
    legend('O_2 Sat Anom', 'LSF ')
    title('Anomaly from daily mean')
    if exportpdf == 1
        exportgraphics(fig1,pdffilename, 'ContentType', 'vector', 'Append', true);
    end
    %% Plot and LSF DIC diel cycle
    if exportpdf == 1
        fig1 = figure('Visible','off');
    else
        fig1 = figure('Visible','on');
    end
    tl = tiledlayout(1,1);
    
    nexttile
    plot(dailyAnoms.date,dailyAnoms.DIC,'.',Color='r')
    if FitModel == 1
        hold on
        plot(dailyAnoms.date,fit_DIC,'k',LineWidth=1)
        hold off
    end
    ylim auto
    if numel(dailyAverages.date(1):.25:dailyAverages.date(end)+1) < 15
        xticks(dailyAverages.date(1):.25:dailyAverages.date(end)+1)
    else
        xticks(dailyAverages.date(1):.5:dailyAverages.date(end)+1)
    end
    xtickformat('HH:mm')
    xtickangle(45)
    xlabel(tl,'Local time (hours:minutes)')
    ylabel(tl,'DIC Anomaly')
    legend('DIC', 'LSF ')
    title('Anomaly from daily mean')
    if exportpdf == 1
        exportgraphics(fig1,pdffilename, 'ContentType', 'vector', 'Append', true);
    end
    %% Plot O2 and DIC together
    % Need function for nighttime fill
    
    if exportpdf == 1
        fig1 = figure('Visible','off');
    else
        fig1 = figure('Visible','on');
    end
    
    tl = tiledlayout(1,1);
    nexttile
    % fill([RiseTime SetTime SetTime RiseTime], [-100 -100 100 100], [.5 .5 .5], 'FaceAlpha', 0.5,'EdgeColor','none');
    % fill(t, y, [.5 .5 .5], 'FaceAlpha', 0.5,'EdgeColor','none');
    
    plot(dailyAnoms.date,dailyAnoms.meano2anom,'.',Color='r')
    hold on
    plot(dailyAnoms.date,dailyAnoms.DIC,'.',Color='b')
    plot(dailyAnoms.date,zeros(length(dailyAnoms.date)),'--',Color='k')
    if FitModel == 1
        plot(dailyAnoms.date,fit_o2,'r',LineWidth=1)
        plot(dailyAnoms.date,fit_DIC,'b',LineWidth=1)
    end
    legend('O_2 Anom','DIC',Location='northwest')
    ylim auto
    
    if numel(dailyAverages.date(1):.25:dailyAverages.date(end)+1) < 15
        xticks(dailyAverages.date(1):.25:dailyAverages.date(end)+1)
    else
        xticks(dailyAverages.date(1):.5:dailyAverages.date(end)+1)
    end
    xtickformat('HH:mm')
    xtickangle(45)
    xlabel(tl,'Local time (hours:minutes)')
    title('Anomaly from daily mean')
    ylabel('\mumol/kg')
    
    if exportpdf == 1
        exportgraphics(fig1,pdffilename, 'ContentType', 'vector', 'Append', true);
    end
    %% Plot hourly averaged O2 and DIC together
    % Need function for nighttime fill
    
    if exportpdf == 1
        fig1 = figure('Visible','off');
    else
        fig1 = figure('Visible','on');
    end

    tl = tiledlayout(1,1);
    nexttile
    % fill([RiseTime SetTime SetTime RiseTime], [-100 -100 100 100], [.5 .5 .5], 'FaceAlpha', 0.5,'EdgeColor','none');
    % fill(t, y, [.5 .5 .5], 'FaceAlpha', 0.5,'EdgeColor','none');

    plot(hourlyAnoms.date,hourlyAnoms.meano2anom,'.',Color='r')
    hold on
    plot(hourlyAnoms.date,hourlyAnoms.DIC,'.',Color='b')
    plot(hourlyAnoms.date,zeros(length(hourlyAnoms.date)),'--',Color='k')
    % if FitModel == 1
    %     plot(dailyAnomshourly.date,fit_o2,'r',LineWidth=1)
    %     plot(dailyAnomshourly.date,fit_DIC,'b',LineWidth=1)
    % end
    legend('O_2 Anom','DIC',Location='northwest')
    ylim auto

    if numel(dailyAverages.date(1):.25:dailyAverages.date(end)+1) < 15
        xticks(dailyAverages.date(1):.25:dailyAverages.date(end)+1)
    else
        xticks(dailyAverages.date(1):.5:dailyAverages.date(end)+1)
    end
    xtickformat('HH:mm')
    xtickangle(45)
    xlabel(tl,'Local time (hours:minutes)')
    title('Anomaly from daily mean, hourly avg')
    ylabel('\mumol/kg')

    if exportpdf == 1
        exportgraphics(fig1,pdffilename, 'ContentType', 'vector', 'Append', true);
    end
    %% Composite Diel
    %fix axes etc... this plot should be same as above except the composite for
    %the three day period
    
    if exportpdf == 1
        fig1 = figure('Visible','off');
    else
        fig1 = figure('Visible','on');
    end
    
    tl = tiledlayout(1,1);
    nexttile
    % fill([0 6.75 6.75 0],[-100 -100 100 100],[.5 .5 .5],'FaceAlpha', 0.5,'EdgeColor','none')
    % hold on
    % fill([19.75 23 23 19.75],[-100 -100 100 100],[.5 .5 .5],'FaceAlpha', 0.5,'EdgeColor','none')
    
    plot(hourlyAverage.hourOfDay,hourlyAverage.nanmedian_meano2anom,'-o',Color='r')
    hold on
    % plot(hourlyAverage.hourOfDay,fO2(hourlyAverage.hourOfDay),Color='r',LineStyle='--')
    plot(hourlyAverage.hourOfDay,hourlyAverage.nanmedian_DIC,'-o',Color='b')
    % plot(hourlyAverage.hourOfDay,fDIC(hourlyAverage.hourOfDay),Color='b',LineStyle='--')
    
    plot(hourlyAverage.hourOfDay,zeros(length(hourlyAverage.hourOfDay)),'--',Color='k')

    lgd = legend('O_2 Anom Composite', '','DIC Anom Composite','','',Location='northwest');
    lgd.FontSize = 5;
    xticks([0, 6, 12, 18, 23])
    xticklabels({'00:00','06:00','12:00','18:00','00:00'})
    xtickangle(45)
    xlabel(tl,'Local time (hours:minutes)')
    title([num2str(numel(dailyAverages.date)), ' day mean diel cycle'])
    subtitle('medians, daily mean removed')
    ylabel('\mumol/kg')
    xlim([0 23])
    ylim auto
    if exportpdf == 1
        exportgraphics(fig1,pdffilename, 'ContentType', 'vector', 'Append', true);
    end
    %% Calculate GPP and CR from O2 with Barone algorithm
    % dates = sdn;
    % lat = mean(lat);
    % lon = mean(lon);
    % timezone = -7;
    % 
    % conc = meano2anom;
    % [rates_O2,fig1] = dielFit(dates,conc,lat,lon,timezone);
    % title('O_2 Anom')
    % subtitle(['Sin fit, GPP: ', num2str(rates_O2.GPP(2)), ', CR: ', num2str(rates_O2.CR(2))])
    % if exportpdf == 1
    %     exportgraphics(fig1,pdffilename, 'ContentType', 'vector', 'Append', true);
    % end
    % conc = DIC';
    % [rates_DIC,fig1] = dielFit(dates,conc,lat,lon,timezone);
    % title('DIC')
    % subtitle(['Sin fit, GPP: ', num2str(rates_DIC.GPP(2)), ', CR: ', num2str(rates_DIC.CR(2))])
    % if exportpdf == 1
    %     exportgraphics(fig1,pdffilename, 'ContentType', 'vector', 'Append', true);
    % end
    % conc = dailyAnoms.meano2anom';
    % [rates_O2_dailyanom,fig1] = dielFit(dates,conc,lat,lon,timezone);
    % title('O_2 Anom daily mean removed')
    % subtitle(['Sin fit, GPP: ', num2str(rates_O2_dailyanom.GPP(2)), ', CR: ', num2str(rates_O2_dailyanom.CR(2))])
    % if exportpdf == 1
    %     exportgraphics(fig1,pdffilename, 'ContentType', 'vector', 'Append', true);
    % end
    % conc = dailyAnoms.DIC';
    % [rates_DIC_dailyanom,fig1] = dielFit(dates,conc,lat,lon,timezone);
    % title('DIC daily mean removed')
    % subtitle(['Sin fit, GPP: ', num2str(-rates_DIC_dailyanom.GPP(2)), ', CR: ', num2str(-rates_DIC_dailyanom.CR(2))])
    % 
    % if exportpdf == 1
    %     exportgraphics(fig1,pdffilename, 'ContentType', 'vector', 'Append', true);
    % end
    %%
    if exportpdf == 0
        fig1 = figure('Visible','on');
    end
    conc = hourlyAverage.nanmedian_meano2anom';
    dates = datenum(hourlyAverage.date)';
    lat = mean(lat);
    lon = mean(lon);
    timezone = -7;
    [rates_O2,fig1] = dielFit(dates,conc,lat,lon,timezone);
    sc = plot(hourlyAnoms.hourOfDay/24,hourlyAnoms.meano2anom,'o',Color=[.7 .7 .7]);
    uistack(sc,"bottom");
    % title('O_2 Anomaly, daily mean removed, median composite')
    subtitle(['Linear fit, GOP: ', num2str(abs(rates_O2.GPP(2))), ', CRO: ', num2str(abs(rates_O2.CR(2)))])
    lg = legend('Hourly','Median','Sin Fit');
    set(lg,'Fontsize',8)
    ylabel('O_2 anom \mumol/kg')
    if exportpdf == 1
        exportgraphics(fig1,pdffilename, 'ContentType', 'vector', 'Append', true);
    elseif exportpng == 1
        pngfilename = ['Mission_',num2str(gdata.MissionID(mission)),'_Profiles_',num2str(gdata.StartProf(mission)),'_',num2str(gdata.EndProf(mission)),'_Diel_O2.png'];
        exportgraphics(fig1,pngfilename, "Resolution",600);
    end
    
    if exportpdf == 1
        fig1 = figure('Visible','off');
    else
        fig1 = figure('Visible','on');
    end
    conc = hourlyAverage.nanmedian_DIC';
    [rates_DIC,fig1] = dielFit(dates,conc,lat,lon,timezone);
    sc = plot(hourlyAnoms.hourOfDay/24,hourlyAnoms.DIC,'o',Color=[.7 .7 .7]);
    uistack(sc,"bottom");
    subtitle(['Sin fit, GCP: ', num2str(abs(rates_DIC.GPP(2))), ', CRC: ', num2str(abs(rates_DIC.CR(2)))])
    lg = legend('Hourly','Median','Sin Fit');
    set(lg,'Fontsize',8)
    ylabel('DIC anom \mumol/kg')
    if exportpdf == 1
        exportgraphics(fig1,pdffilename, 'ContentType', 'vector', 'Append', true);
    elseif exportpng == 1
        pngfilename = ['Mission_',num2str(gdata.MissionID(mission)),'_Profiles_',num2str(gdata.StartProf(mission)),'_',num2str(gdata.EndProf(mission)),'_Diel_DIC.png'];
        exportgraphics(fig1,pngfilename, "Resolution",600);
    end
   
    %% Redfield ratio
    % redfield = rates_O2.GPP(2)/rates_DIC.GPP(2); % O2:CO2 | 138:106
    %sum([rates_O2.GPP(2),rates_DIC.CR(2)])
    
    %% Plot lat lon from data
    if exportpdf == 1
        fig1 = figure('Visible','off');
    else
        fig1 = figure('Visible','on');
    end
    
    geobasemap satellite
    geoscatter(s.lat(idiel),s.lon(idiel))
    mosslat = 36.80594819607352;
    mosslon = -121.786377002975;
    ellipsoid = wgs84Ellipsoid('miles');
    dist = abs(distance(mosslat, mosslon, lat, lon,ellipsoid));
    if dist < 30
        n = 1;
        geolimits([lat-n lat+n],[lon-n lon+n])
    else
        n = 2;
        geolimits([lat-n lat+n],[lon-n lon+n])
    end
    title('Location of Profiles')
    subtitle(['Distance from Moss Landing: ', num2str(dist),' Miles'])
    % txt = {['GPP: ', num2str(rates_O2_dailyanom.GPP(2))],['CR: ', num2str(rates_DIC_dailyanom.GPP(2))],['Redfield: ', num2str(redfield)]};
    % text(lat+.5,lon-.5,txt)
    if exportpdf == 1
        exportgraphics(fig1,pdffilename, 'ContentType', 'vector', 'Append', true);
    end
end