%% gulf_stream_timeseries.m  –  Full time-series pH section for Gulf Stream

%% ── Inputs ───────────────────────────────────────────────────────────────
data_files = { ...
    "C:\Users\bwerb\Documents\ProposalFigures\data\25B20901sat.mat"; ...
    "C:\Users\bwerb\Documents\ProposalFigures\data\26321301sat.mat" };
pH_field   = 'pHin';
pH_clim    = [7.85  8.05];
bin_size   = 50;
gap_thresh = 2;    % days — gaps longer than this render as white

%% ── Load, filter, and bin first glider ──────────────────────────────────
d1  = load(data_files{1});
s1  = d1.s;
sdn_g1 = s1.sdn;
dep_g1 = s1.depth(:, 1);
pH_g1  = s1.(pH_field);

date_start = min(sdn_g1);
date_end   = max(sdn_g1);

bin_edges = 0 : bin_size : max(dep_g1);
bin_ctrs  = bin_edges(1:end-1) + bin_size/2;
n_bins    = numel(bin_ctrs);
pH_binned1 = NaN(n_bins, numel(sdn_g1));
for b = 1:n_bins
    rows = dep_g1 >= bin_edges(b) & dep_g1 < bin_edges(b+1);
    pH_binned1(b, :) = mean(pH_g1(rows, :), 1, 'omitnan');
end

fprintf('GS glider 1 range: %s  to  %s  (%d profiles)\n', ...
    char(datetime(date_start,'ConvertFrom','datenum','Format','yyyy-MM-dd')), ...
    char(datetime(date_end,  'ConvertFrom','datenum','Format','yyyy-MM-dd')), numel(sdn_g1));

%% ── Load and bin second glider ───────────────────────────────────────────
d2  = load(data_files{2});
s2  = d2.s;
sdn_g2 = s2.sdn;
dep_g2 = s2.depth(:, 1);
pH_g2  = s2.(pH_field);

fprintf('GS glider 2 range: %s  to  %s  (%d profiles)\n', ...
    char(datetime(min(sdn_g2),'ConvertFrom','datenum','Format','yyyy-MM-dd')), ...
    char(datetime(max(sdn_g2),'ConvertFrom','datenum','Format','yyyy-MM-dd')), numel(sdn_g2));

pH_binned2 = NaN(n_bins, numel(sdn_g2));
for b = 1:n_bins
    rows = dep_g2 >= (bin_ctrs(b) - bin_size/2) & dep_g2 < (bin_ctrs(b) + bin_size/2);
    pH_binned2(b, :) = mean(pH_g2(rows, :), 1, 'omitnan');
end

%% ── Merge and update date range ──────────────────────────────────────────
[sdn_merged, sort_idx] = sort([sdn_g1(:); sdn_g2(:)]);
pH_merged = [pH_binned1, pH_binned2];
pH_merged = pH_merged(:, sort_idx);

date_start = min(sdn_merged);
date_end   = max(sdn_merged);

%% ── Insert NaN columns at data gaps ─────────────────────────────────────
sdn_p = sdn_merged(:)';
pH_p  = pH_merged;
gaps  = find(diff(sdn_p) > gap_thresh);
for k = numel(gaps):-1:1
    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(n_bins, 1),          pH_p(:,gaps(k)+1:end)];
end

%% ── Figure ───────────────────────────────────────────────────────────────
fig = figure(1);  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(:), pH_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, 'Gulf Stream')
ylabel(ax, 'Depth (m)')

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

pos = ax.Position;
cb  = colorbar(ax, 'eastoutside');
cb.Label.String = 'pH_{in situ} [total]';
cb.FontSize = 10;
ax.Position = pos;

%% ── Export ───────────────────────────────────────────────────────────────
out_file = sprintf('gulf_stream_timeseries_%s_%s.png', ...
    char(datetime(date_start, 'ConvertFrom','datenum', 'Format','yyyyMMdd')), ...
    char(datetime(date_end,   'ConvertFrom','datenum', 'Format','yyyyMMdd')));
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)
