%% pH_sections_3panel.m  –  Three stacked pH section plots

%% ── Inputs ───────────────────────────────────────────────────────────────
data_files = { ...
    "C:\Users\bwerb\Documents\ProposalFigures\data\25B68601sat.mat"; ...
    "C:\Users\bwerb\Documents\ProposalFigures\data\26220301sat.mat"; ...
    "C:\Users\bwerb\Documents\ProposalFigures\data\26120501sat.mat" };
panel_titles = {'Trinidad Head, OR', 'Monterey Bay, CA', 'San Diego, CA'};
pH_field   = 'pHin';
date_start = datenum(2025, 8, 1);
date_end   = datenum(2026, 4, 14);
pH_clim    = [7.5  8.05];
bin_size   = 50;
gap_thresh = 2;    % days — gaps longer than this render as white

%% ── Load, filter, and bin each glider ───────────────────────────────────
n = numel(data_files);
sdn_all = cell(n,1);
dep_all = cell(n,1);
pH_all  = cell(n,1);

for g = 1:n
    d = load(data_files{g});
    s = d.s;
    m = s.sdn >= date_start & s.sdn <= date_end;

    sdn_g    = s.sdn(m);
    dep_g    = s.depth(:, 1);
    pH_g     = s.(pH_field)(:, m);

    % 50 m binning
    bin_edges = 0 : bin_size : max(dep_g);
    bin_ctrs  = bin_edges(1:end-1) + bin_size/2;
    n_bins    = numel(bin_ctrs);
    pH_binned = NaN(n_bins, numel(sdn_g));
    for b = 1:n_bins
        rows = dep_g >= bin_edges(b) & dep_g < bin_edges(b+1);
        pH_binned(b, :) = mean(pH_g(rows, :), 1, 'omitnan');
    end

    sdn_all{g} = sdn_g;
    dep_all{g} = bin_ctrs(:);
    pH_all{g}  = pH_binned;
    fprintf('%s full range: %s  to  %s  |  in-range profiles: %d\n', ...
        panel_titles{g}, datestr(min(s.sdn),'yyyy-mm-dd'), datestr(max(s.sdn),'yyyy-mm-dd'), numel(sdn_g));
end

% Offset correction for Trinidad Head glider
pH_all{1} = pH_all{1} - 0.16;

% Merge second San Diego glider into panel 3
d2  = load("C:\Users\bwerb\Documents\ProposalFigures\data\25820301sat.mat");
s2  = d2.s;
fprintf('SD glider 2 full range: %s  to  %s  (%d profiles)\n', ...
    datestr(min(s2.sdn),'yyyy-mm-dd'), datestr(max(s2.sdn),'yyyy-mm-dd'), numel(s2.sdn));
fprintf('SD glider 2 fields: %s\n', strjoin(fieldnames(s2)', ', '));
m2  = s2.sdn >= date_start & s2.sdn <= date_end;
sdn_g2  = s2.sdn(m2);
fprintf('SD glider 2 profiles in range: %d\n', numel(sdn_g2));
dep_g2  = s2.depth(:, 1);
pH_g2   = s2.(pH_field)(:, m2);

bin_ctrs3  = dep_all{3};
n_bins3    = numel(bin_ctrs3);
pH_binned2 = NaN(n_bins3, numel(sdn_g2));
for b = 1:n_bins3
    rows = dep_g2 >= (bin_ctrs3(b) - bin_size/2) & dep_g2 < (bin_ctrs3(b) + bin_size/2);
    pH_binned2(b, :) = mean(pH_g2(rows, :), 1, 'omitnan');
end

[sdn_all{3}, sort_idx] = sort([sdn_all{3}(:); sdn_g2(:)]);
pH_all{3} = [pH_all{3}, pH_binned2];
pH_all{3} = pH_all{3}(:, sort_idx);

% Insert NaN columns at data gaps so contourf shows white instead of interpolating
for p = 1:n
    sdn_p = sdn_all{p}(:)';   % ensure row vector
    pH_p  = pH_all{p};
    gaps  = find(diff(sdn_p) > gap_thresh);
    for k = numel(gaps):-1:1   % back-to-front keeps indices valid
        t_mid = (sdn_p(gaps(k)) + sdn_p(gaps(k)+1)) / 2;
        sdn_p = [sdn_p(1:gaps(k)),   t_mid,                  sdn_p(gaps(k)+1:end)];
        pH_p  = [pH_p(:,1:gaps(k)),  NaN(size(pH_p,1), 1),   pH_p(:,gaps(k)+1:end)];
    end
    sdn_all{p} = sdn_p;
    pH_all{p}  = pH_p;
end

out_file = sprintf('pH_sections_3panel_%s_%s.png', ...
    char(datetime(date_start, 'ConvertFrom','datenum', 'Format','yyyyMMdd')), ...
    char(datetime(date_end,   'ConvertFrom','datenum', 'Format','yyyyMMdd')));

%% ── Figure ───────────────────────────────────────────────────────────────
fig = figure(1);  clf(fig);
set(fig, 'Units','inches', 'Position',[0.5 0.5 13.33 7.5]);
set(fig, 'DefaultAxesFontSize', 11);

contour_ticks      = linspace(date_start, date_end, 5);
contour_ticklabels = cellstr(datestr(contour_ticks, 'mmm dd'));

% Panel layout
fig_left  = 0.10;
fig_right = 0.88;
fig_top   = 0.95;
fig_bot   = 0.15;   % space for colorbar below bottom panel
gap       = 0.02;
panel_w   = fig_right - fig_left;
panel_h   = (fig_top - fig_bot - gap*(n-1)) / n;

for p = 1:n
    y_b = fig_top - p*panel_h - (p-1)*gap;
    ax  = axes(fig, 'Position', [fig_left  y_b  panel_w  panel_h]);

    contourf(ax, sdn_all{p}, dep_all{p}, pH_all{p}, 50, 'LineStyle', 'none')
    shading(ax, 'interp')
    axis(ax, 'ij')
    clim(ax, pH_clim)
    xlim(ax, [date_start  date_end])
    ylim(ax, [0  1000])
    colormap(ax, turbo)
    grid(ax, 'on')
    title(ax, panel_titles{p})
    colorbar(ax, 'off')

    if p == 2
        ylabel(ax, 'Depth (m)')
    end

    if p < n
        xticklabels(ax, {})
    else
        xticks(ax, contour_ticks)
        xticklabels(ax, contour_ticklabels)
        % Shared colorbar below bottom panel — restore Position after adding
        pos = ax.Position;
        cb  = colorbar(ax, 'southoutside');
        cb.Label.String = 'pH_{in situ} [total]';
        cb.FontSize = 11;
        ax.Position = pos;
    end
end

annotation(fig, 'textbox', [0 0.96 1 0.04], ...
    'HorizontalAlignment','center', 'VerticalAlignment','middle', ...
    'FontSize', 13, 'FontWeight','bold', 'EdgeColor','none')

%% ── Export ───────────────────────────────────────────────────────────────
set(fig, 'PaperUnits','inches', 'PaperSize',[13.33 7.5], ...
         'PaperPositionMode','manual', 'PaperPosition',[0 0 13.33 7.5]);
print(fig, '-dpng', '-r300', out_file)
fprintf('Saved %s\n', out_file)
