%% CUGN_pH.m  –  CUGN glider pH section comparison
%  Single landscape page: map (left) + stacked pHin sections (right)

%% ── Mission selection ────────────────────────────────────────────────────
% Set each entry to true/false to include or exclude that glider
%                  [25B68601,       26220301,         26120501  ]
%                  [Trinidad Head,  Monterey Bay,     San Diego ]
use_glider       = [true,           true,             true      ];

%% ── Date range control ───────────────────────────────────────────────────
% Leave as [] to use the full combined span of included deployments.
% Override with datenum values, e.g.:  datenum(2026, 3, 1)
date_start = datenum(2026, 1, 1);    % earliest date to show
date_end   = datenum(2026, 4, 14);    % latest   date to show

%% ── Glider / line metadata ───────────────────────────────────────────────
glider_labels  = {'Spray?: SN???',    'Spray2: SN203',    'Spray2: SN205'  };  % ← verify SN for 25B68601
legend_labels  = {'Trinidad Head, OR', 'Monterey Bay, CA', 'San Diego, CA' };
panel_titles   = {'Trinidad Head, OR', 'Monterey Bay, CA', 'San Diego, CA' };
markers        = {'s',                 'o',                '^'             };
track_colors   = {[1.0 0.5 0.1],       [1.0 1.0 0.2],      [0.3 0.9 1.0]  };  % orange, yellow, cyan

%% ── Load data ────────────────────────────────────────────────────────────
d1 = load("\\thalassa\ProjectLibrary\901805_Coastal_Biogeochemical_Sensing\" + ...
          "Spray_Data\25B68601\real_time_binned\25B68601sat.mat");

d2 = load("\\thalassa\ProjectLibrary\901805_Coastal_Biogeochemical_Sensing\" + ...
          "Spray_Data\26220301\real_time_binned\26220301sat.mat");

d3 = load("\\thalassa\ProjectLibrary\901805_Coastal_Biogeochemical_Sensing\" + ...
          "Spray_Data\26120501\real_time_binned\26120501sat.mat");

structs = {d1.s, d2.s, d3.s};

%% ── Determine unified date range ─────────────────────────────────────────
active_idx = find(use_glider);
sdn_cells  = cellfun(@(s) s.sdn(:), structs(active_idx), 'UniformOutput', false);
all_sdn    = cat(1, sdn_cells{:});
fin        = isfinite(all_sdn);
if isempty(date_start),  date_start = min(all_sdn(fin));  end
if isempty(date_end),    date_end   = max(all_sdn(fin));  end

%% ── Filter profiles and extract per-glider data ──────────────────────────
sdn  = cell(1,3);  lat  = cell(1,3);  lon  = cell(1,3);
dep  = cell(1,3);

for k = 1:3
    s    = structs{k};
    m    = s.sdn >= date_start & s.sdn <= date_end;
    sdn{k} = s.sdn(m);
    lat{k} = s.lat(m);
    lon{k} = s.lon(m);
    dep{k} = s.depth(:,1);
end

%% ── Shared date tick labels ──────────────────────────────────────────────
date_ticks         = linspace(date_start, date_end, 10);
date_ticklabels    = cellstr(datestr(date_ticks, 'yyyy-mm-dd'));
contour_ticks      = linspace(date_start, date_end, 4);
contour_ticklabels = cellstr(datestr(contour_ticks, 'mmm dd'));

%% ── Figure title ─────────────────────────────────────────────────────────
if all(use_glider)
    fig_title = 'California Current System Glider pH Transects 2026';
else
    fig_title = 'California Underwater Glider Network pH Transects';
end

%% ── Panel layout (computed from number of active gliders) ────────────────
n_active  = numel(active_idx);
fig_top   = 0.89;   % normalized y just below supertitle
fig_bot   = 0.17;   % normalized y above horizontal colorbar
gap       = 0.035;  % gap between panels
panel_h   = (fig_top - fig_bot - gap*(n_active-1)) / n_active;
panel_x   = 0.57;
panel_w   = 0.37;

%% ── Extract pH and compute shared color limits ───────────────────────────
pH_field = 'pHin';
cb_label = 'pH_{in situ} [total]';
out_file = 'CUGN_pHin.png';

pH = cell(1,3);
for k = 1:3
    s    = structs{k};
    m    = s.sdn >= date_start & s.sdn <= date_end;
    pH{k} = s.(pH_field)(:, m);
end
pH_cells = cellfun(@(x) x(:), pH(active_idx), 'UniformOutput', false);
all_pH   = cat(1, pH_cells{:});
pH_clim = [7.5  8.1];

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

% ── MAP  (left half, full height) ─────────────────────────────────────
ax_map = geoaxes(fig, 'Position', [0.03 0.08 0.42 0.84]);
hold(ax_map, 'on')
for k = active_idx
    geoscatter(ax_map, lat{k}, lon{k}, 20, sdn{k}, 'filled', ...
               'Marker', markers{k}, 'DisplayName', glider_labels{k});
end
colormap(ax_map, jet)
clim(ax_map, [date_start  date_end])
cb_map            = colorbar(ax_map);
cb_map.Ticks      = date_ticks;
cb_map.TickLabels = date_ticklabels;
cb_map.FontSize   = 7;
ax_map.LatitudeLabel.String  = '';
ax_map.LongitudeLabel.String = '';
ax_map.FontSize = 8;
legend(ax_map, legend_labels(active_idx), 'Location','northeast')
try; geobasemap('satellite'); catch, end

% ── Contourf panels (one per active glider, top to bottom) ────────────
for p = 1:n_active
    k   = active_idx(p);
    y_b = fig_top - p*panel_h - (p-1)*gap;   % bottom of this panel

    ax = axes(fig, 'Position', [panel_x  y_b  panel_w  panel_h]);
    contourf(ax, sdn{k}, dep{k}, pH{k}, 5)
    shading(ax, 'interp')
    axis(ax, 'ij')
    clim(ax, pH_clim)
    xlim(ax, [date_start  date_end])
    ylim(ax, [0  1000])
    xticks(ax, contour_ticks)
    grid(ax, 'on')
    ax.GridColor      = 'k';
    ax.GridAlpha      = 0.3;
    ax.Layer          = 'top';
    title(ax, panel_titles{k})

    if p < n_active
        % Upper panels: hide x tick labels, no colorbar
        xticklabels(ax, {})
        colorbar(ax, 'off')
    else
        % Bottom panel: show date labels and horizontal colorbar
        xticklabels(ax, contour_ticklabels)
        pos = ax.Position;
        cb  = colorbar(ax, 'southoutside');
        cb.Label.String = cb_label;
        cb.FontSize = 7;
        ax.Position = pos;
    end
end

% ── Super-title ────────────────────────────────────────────────────────
annotation(fig, 'textbox', [0 0.93 1 0.06], ...
    'String', fig_title, ...
    'HorizontalAlignment','center', 'VerticalAlignment','middle', ...
    'FontSize', 14, 'FontWeight','bold', 'EdgeColor','none')

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