%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Supplemntary Figure 6 - chematic profiles for decomposition illustration %%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 1x6 panels (adjust for layout in paper)
% Panels:
% 1) 2015 
% 2) 2100 
% 3) 2015 with Total Phyto scaled to equal 2100 Total Phyto (Zoo unchanged)
% 4) 2015 with Total Zoo   scaled to equal 2100 Total Zoo   (Phyto unchanged)
% 5) 2015 Total-P and Total-Z SHAPES conserved; at each depth, redistribute
%    Small/Large fractions to match 2100 fractions (per-depth) [Totals conserved]
% 6) 2015 totals & Small/Large RATIOS conserved; replace Total-P and Total-Z
%    vertical SHAPES with 2100 (normalized) shapes
%
% All panels: same styles, axis bounds, MLD lines, annotations.
cd /Users/twrohr/Desktop/Matlab/Matlab/Grazing/CMIP6_Comp/NCC_Figs
close all

%% ---------------- Parameters ----------------
zmax     = 150; dz = 1; z = (0:dz:zmax)'; ydir_rev = true;

% Mixed layer depths
MLD_2015 = 60;
MLD_2100 = 30;

% -------- Biomass shapes (keep your current choices) --------
% 2015: monotonic decreasing with depth
SP15_scale = 20;   LP15_scale = 20;    % phyto
SZ15_scale = 20;   LZ15_scale = 20;    % zoo
A_SP_15 = 20;      A_LP_15 = 30;       % amplitudes
A_SZ_15 = 8;       A_LZ_15 = 15;

% 2100: DCMs (your parameters)
A_SP_21 = A_SP_15 * .1;   mu_SP_21 = 45;  sig_SP_21 = 10;
A_LP_21 = A_LP_15 * 0.02; mu_LP_21 = 45;  sig_LP_21 = 10;
A_SZ_21 = A_SZ_15 * 0.1;  mu_SZ_21 = 45;  sig_SZ_21 = 10;
A_LZ_21 = A_LZ_15 * 0.02; mu_LZ_21 = 45;  sig_LZ_21 = 10;

% Styling (preserved)
col_SP = [0.30 0.65 0.35];
col_LP = [0.10 0.40 0.15];
col_SZ = [0.55 0.40 0.70];
col_LZ = [0.35 0.20 0.55];
lw     = 4.0; fs = 15;

%% ---------------- Build base profiles (preserved) ----------------
% 2015 biomass: monotonic decays
SP2015 = A_SP_15 * exp_decay(z, SP15_scale);
LP2015 = A_LP_15 * exp_decay(z, LP15_scale);
SZ2015 = A_SZ_15 * exp_decay(z, SZ15_scale);
LZ2015 = A_LZ_15 * exp_decay(z, LZ15_scale);

% Add Mixed layer (preserved smoothing approach)
smkn = 30;
SP2015(1:MLD_2015) = SP2015(MLD_2015); SP2015 = movmean(SP2015, smkn);
LP2015(1:MLD_2015) = LP2015(MLD_2015); LP2015 = movmean(LP2015, smkn);
SZ2015(1:MLD_2015) = SZ2015(MLD_2015); SZ2015 = movmean(SZ2015, smkn);
LZ2015(1:MLD_2015) = LZ2015(MLD_2015); LZ2015 = movmean(LZ2015, smkn);

% 2100 biomass: Gaussian DCMs
SP2100 = A_SP_21 * gaussz(z, mu_SP_21, sig_SP_21);
LP2100 = A_LP_21 * gaussz(z, mu_LP_21, sig_LP_21);
SZ2100 = A_SZ_21 * gaussz(z, mu_SZ_21, sig_SZ_21);
LZ2100 = A_LZ_21 * gaussz(z, mu_LZ_21, sig_LZ_21);

% Add Mixed layer (preserved)
SP2100(1:MLD_2100) = SP2100(MLD_2100); SP2100 = movmean(SP2100, smkn);
LP2100(1:MLD_2100) = LP2100(MLD_2100); LP2100 = movmean(LP2100, smkn);
SZ2100(1:MLD_2100) = SZ2100(MLD_2100); SZ2100 = movmean(SZ2100, smkn);
LZ2100(1:MLD_2100) = LZ2100(MLD_2100); LZ2100 = movmean(LZ2100, smkn);

% Depth-integrated totals & overall small fractions (preserved)
TotP_15 = trapz(z, SP2015+LP2015);   TotP_21 = trapz(z, SP2100+LP2100);
TotZ_15 = trapz(z, SZ2015+LZ2015);   TotZ_21 = trapz(z, SZ2100+LZ2100);
fracSP_15 = trapz(z, SP2015)/TotP_15; fracSP_21 = trapz(z, SP2100)/TotP_21;
fracSZ_15 = trapz(z, SZ2015)/TotZ_15; fracSZ_21 = trapz(z, SZ2100)/TotZ_21;

% Depth-wise helper fractions & totals
P15_tot_z = SP2015 + LP2015;   P21_tot_z = SP2100 + LP2100;
Z15_tot_z = SZ2015 + LZ2015;   Z21_tot_z = SZ2100 + LZ2100;

fSP15_z = safediv(SP2015, P15_tot_z, fracSP_15);      % 2015 small-P fraction by depth
fSP21_z = safediv(SP2100, P21_tot_z, fracSP_21);      % 2100 small-P fraction by depth
fSZ15_z = safediv(SZ2015, Z15_tot_z, fracSZ_15);      % 2015 small-Z fraction by depth
fSZ21_z = safediv(SZ2100, Z21_tot_z, fracSZ_21);      % 2100 small-Z fraction by depth

% Normalized total shapes (for Panel 6)
wP21_z = safediv(P21_tot_z, trapz(z, P21_tot_z), 0);  % integrates to 1
wZ21_z = safediv(Z21_tot_z, trapz(z, Z21_tot_z), 0);  % integrates to 1

%% ---------------- Construct diagnostic panels ----------------
% Panel 3: Scale 2015 PHYTO totals to 2100 total PHY, zoo unchanged
sP = TotP_21 / TotP_15;
SP_p3 = sP * SP2015;  LP_p3 = sP * LP2015;
SZ_p3 = SZ2015;       LZ_p3 = LZ2015;

% Panel 4: Scale 2015 ZOO totals to 2100 total ZOO, phyto unchanged
sZ = TotZ_21 / TotZ_15;
SP_p4 = SP2015;       LP_p4 = LP2015;
SZ_p4 = sZ * SZ2015;  LZ_p4 = sZ * LZ2015;

% Panel 5: Keep 2015 Total-P,Z SHAPES & totals; apply 2100 small/large FRACTIONS by depth
SP_p5 = fSP21_z .* P15_tot_z;           LP_p5 = (1 - fSP21_z) .* P15_tot_z;
SZ_p5 = fSZ21_z .* Z15_tot_z;           LZ_p5 = (1 - fSZ21_z) .* Z15_tot_z;

% Panel 6: Keep 2015 totals & 2015 small/large RATIOS; use 2100 Total-P,Z SHAPES
P_tot_p6 = TotP_15 * wP21_z;            Z_tot_p6 = TotZ_15 * wZ21_z;
SP_p6 = fSP15_z .* P_tot_p6;            LP_p6 = (1 - fSP15_z) .* P_tot_p6;
SZ_p6 = fSZ15_z .* Z_tot_p6;            LZ_p6 = (1 - fSZ15_z) .* Z_tot_p6;

%% ---------------- Shared axis bounds (same for all 6 panels) ----------------
all_vals = [SP2015; LP2015; SZ2015; LZ2015; ...
            SP2100; LP2100; SZ2100; LZ2100; ...
            SP_p3;  LP_p3;  SZ_p3;  LZ_p3;  ...
            SP_p4;  LP_p4;  SZ_p4;  LZ_p4;  ...
            SP_p5;  LP_p5;  SZ_p5;  LZ_p5;  ...
            SP_p6;  LP_p6;  SZ_p6;  LZ_p6];
xmax_all = max(all_vals(:));
xlims = [0, 1.05 * xmax_all];
ylims = [0, 150];

%% ---------------- Layout: 1 x 6 ----------------
f = figure('Units','normalized','Position',[0 0 1 .3], 'Color','w');
t = tiledlayout(f, 1, 6, 'Padding','compact','TileSpacing','compact');

% Utility handle for repetitive plotting (styles preserved)
plot_panel = @(ax, SP, LP, SZ, LZ, mld, ttl) ...
    local_plot_panel(ax, z, SP, LP, SZ, LZ, mld, ttl, ...
        col_SP, col_LP, col_SZ, col_LZ, lw, fs, xlims, ylims, ydir_rev);

annotate = @(ax, SP, LP, SZ, LZ) local_annotate(ax, z, SP, LP, SZ, LZ);

%% ---- Panel 1: 2015 (as is) ----
ax1 = nexttile(1); 
plot_panel(ax1, SP2015, LP2015, SZ2015, LZ2015, MLD_2015, '$\overline{2015-2020}$ Profiles');
annotate(ax1, SP2015, LP2015, SZ2015, LZ2015);

%% ---- Panel 2: 2100 (as is) ----
ax2 = nexttile(2);
plot_panel(ax2, SP2100, LP2100, SZ2100, LZ2100, MLD_2100, '$\overline{2095-2100}$ Profiles');
annotate(ax2, SP2100, LP2100, SZ2100, LZ2100);

%% ---- Panel 3: 2015 with Total Phyto = 2100 Total Phyto ----
ax3 = nexttile(3);
plot_panel(ax3, SP_p3, LP_p3, SZ_p3, LZ_p3, MLD_2015, '$\Delta \Sigma$ Phy');
annotate(ax3, SP_p3, LP_p3, SZ_p3, LZ_p3);

%% ---- Panel 4: 2015 with Total Zoo = 2100 Total Zoo ----
ax4 = nexttile(4);
plot_panel(ax4, SP_p4, LP_p4, SZ_p4, LZ_p4, MLD_2015, '$\Delta \Sigma$ Zoo');
annotate(ax4, SP_p4, LP_p4, SZ_p4, LZ_p4);

%% ---- Panel 5: 2015 totals & shapes; 2100 Small/Large fractions by depth ----
ax5 = nexttile(5);
plot_panel(ax5, SP_p5, LP_p5, SZ_p5, LZ_p5, MLD_2015, '$\Delta$ Vertical Distrubution');
annotate(ax5, SP_p5, LP_p5, SZ_p5, LZ_p5);

%% ---- Panel 6: 2015 totals & S:L ratios; 2100 total vertical SHAPE ----
ax6 = nexttile(6);
plot_panel(ax6, SP_p6, LP_p6, SZ_p6, LZ_p6, MLD_2015, '$\Delta$ Community Compositiion');
annotate(ax6, SP_p6, LP_p6, SZ_p6, LZ_p6);

%% ---------------- Helpers (preserved + new local utilities) ----------------
function g = gaussz(z, mu, sigma)
    g = exp(-0.5*((z - mu)./sigma).^2);
    g(g<0) = 0;
end

function g = exp_decay(z, zscale)
    % Simple positive, monotonic decay with depth
    g = exp(-z ./ zscale);
end

function set_common_axes(ax, fs, yrev)
    set(ax,'FontSize',fs,'TickLabelInterpreter','latex');
    if yrev, set(ax,'YDir','reverse'); end
    grid(ax,'on');
end

function r = safediv(a, b, defaultVal)
    % Safe division that supports scalar or array denominators.
    % - If b is scalar and valid, returns a./b.
    % - If b is scalar and invalid (~finite or ~>eps), returns defaultVal (broadcast).
    % - If b is array, divides element-wise where valid; elsewhere defaultVal.

    if isscalar(b)
        if isfinite(b) && abs(b) > eps
            r = a ./ b;                         % broadcast scalar
        else
            r = defaultVal + zeros(size(a));    % fill with default
        end
        return
    end

    % Array denominator path (sizes should match a)
    r = zeros(size(a));
    mask = isfinite(a) & isfinite(b) & (abs(b) > eps);
    r(mask)  = a(mask) ./ b(mask);
    r(~mask) = defaultVal;                      % defaultVal should be scalar
end


function local_plot_panel(ax, z, SP, LP, SZ, LZ, MLD, ttl, ...
            colSP, colLP, colSZ, colLZ, lw, fs, xlims, ylims, yrev)
    axes(ax); hold on
    % Preserve your line styles: small groups dashed, large groups solid
    plot(SP, z, 'Color',colSP, 'LineWidth',lw, 'LineStyle','--');
    plot(LP, z, 'Color',colLP, 'LineWidth',lw, 'LineStyle','-');
    plot(SZ, z, 'Color',colSZ, 'LineWidth',lw, 'LineStyle','--');
    plot(LZ, z, 'Color',colLZ, 'LineWidth',lw, 'LineStyle','-');

    yline(ax, MLD, '-', 'Color',[0.3 0.3 0.3], 'LineWidth',1.2)

    set_common_axes(ax, fs, yrev)
    xlabel(ax,'Biomass (mmolC m$^{-3}$)','Interpreter','latex')
    ylabel(ax,'Depth (m)','Interpreter','latex')
    title(ax, ttl, 'Interpreter','latex','FontWeight','bold','FontSize',18)

    xlim(ax, xlims); ylim(ax, ylims);

    % (Legend intentionally omitted to preserve your current layout)
end

function local_annotate(ax, z, SP, LP, SZ, LZ)
    % Panel-specific totals and overall small fractions (LaTeX-safe text)
    TotP = trapz(z, SP + LP);
    TotZ = trapz(z, SZ + LZ);
    fSP  = trapz(z, SP) / max(eps, TotP);
    fSZ  = trapz(z, SZ) / max(eps, TotZ);

    txt1 = strcat(['$\Sigma$ Phy = ', num2str(round(TotP))]);
    txt2 = strcat(['$\Sigma$ Zoo = ', num2str(round(TotZ))]);
    txt3 = strcat(['$f_{\textrm{SP}} $= ', num2str(round(100*fSP))]);
    txt4 = strcat(['$f_{\textrm{SZ}} $ = ', num2str(round(100*fSZ))]);

    text(ax, 0.965, 0.015, {txt1, txt2, txt3, txt4}, ...
        'Units','normalized', ...
        'Interpreter','latex', ...
        'Color',[0.1 0.1 0.1], ...
        'FontSize',15, ...
        'VerticalAlignment','bottom', ...
        'HorizontalAlignment','right');
end
