%% gulf_stream_section.m – Contour and scatter plots for Gulf Stream glider
%  (glider 25B20901sat, profiles 41-91)

%% ── Settings ─────────────────────────────────────────────────────────────────
data_file  = "C:\Users\bwerb\Documents\ProposalFigures\data\26321301sat.mat";
out_dir    = "C:\Users\bwerb\Documents\ProposalFigures\gulf_stream_section";
if ~isfolder(out_dir), mkdir(out_dir); end

prof_idx   = 41:92;
pH_field   = 'pHin';   % switch to 'pHin' to use the QC'd field
bin_size   = 10;     % m
gap_thresh = 10;      % days
depth_max  = 1000;   % m

%% ── Load ─────────────────────────────────────────────────────────────────────
d = load(data_file);
s = d.s;

sdn = s.sdn(prof_idx);
dep = s.depth(:, 1);

fields = {pH_field,'doxy','tc','psal','pHin_esper','pHin_canb','pco2in','dic_canb'};
raw = struct();
for k = 1:numel(fields)
    f   = fields{k};
    arr = s.(f);
    if ismatrix(arr) && size(arr,2) > 1
        raw.(f) = arr(:, prof_idx);
    else
        raw.(f) = arr(prof_idx);
    end
end

n_profs = numel(sdn);

%% ── Depth binning (50 m) ─────────────────────────────────────────────────────
bin_edges = 0 : bin_size : depth_max;
bin_ctrs  = bin_edges(1:end-1) + bin_size/2;
n_bins    = numel(bin_ctrs);

binned = struct();
for k = 1:numel(fields)
    f   = fields{k};
    arr = raw.(f);
    if ~ismatrix(arr) || size(arr,2) == 1, continue; end
    out = NaN(n_bins, n_profs);
    for b = 1:n_bins
        rows = dep >= bin_edges(b) & dep < bin_edges(b+1);
        out(b,:) = mean(arr(rows,:), 1, 'omitnan');
    end
    binned.(f) = out;
end

%% ── Contour plots ────────────────────────────────────────────────────────────
save_contour(sdn, binned.(pH_field),   bin_ctrs, turbo,   [],          gap_thresh, ...
    'Gulf Stream – pH_{in situ}', 'pH_{in situ} [total]', ...
    fullfile(out_dir, 'pHin_section.png'));

save_contour(sdn, binned.doxy,     bin_ctrs, parula,  [100 220],   gap_thresh, ...
    'Gulf Stream – Dissolved Oxygen', 'DO (\mumol kg^{-1})', ...
    fullfile(out_dir, 'doxy_section.png'));

save_contour(sdn, binned.tc,       bin_ctrs, jet,     [5 30],      gap_thresh, ...
    'Gulf Stream – Temperature', 'Temperature (^{\circ}C)', ...
    fullfile(out_dir, 'tc_section.png'));

save_contour(sdn, binned.psal,     bin_ctrs, parula,  [33.5 37],   gap_thresh, ...
    'Gulf Stream – Salinity', 'Salinity (PSS)', ...
    fullfile(out_dir, 'psal_section.png'));

save_contour(sdn, binned.pco2in,   bin_ctrs, hot,     [200 700],   gap_thresh, ...
    'Gulf Stream – pCO_2', 'pCO_2 (\muatm)', ...
    fullfile(out_dir, 'pco2_section.png'));

save_contour(sdn, binned.dic_canb, bin_ctrs, parula,  [1800 2200], gap_thresh, ...
    'Gulf Stream – DIC (CANYON-B)', 'DIC (\mumol kg^{-1})', ...
    fullfile(out_dir, 'dic_section.png'));

%% ── ΔpH diverging contour ────────────────────────────────────────────────────
delta_esper_bin = binned.(pH_field) - binned.pHin_esper;
save_diverging(sdn, delta_esper_bin, bin_ctrs, gap_thresh, ...
    'Gulf Stream – \DeltapH (measured − ESPER)', '\DeltapH', ...
    fullfile(out_dir, 'deltapH_esper_section.png'));

delta_canb_bin = binned.(pH_field) - binned.pHin_canb;
save_diverging(sdn, delta_canb_bin, bin_ctrs, gap_thresh, ...
    'Gulf Stream – \DeltapH (measured − CANYON-B)', '\DeltapH', ...
    fullfile(out_dir, 'deltapH_canb_section.png'));

%% ── Scatter plots ────────────────────────────────────────────────────────────
delta_esper_raw = raw.(pH_field) - raw.pHin_esper;
delta_canb_raw  = raw.(pH_field) - raw.pHin_canb;

save_scatter_delta(delta_esper_raw, delta_canb_raw, dep, depth_max, ...
    '\DeltapH_{in} (pH_{in} - Algorithm)', ...
    'Gulf Stream – \DeltapH vs Depth', ...
    fullfile(out_dir, 'scatter_deltapH.png'));

fprintf('\nAll figures saved to: %s\n', out_dir);

%% ════════════════════════════════════════════════════════════════════════════
%  Local functions
%  ════════════════════════════════════════════════════════════════════════════

function [t_out, Z_out] = add_gaps(t, Z, thresh)
    t_out = t(:)';
    Z_out = Z;
    gaps  = find(diff(t_out) > thresh);
    for k = numel(gaps):-1:1
        i     = gaps(k);
        t_mid = (t_out(i) + t_out(i+1)) / 2;
        t_out = [t_out(1:i),   t_mid,                t_out(i+1:end)];
        Z_out = [Z_out(:,1:i), NaN(size(Z_out,1),1), Z_out(:,i+1:end)];
    end
end

function cmap = redblue(n)
    if nargin < 1, n = 256; end
    half  = floor(n/2);
    blues = [linspace(0,1,half)', linspace(0,1,half)', ones(half,1)];
    reds  = [ones(n-half,1), linspace(1,0,n-half)', linspace(1,0,n-half)'];
    cmap  = [blues; reds];
end

function save_contour(sdn, Z, bin_ctrs, cmap, clim_val, gap_thresh, ttl, cbar_lbl, out_path)
    [sdn_p, Z_p] = add_gaps(sdn, Z, gap_thresh);

    date_start = min(sdn_p(isfinite(sdn_p)));
    date_end   = max(sdn_p(isfinite(sdn_p)));

    fig = figure('Visible', 'off');  clf(fig);
    set(fig, 'Units','inches', 'Position',[0.5 0.5 13.33 7.5]);
    set(fig, 'DefaultAxesFontSize', 14);

    ax = axes(fig);
    contourf(ax, sdn_p, bin_ctrs(:), Z_p, 50, 'LineStyle','none')
    shading(ax, 'interp')
    axis(ax, 'ij')
    if ~isempty(clim_val), clim(ax, clim_val); end
    xlim(ax, [date_start  date_end])
    ylim(ax, [0  max(bin_ctrs)])
    colormap(ax, cmap)
    grid(ax, 'on')
    title(ax, ttl)
    ylabel(ax, 'Depth (m)')

    n_ticks = 6;
    ticks   = linspace(date_start, date_end, n_ticks);
    xticks(ax, ticks)
    xticklabels(ax, cellstr(datetime(ticks, 'ConvertFrom','datenum', 'Format','MMM dd yyyy')))

    pos = ax.Position;
    cb  = colorbar(ax, 'eastoutside');
    cb.Label.String = cbar_lbl;
    cb.FontSize = 10;
    ax.Position = pos;

    set(fig, 'PaperUnits','inches', 'PaperSize',[13.33 7.5], ...
             'PaperPositionMode','manual', 'PaperPosition',[0 0 13.33 7.5]);
    print(fig, '-dpng', '-r300', out_path)
    close(fig);
    [~, fname, ext] = fileparts(out_path);
    fprintf('Saved  %s\n', [fname ext]);
end

function save_diverging(sdn, Z, bin_ctrs, gap_thresh, ttl, cbar_lbl, out_path)
    finite_vals = Z(isfinite(Z));
    vmax = max(prctile(abs(finite_vals), 98), 1e-3);
    save_contour(sdn, Z, bin_ctrs, redblue(256), [-vmax vmax], gap_thresh, ...
                 ttl, cbar_lbl, out_path);
end

function save_scatter_delta(Z1, Z2, dep, depth_max, xlabel_str, ttl, out_path)
    n_profs = size(Z1, 2);
    d       = dep(:);

    fig = figure('Visible', 'off');  clf(fig);
    set(fig, 'Units','inches', 'Position',[0.5 0.5 6 8]);
    set(fig, 'DefaultAxesFontSize', 14);

    ax = axes(fig);
    hold(ax, 'on')

    c1 = [0.27 0.51 0.71];
    c2 = [0.85 0.33 0.10];

    for p = 1:n_profs
        z1 = Z1(:, p);  mask1 = isfinite(d) & isfinite(z1) & (d <= depth_max);
        z2 = Z2(:, p);  mask2 = isfinite(d) & isfinite(z2) & (d <= depth_max);

        h1 = plot(ax, z1(mask1), d(mask1), '-o', 'Color', [c1 0.4], ...
                  'MarkerSize', 2, 'MarkerFaceColor', c1, 'LineWidth', 0.6);
        h2 = plot(ax, z2(mask2), d(mask2), '-o', 'Color', [c2 0.4], ...
                  'MarkerSize', 2, 'MarkerFaceColor', c2, 'LineWidth', 0.6);

        if p == 1
            h1.DisplayName = 'ESPER';
            h2.DisplayName = 'CANYON-B';
        else
            h1.HandleVisibility = 'off';
            h2.HandleVisibility = 'off';
        end
    end

    xline(ax, 0, 'k--', 'LineWidth', 0.8, 'HandleVisibility','off')
    legend(ax, 'Location','best', 'FontSize', 11)
    axis(ax, 'ij')
    ylim(ax, [0  depth_max])
    xlabel(ax, xlabel_str)
    ylabel(ax, 'Depth (m)')
    title(ax, ttl)
    grid(ax, 'on')

    set(fig, 'PaperUnits','inches', 'PaperSize',[6 8], ...
             'PaperPositionMode','manual', 'PaperPosition',[0 0 6 8]);
    print(fig, '-dpng', '-r300', out_path)
    close(fig);
    [~, fname, ext] = fileparts(out_path);
    fprintf('Saved  %s\n', [fname ext]);
end
