%% Here's a changelog of everything that's happened with the pH data:
% Dive 90: pH sample climb-only
% Dive 123: pumptime changed from 5000 to 2500 ms
% Dive 124-131: pH sensor off
% Dive 132: pH sensor on
% Dive 158: Oxygen, fluorescence, pH sampling only on climb
%% Create struct for ascent and descent profiles
id_25B68601 = Spray2Data('25B68601');
id_25B68601.ascend = 1;
id_25B68601.descend = 1;
id_25B68601.both = 0;
id_25B68601.copyNc;
id_25B68601.extractNc;
id_25B68601.buildStruct;
id_25B68601.addCarbonateParameters;
id_25B68601.buildTable;
s = id_25B68601.s;
s.phase = id_25B68601.bindata.phase';
s.divenumber = id_25B68601.bindata.divenumber';
data = id_25B68601.data;
save('Ascent_Descent_Compare_data.mat',"data",'-mat')
save('Ascent_Descent_Compare_s.mat',"s",'-mat')
%% Alternatively load previous s, data structs
load("Ascent_Descent_Compare_s.mat")
load("Ascent_Descent_Compare_data.mat")
%% Verify t/s are interpolating correctly onto pH time grid
phase = data.ctd.phase{end,:};
iascend = phase == 0;
fig = figure;
subplot(121)
plot(data.ph.s{:,end},data.ph.depth{end,:},'ko')
hold on
plot(data.ctd.s{end,1}(iascend),data.ctd.depth{end,1}(iascend),'r.')
set(gca,'ydir','reverse')
title('psal: pH vs ctd')
subplot(122)
plot(data.ph.t{:,end},data.ph.depth{end,:},'ko')
hold on
plot(data.ctd.t{end,1}(iascend),data.ctd.depth{end,1}(iascend),'r.')
set(gca,'ydir','reverse')
title('t: pH vs ctd')
%% Plot ascent vs descent
iascend = s.phase == 0;
idescend = s.phase == 1;
idive90 = s.divenumber >= 90; % Changed to ascend only sampling for pH
idive123 = s.divenumber >= 123; % Changed pump time 5s to 2.5s
idive124_131 = s.divenumber >= 124 & s.divenumber <= 131; % pH sensor off
idive158 = s.divenumber >= 158; % Changed to ascend only sampling for O2, chl, pH
%% Compare Ascent dives before and after ascent only
fig = figure;
tl = tiledlayout(1,1);
nexttile
p1 = plot(s.pHin(:,iascend & ~idive158),s.depth(:,iascend & ~idive158),'k>');
hold on;
p2 = plot(s.pHin(:,iascend & idive158),s.depth(:,iascend & idive158),'r<');
legend([p1(1), p2(1)],{'dnum<158', 'dnum>=158'},Location="southeast");
title('Ascend')
set(gca,'ydir','reverse')
xlabel('pHin')
% No clear difference
%% Compare 2.5 to 5 second pump times independent of ascend/descend
fig = figure;
tl = tiledlayout(3,2);
title(tl,'Pump on Time: 5 vs 2.5 s')
nexttile
p1 = plot(s.pHin(:,idescend & ~idive123),s.depth(:,idescend & ~idive123),'k-');
hold on;
p2 = plot(s.pHin(:,idescend & idive123),s.depth(:,idescend & idive123),'r-');
legend([p1(1), p2(1)],{'5 s', '2.5 s'},Location="southeast");
title('Descend')
xlabel('pHin')
set(gca,'ydir','reverse')
nexttile
p1 = plot(s.pHin(:,iascend & ~idive123),s.depth(:,iascend & ~idive123),'k-');
hold on;
p2 = plot(s.pHin(:,iascend & idive123),s.depth(:,iascend & idive123),'r-');
legend([p1(1), p2(1)],{'5 s', '2.5 s'},Location="southeast");
title('Ascend')
xlabel('pHin')
set(gca,'ydir','reverse')
nexttile
p1 = plot(s.psal(:,idescend & ~idive123),s.depth(:,idescend & ~idive123),'k-');
hold on;
p2 = plot(s.psal(:,idescend & idive123),s.depth(:,idescend & idive123),'r-');
legend([p1(1), p2(1)],{'5 s', '2.5 s'},Location="southeast");
title('Descend')
xlabel('psal')
set(gca,'ydir','reverse')
nexttile
p1 = plot(s.psal(:,iascend & ~idive123),s.depth(:,iascend & ~idive123),'k-');
hold on;
p2 = plot(s.psal(:,iascend & idive123),s.depth(:,iascend & idive123),'r-');
legend([p1(1), p2(1)],{'5 s', '2.5 s'},Location="southeast");
title('Ascend')
xlabel('psal')
set(gca,'ydir','reverse')
% scatter
nexttile
p1 = scatter(s.psal(:,idescend & ~idive123),s.pHin(:,idescend & ~idive123),'k.');
hold on;
p2 = scatter(s.psal(:,idescend & idive123),s.pHin(:,idescend & idive123),'r.');
legend([p1(1), p2(1)],{'5 s', '2.5 s'},Location="southwest");
title('Descend')
xlabel('pHin vs psal')
% set(gca,'ydir','reverse')
nexttile
p1 = scatter(s.psal(:,iascend & ~idive123),s.pHin(:,iascend & ~idive123),'k.');
hold on;
p2 = scatter(s.psal(:,iascend & idive123),s.pHin(:,iascend & idive123),'r.');
legend([p1(1), p2(1)],{'5 s', '2.5 s'},Location="southwest");
title('Ascend')
xlabel('pHin vs psal')
% set(gca,'ydir','reverse')
% Something is going on with descend dives but not sure if it's related to
% pump on time
%% Look at vrse not pHin. psal might be the issue
fig = figure;
tl = tiledlayout(1,2);
title(tl,'Pump on Time: 5 vs 2.5 s')
nexttile
p1 = plot(s.vrse(:,idescend & ~idive123),s.depth(:,idescend & ~idive123),'k-');
hold on;
p2 = plot(s.vrse(:,idescend & idive123),s.depth(:,idescend & idive123),'r-');
legend([p1(1), p2(1)],{'5 s', '2.5 s'},Location="southeast");
title('Descend')
xlabel('vrse')
set(gca,'ydir','reverse')
nexttile
p1 = plot(s.vrse(:,iascend & ~idive123),s.depth(:,iascend & ~idive123),'k-');
hold on;
p2 = plot(s.vrse(:,iascend & idive123),s.depth(:,iascend & idive123),'r-');
legend([p1(1), p2(1)],{'5 s', '2.5 s'},Location="southeast");
title('Ascend')
xlabel('vrse')
set(gca,'ydir','reverse')
%% Compare the distributions of descend data. Upper 200 seems bimodal
fig = figure;
clf
tl = tiledlayout(1,2);
nexttile

% Get unique dive numbers for descents
unique_dives = unique(s.divenumber(idescend));
n_dives = length(unique_dives);

% Create a colormap
colors = jet(n_dives);  % or use parula, turbo, etc.

hold on;
% Plot each dive separately
for i = 1:158
    dive_num = unique_dives(i);
    idive = (s.divenumber == dive_num) & idescend;
    
    plot(s.pHin(:,idive), s.depth(:,idive), '-', ...
         'Color', colors(i,:), 'MarkerSize', 10);
end

title('Descend')
xlabel('pHin')
ylabel('Depth')
set(gca,'ydir','reverse')

% Add colorbar
colormap(colors);
c = colorbar;
c.Label.String = 'Dive Number';
caxis([min(unique_dives) max(unique_dives)]);

nexttile

% Get unique dive numbers for descents
unique_dives = unique(s.divenumber(iascend));
n_dives = length(unique_dives);

% Create a colormap
colors = jet(n_dives);  % or use parula, turbo, etc.

hold on;
% Plot each dive separately
for i = 1:158
    dive_num = unique_dives(i);
    idive = (s.divenumber == dive_num) & iascend;
    
    plot(s.pHin(:,idive), s.depth(:,idive), '-', ...
         'Color', colors(i,:), 'MarkerSize', 10);
end

title('Ascend')
xlabel('pHin')
ylabel('Depth')
set(gca,'ydir','reverse')

% Add colorbar
colormap(colors);
c = colorbar;
c.Label.String = 'Dive Number';
caxis([min(unique_dives) max(unique_dives)]);
%% Something happens to the ascend profiles after dive 58
idive48_58 = s.divenumber >= 48 & s.divenumber <= 58;
idive59_69 = s.divenumber >= 59 & s.divenumber <= 69;
fig = figure;
tl = tiledlayout(3,2);
nexttile
p1 = plot(s.pHin(:,iascend & idive48_58),s.depth(:,iascend & idive48_58),'k>');
hold on;
p2 = plot(s.pHin(:,idescend & idive48_58),s.depth(:,idescend & idive48_58),'r*');
legend([p1(1), p2(1)],{'ascend', 'descend'},Location="southeast");
title('divenumber: 48-58')
set(gca,'ydir','reverse')
xlabel('pHin')
nexttile
p1 = plot(s.pHin(:,iascend & idive59_69),s.depth(:,iascend & idive59_69),'k>');
hold on;
p2 = plot(s.pHin(:,idescend & idive59_69),s.depth(:,idescend & idive59_69),'r*');
legend([p1(1), p2(1)],{'ascend', 'descend'},Location="southeast");
title('divenumber: 59-69')
set(gca,'ydir','reverse')
xlabel('pHin')
nexttile
p1 = plot(s.tc(:,iascend & idive48_58),s.depth(:,iascend & idive48_58),'k>');
hold on;
p2 = plot(s.tc(:,idescend & idive48_58),s.depth(:,idescend & idive48_58),'r*');
legend([p1(1), p2(1)],{'ascend', 'descend'},Location="southeast");
title('divenumber: 48-58')
set(gca,'ydir','reverse')
xlabel('tc')
% No clear difference
nexttile
p1 = plot(s.tc(:,iascend & idive59_69),s.depth(:,iascend & idive59_69),'k>');
hold on;
p2 = plot(s.tc(:,idescend & idive59_69),s.depth(:,idescend & idive59_69),'r*');
legend([p1(1), p2(1)],{'ascend', 'descend'},Location="southeast");
title('divenumber: 59-69')
set(gca,'ydir','reverse')
xlabel('tc')
nexttile
p1 = plot(s.psal(:,iascend & idive48_58),s.depth(:,iascend & idive48_58),'k>');
hold on;
p2 = plot(s.psal(:,idescend & idive48_58),s.depth(:,idescend & idive48_58),'r*');
legend([p1(1), p2(1)],{'ascend', 'descend'},Location="southeast");
title('divenumber: 48-58')
set(gca,'ydir','reverse')
xlabel('psal')
% No clear difference
nexttile
p1 = plot(s.psal(:,iascend & idive59_69),s.depth(:,iascend & idive59_69),'k>');
hold on;
p2 = plot(s.psal(:,idescend & idive59_69),s.depth(:,idescend & idive59_69),'r*');
legend([p1(1), p2(1)],{'ascend', 'descend'},Location="southeast");
title('divenumber: 59-69')
set(gca,'ydir','reverse')
xlabel('psal')
%% Plot all dive ranges in chunks of 10 (0-90)
fig = figure('Position', [100, 100, 2000, 800]); % Wide figure for horizontal layout
tl = tiledlayout(2, 9, 'TileSpacing', 'compact', 'Padding', 'compact');

% Define dive ranges in chunks of 10
dive_ranges = [0 10; 10 20; 20 30; 30 40; 40 50; 50 60; 60 70; 70 80; 80 90];

for i = 1:size(dive_ranges, 1)
    dive_start = dive_ranges(i, 1);
    dive_end = dive_ranges(i, 2);
    
    % Create logical index for current dive range
    idive_range = s.divenumber >= dive_start & s.divenumber <= dive_end;
    
    % Plot pHin (top row)
    nexttile
    p1 = plot(s.pHin(:, iascend & idive_range), s.depth(:, iascend & idive_range), 'k>');
    hold on;
    p2 = plot(s.pHin(:, idescend & idive_range), s.depth(:, idescend & idive_range), 'r*');
    if i == 1
        legend([p1(1), p2(1)], {'ascend', 'descend'}, 'Location', 'southeast');
    end
    title(sprintf('%d-%d', dive_start, dive_end))
    set(gca, 'ydir', 'reverse')
    xlabel('pHin')
    if i == 1
        ylabel('Depth')
    end
end

% Second row: tc and psal on same plot
for i = 1:size(dive_ranges, 1)
    dive_start = dive_ranges(i, 1);
    dive_end = dive_ranges(i, 2);
    
    % Create logical index for current dive range
    idive_range = s.divenumber >= dive_start & s.divenumber <= dive_end;
    
    % Plot tc and psal together
    nexttile
    yyaxis left
    p1 = plot(s.tc(:, iascend & idive_range), s.depth(:, iascend & idive_range), 'k>');
    hold on;
    p2 = plot(s.tc(:, idescend & idive_range), s.depth(:, idescend & idive_range), 'r*');
    set(gca, 'ydir', 'reverse')
    xlabel('tc')
    ylabel('Depth')
    
    yyaxis right
    p3 = plot(s.psal(:, iascend & idive_range), s.depth(:, iascend & idive_range), 'b>');
    p4 = plot(s.psal(:, idescend & idive_range), s.depth(:, idescend & idive_range), 'm*');
    set(gca, 'ydir', 'reverse')
    xlabel('tc / psal')
    
    if i == 1
        legend([p1(1), p2(1), p3(1), p4(1)], {'tc asc', 'tc desc', 'psal asc', 'psal desc'}, 'Location', 'southeast');
        ylabel('Depth')
    end
end

% Add overall title
title(tl, 'Dive Profiles by Chunks (0-90): pHin (top), tc/psal (bottom)')
%% Do this for the full mission
%% Plot in multiple figures and save
dive_ranges = [0 10; 10 20; 20 30; 30 40; 40 50; 50 60; 60 70; 70 80; 80 90; ...
               90 100; 100 110; 110 120; 120 130; 130 140; 140 150; 150 160; ...
               160 170; 170 180; 180 190];

chunks_per_fig = 5;
num_figs = ceil(size(dive_ranges, 1) / chunks_per_fig);

for fig_num = 1:num_figs
    fig = figure('Position', [100, 100, 1400, 800]);
    tl = tiledlayout(2, chunks_per_fig, 'TileSpacing', 'compact', 'Padding', 'compact');
    
    start_idx = (fig_num - 1) * chunks_per_fig + 1;
    end_idx = min(fig_num * chunks_per_fig, size(dive_ranges, 1));
    
    for i = start_idx:end_idx
        dive_start = dive_ranges(i, 1);
        dive_end = dive_ranges(i, 2);
        idive_range = s.divenumber >= dive_start & s.divenumber <= dive_end;
        
        % Plot pHin (top row)
        nexttile
        p1 = plot(s.pHin(:, iascend & idive_range), s.depth(:, iascend & idive_range), 'k>');
        hold on;
        p2 = plot(s.pHin(:, idescend & idive_range), s.depth(:, idescend & idive_range), 'r*');
        if i == start_idx
            legend([p1(1), p2(1)], {'ascend', 'descend'}, 'Location', 'southeast');
        end
        title(sprintf('%d-%d', dive_start, dive_end))
        set(gca, 'ydir', 'reverse')
        xlabel('pHin')
        if i == start_idx
            ylabel('Depth')
        end
    end
    
    % Second row: tc and psal
    for i = start_idx:end_idx
        dive_start = dive_ranges(i, 1);
        dive_end = dive_ranges(i, 2);
        idive_range = s.divenumber >= dive_start & s.divenumber <= dive_end;
        
        nexttile
        yyaxis left
        p1 = plot(s.tc(:, iascend & idive_range), s.depth(:, iascend & idive_range), 'k>');
        hold on;
        p2 = plot(s.tc(:, idescend & idive_range), s.depth(:, idescend & idive_range), 'r*');
        set(gca, 'ydir', 'reverse')
        xlabel('tc')
        ylabel('Depth')
        
        yyaxis right
        p3 = plot(s.psal(:, iascend & idive_range), s.depth(:, iascend & idive_range), 'b>');
        p4 = plot(s.psal(:, idescend & idive_range), s.depth(:, idescend & idive_range), 'm*');
        set(gca, 'ydir', 'reverse')
        xlabel('tc / psal')
        
        if i == start_idx
            legend([p1(1), p2(1), p3(1), p4(1)], {'tc asc', 'tc desc', 'psal asc', 'psal desc'}, 'Location', 'southeast');
            ylabel('Depth')
        end
    end
    
    title(tl, sprintf('Dive Profiles: Chunks %d-%d', dive_ranges(start_idx, 1), dive_ranges(end_idx, 2)))
    
    % Save figure
    filename = sprintf('dive_profiles_%d_%d.png', dive_ranges(start_idx, 1), dive_ranges(end_idx, 2));
    saveas(fig, filename);
    % Or for higher quality:
    % exportgraphics(fig, filename, 'Resolution', 300);
    
    fprintf('Saved: %s\n', filename);
end
%% Likely issue is a pump malfunction beginning on dive 59
% Email:
% Hi Jace and all, sorry for the delay, Yui and I have been reviewing the deployment data.
% 
% I have attached figures comparing the ascend/descend dives as chunks of 10 over the entire deployment. After dive 59, there is a significant change in the pH sensor measurements. The ascend measurements are consistently biased low, which could indicate that the sensor is reading unflushed deep water. Could the pump be weaker on the ascent or not functioning at all? Does this dive number coincide with the issues affecting the science board?
% 
% For reference, here is the change log that you sent earlier:
% •	Dive 90: pH sample climb-only
% •	Dive 123: pumptime changed from 5000 to 2500 ms
% •	Dive 124-131: pH sensor off
% •	Dive 132: pH sensor on
% •	Dive 158: Oxygen, fluorescence, pH sampling only on climb
% 
% Let me know what you think!
% 
% Thanks,
% Ben W
