DepLog = ReadDeploymentLog();
[missionIDs, idxunique] = unique(DepLog.("Deployment ID"));
missionIDs = cellstr(missionIDs);
depthnegative = NaN(length(missionIDs),1);
AdjustpHCorr = NaN(length(missionIDs),1);
mn_adj = NaN(length(missionIDs),1);
mn_noadj = NaN(length(missionIDs),1);
Cal_file_deltapHfit = NaN(length(missionIDs),1);
filetype = cell(length(missionIDs),1);
filetype(:) = {NaN};
depsum = DepLog.("Deployment summary")(idxunique);
for i = 1:length(missionIDs)
    %% Load 's' and 'Cal' from a deployment
    basepath = '\\atlas.shore.mbari.org\ProjectLibrary\901805_Coastal_Biogeochemical_Sensing\Spray_Data';
    % try to read in the QC file. If that fails, try to read in the
    % calibrated file. If that fails just continue
    try
        fileName = fullfile(basepath,char(missionIDs(i)),'QC',[char(missionIDs(i)),'_QC.mat']);
        load(fileName);
        filetype{i} = 'QC';
        Cal_file_deltapHfit(i) = isfield(cal,'deltapHfit') | isfield(Cal,'deltapHfit') | isfield(s,'deltapHfit'); % Does the delta pH var exist?
    catch
        try
            fileName = fullfile(basepath,char(missionIDs(i)),'QC',[char(missionIDs(i)),'_calibrated.mat']);
            load(fileName);
            filetype{i} = 'Calibrated';
            Cal_file_deltapHfit(i) = isfield(cal,'deltapHfit') | isfield(Cal,'deltapHfit');
        catch
            continue
        end
    end
    
    %%
    if mean(s.depth,'all','omitmissing') < 0
        depthnegative(i) = true;
        s.depth = -s.depth;
    elseif mean(s.depth,'all','omitmissing') > 0
        depthnegative(i) = false;
    end
    %% Recalculate pHCorr using the corrected equation
    try
        pHinCorrAdj = (-Cal.deltapHfit.*s.Tadj) + s.pHin; % flipped sign of Cal.deltapHfit
    catch
        try
            pHinCorrAdj = (-s.deltapHfit.*s.Tadj) + s.pHin; % flipped sign of Cal.deltapHfit
        catch
            AdjustpHCorr(i) = 1;
            continue
        end
    end
    %% Create variables
    time = s.sdn;
    pHin = s.pHcorr; % Drift adjustment applied
    pHin(pHin == -999) = nan;
    % pHin(s.pHin_QC ~= 1) = nan; % 1 means good data
    
    pHin_CanB = s.pHin_canb;
    pHin_CanB(pHin_CanB == -999) = nan;
    % pHin_CanB(s.pHin_QC ~= 1) = nan;
    
    pHcanB_residuals = pHin - pHin_CanB;
    pHcanB_residuals_adj = pHinCorrAdj - pHin_CanB;
    
    ncolorblocks = 60;
    ncolorblocksredblue = 30; % cbar range is 0.3. Each block is .01. White space is .02.
    
    depthidx = s.depth >= 425 & s.depth <= 450;
    mn_noadj(i) = mean(pHcanB_residuals(depthidx),'all','omitmissing');
    mn_adj(i) = mean(pHcanB_residuals_adj(depthidx),'all','omitmissing');
    
    % Threshold
    if abs(mn_noadj(i)) > 0.001
        AdjustpHCorr(i) = true;
    elseif abs(mn_noadj(i)) <= 0.001
        AdjustpHCorr(i) = false;
    end
    
    %% Plotting: Residuals vs depth
    
    cla; clf; close('all')
    fig = figure;
    tl = tiledlayout(1,2);
    title(tl,{['pHin Corrected - pHin CanyonB'],['Adjustment is flipping sign of pH correction'] ,['MissionID: ',char(missionIDs(i))]});
    if depthnegative(i) == true
        sbtxt = 'Depth is negative';
    elseif depthnegative(i) == false
        sbtxt = 'Depth is positive';
    end
    subtitle(tl,sbtxt)
    ylabel(tl,'Depth (m)')
    xlabel(tl,'pH Residuals')
    
    nexttile
    hold on;
    histogram2(pHcanB_residuals,s.depth,"Normalization","pdf","DisplayStyle","tile",'FaceColor','flat','BinWidth',[.01 10])
    set(gca,'Ydir','reverse')
    colorbar
    grid on
    title('Pre-adjustment')
    xline(0,'r-','LineWidth',2)
    yline(450,'r-','LineWidth',2)
    shading("interp")
    clim([0 .04])
    scatter(0,450,'o','red')
    scatter(mn_noadj(i),450,'o','green')
    
    nexttile
    hold on;
    histogram2(pHcanB_residuals_adj,s.depth,"Normalization","pdf","DisplayStyle","tile",'FaceColor','flat','BinWidth',[.01 10])
    set(gca,'Ydir','reverse')
    colorbar
    grid on
    title('Post-adjustment')
    xline(0,'r-','LineWidth',2)
    yline(450,'r-','LineWidth',2)
    shading("interp")
    clim([0 .04])
    scatter(0,450,'o','red')
    scatter(mn_adj(i),450,'o','green')
    
    % Link axes for same x and y limits
    ax = findall(fig, 'Type', 'axes');
    linkaxes(ax, 'xy');
    xMax = max(abs([xlim(ax(1))]));
    xlim(ax, [-xMax, xMax]);

    fsavename = fullfile('C:','Users','bwerb','Documents','SprayDataCorrections','Figures',['missionID_',char(missionIDs(i))]);
    set(gcf, 'PaperOrientation', 'landscape');
    set(gcf, 'PaperUnits', 'inches');
    set(gcf, 'PaperSize', [11 8.5]); % Standard letter size in landscape
    set(gcf, 'PaperPosition', [0 0 11 8.5]); % Fill the whole paper
    saveas(fig,fsavename,'png')
    disp(['Finished ',num2str(i), ' / ', num2str(length(missionIDs))])
    close('all')
end
%%
Info = table(missionIDs, AdjustpHCorr, depthnegative, mn_noadj, mn_adj, filetype,Cal_file_deltapHfit, depsum,...
          'VariableNames', {'MissionID', 'Need_pHCorr_Adjustment', 'Depth_negative', 'mean_residual450_noadj','mean_residual450_adj','filetype','Cal_file_deltapHfit','Deployment_Summary'});
% Save the table as a CSV file
writetable(Info,'Info.csv');