%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Figure 1e Insets - schametics %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cd /Users/twrohr/Desktop/Matlab/Matlab/Grazing/CMIP6_Comp/NCC_Figs/Final_Figs/Revised/

% case_to_plot controls the perturbation plotted against case 1:
%
%   2 = Coupled increase     (purple)
%   3 = Decoupled decrease   (green)
%   4 = Coupled decrease     (orange)
%   5 = Decoupled increase   (blue)
%
% All three main panels always show:
%   - Case 1 in black
%   - The selected case in its assigned colour
%
% The inset in the top-right of the biomass panel shows annual-mean
% percentage changes in mu, loss, and phytoplankton biomass relative
% to case 1.

clear; close all; clc

%% ========================================================================
% User controls
% =========================================================================

% Choose case 2, 3, 4, or 5
case_to_plot = 2;

if ~ismember(case_to_plot, 2:5)
    error('case_to_plot must be 2, 3, 4, or 5.')
end

% ===== Baseline system: Case 1 =====
mu_amp1      = 0.85;
loss_amp1    = 0.85;

mu_center1   = 0;
loss_lag1    = 1.5;
loss_center1 = mu_center1 + loss_lag1;

sigma = 1.55;

% ===== Case 2: Coupled increase =====
amp_factor2 = 1.5;
loss_lag2   = 1.5;

% ===== Case 3: Decoupled decrease =====
amp_factor3 = 1.5;
loss_lag3   = 0.5;

% ===== Case 4: Coupled decrease =====
amp_factor4 = 0.5;
loss_lag4   = 1.5;

% ===== Case 5: Decoupled increase =====
amp_factor5 = 0.7;
loss_lag5   = 3.0;

% Domain control
tail_widths = 3.4;

% Biomass initial condition
B0 = 1;

%% ========================================================================
% Plot settings
% =========================================================================

font_size = 14;
label_size = 22;
title_size = 22;
legend_size = 11;

lw_main = 4.4;
lw_zero = 1.2;

% Baseline colour
baseline_rgb = [0.00, 0.00, 0.00];

% Assigned case colours
% Bottom-left: blue branch
blue_1_rgb = [0.70, 0.83, 0.98];
blue_rgb = [0.43, 0.65, 0.90];
blue_3_rgb = [0.22, 0.45, 0.76];
blue_4_rgb = [0.08, 0.24, 0.55];

% Top-right: green branch
green_1_rgb = [0.82, 0.89, 0.57];
green_rgb = [0.60, 0.76, 0.35];
green_3_rgb = [0.38, 0.58, 0.19];
green_4_rgb = [0.18, 0.38, 0.08];

% Top-left: purple branch
purple_1_rgb = [0.86, 0.73, 0.94];
purple_rgb = [0.69, 0.50, 0.82];
purple_3_rgb = [0.51, 0.31, 0.68];
purple_4_rgb = [0.32, 0.12, 0.48];

% Bottom-right: orange branch
orange_1_rgb = [0.98, 0.76, 0.38];
orange_rgb = [0.91, 0.56, 0.16];
orange_3_rgb = [0.75, 0.37, 0.06];
orange_4_rgb = [0.52, 0.20, 0.02];

case_colors = [ ...
    baseline_rgb; ...
    green_rgb; ...
    orange_rgb; ...
    blue_rgb; ...
    purple_rgb];

case_names = { ...
    'Baseline'; ...
    'Coupled increase'; ...
    'Decoupled decrease'; ...
    'Coupled decrease'; ...
    'Decoupled increase'};

selected_color = case_colors(case_to_plot,:);
selected_name  = case_names{case_to_plot};

baseline_legend = 'Case 1: Baseline';
selected_legend = sprintf('Case %d: %s', ...
    case_to_plot, selected_name);

% Growth and loss line styles
mu_line_style   = '-';
loss_line_style = ':';

% Month labels
month_labels = {'J','F','M','A','M','J','J','A','S','O','N','D'};

%% ========================================================================
% Derived scenario parameters
% =========================================================================

n_cases = 5;

% Multipliers relative to baseline amplitudes
amp_factor = [ ...
    1.0, ...
    amp_factor2, ...
    amp_factor3, ...
    amp_factor4, ...
    amp_factor5];

% Loss lags for all cases
loss_lag = [ ...
    loss_lag1, ...
    loss_lag2, ...
    loss_lag3, ...
    loss_lag4, ...
    loss_lag5];

% Growth and loss amplitudes
mu_amp   = mu_amp1   .* amp_factor;
loss_amp = loss_amp1 .* amp_factor;

% Growth timing remains fixed for all cases
mu_center = repmat(mu_center1, 1, n_cases);

% Loss timing differs among cases
loss_center = mu_center + loss_lag;

%% ========================================================================
% Common x-domain
% =========================================================================

all_left_edges = [ ...
    mu_center   - tail_widths .* sigma, ...
    loss_center - tail_widths .* sigma];

all_right_edges = [ ...
    mu_center   + tail_widths .* sigma, ...
    loss_center + tail_widths .* sigma];

left_edge  = min(all_left_edges);
right_edge = max(all_right_edges);

x = linspace(left_edge, right_edge, 1600);

% Twelve evenly spaced month ticks across the annual schematic domain
month_ticks = linspace(left_edge, right_edge, 12);

%% ========================================================================
% Construct all cases
% =========================================================================

mu = zeros(n_cases, numel(x));
loss = zeros(n_cases, numel(x));
r = zeros(n_cases, numel(x));
B = zeros(n_cases, numel(x));

for c = 1:n_cases

    mu(c,:) = mu_amp(c) .* ...
        exp(-0.5 .* ((x - mu_center(c)) ./ sigma).^2);

    loss(c,:) = loss_amp(c) .* ...
        exp(-0.5 .* ((x - loss_center(c)) ./ sigma).^2);

    r(c,:) = mu(c,:) - loss(c,:);

    B(c,:) = B0 .* exp(cumtrapz(x, r(c,:)));

    % Normalize each biomass trajectory to its initial value
    B(c,:) = B(c,:) ./ B(c,1);
end

%% ========================================================================
% Select baseline and perturbation variables
% =========================================================================

% Baseline: Case 1
mu1   = mu(1,:);
loss1 = loss(1,:);
r1    = r(1,:);
B1    = B(1,:);

% Selected case
mu_selected   = mu(case_to_plot,:);
loss_selected = loss(case_to_plot,:);
r_selected    = r(case_to_plot,:);
B_selected    = B(case_to_plot,:);

%% ========================================================================
% Annual-mean changes relative to Case 1
% =========================================================================
%
% Percentage changes are used because biomass and rates have different
% physical units and therefore should not be shown together as absolute
% differences on one common bar axis.

annual_duration = x(end) - x(1);

mean_mu = trapz(x, mu, 2) ./ annual_duration;
mean_loss = trapz(x, loss, 2) ./ annual_duration;
mean_B = trapz(x, B, 2) ./ annual_duration;

delta_mu_pct = 100 .* ...
    (mean_mu(case_to_plot) - mean_mu(1)) ./ mean_mu(1);

delta_loss_pct = 100 .* ...
    (mean_loss(case_to_plot) - mean_loss(1)) ./ mean_loss(1);

delta_B_pct = 100 .* ...
    (mean_B(case_to_plot) - mean_B(1)) ./ mean_B(1);

annual_change_pct = [ ...
    delta_mu_pct, ...
    delta_loss_pct, ...
    delta_B_pct];

%% ========================================================================
% Peak diagnostics
% =========================================================================

[~, i_B_peak] = max(B, [], 2);
x_B_peak = x(i_B_peak);

[~, i_r_max] = max(r, [], 2);
x_r_max = x(i_r_max);

mu_at_r_max = zeros(n_cases,1);

for c = 1:n_cases
    mu_at_r_max(c) = mu(c, i_r_max(c));
end

%% ========================================================================
% Dynamic y-limits based on Case 1 and selected case
% =========================================================================

B_plot_values = [B1(:); B_selected(:)];

ylim_B = [ ...
    min(B_plot_values) .* 0.98, ...
    max(B_plot_values) .* 1.05];

r_abs = max(abs([r1(:); r_selected(:)]));

if r_abs == 0
    r_abs = 1;
end

ylim_r = [-1.08, 1.08] .* r_abs;

rate_max = max([ ...
    mu1(:); ...
    loss1(:); ...
    mu_selected(:); ...
    loss_selected(:)]);

ylim_rate = [0, rate_max .* 1.12];

%% ========================================================================
% Figure
% =========================================================================

fig = figure( ...
    'Color','w', ...
    'Units','centimeters', ...
    'Position',[0.5, 2, 15, 40]);

t = tiledlayout(fig, 3, 1);

t.TileSpacing = 'loose';
t.Padding = 'compact';

%% ========================================================================
% Top panel: phytoplankton biomass
% =========================================================================

ax3 = nexttile(t);
hold(ax3,'on')

plot(ax3, x, B1, ...
    '-', ...
    'Color',baseline_rgb, ...
    'LineWidth',lw_main, ...
    'DisplayName',baseline_legend);

plot(ax3, x, B_selected, ...
    '-', ...
    'Color',selected_color, ...
    'LineWidth',lw_main, ...
    'DisplayName',selected_legend);

box(ax3,'on')

set(ax3, ...
    'FontSize',font_size, ...
    'XTick',month_ticks, ...
    'XTickLabel',month_labels, ...
    'TickDir','in', ...
    'TickLabelInterpreter','latex');

xlim(ax3,[left_edge right_edge])
ylim(ax3,ylim_B)

ylabel(ax3, ...
    'mmol m$^{-3}$', ...
    'FontSize',label_size, ...
    'Interpreter','latex');

title(ax3, ...
    'Phytoplankton Biomass', ...
    'FontSize',title_size, ...
    'Interpreter','latex');


%% ========================================================================
% Middle panel: accumulation rate r = mu - loss
% =========================================================================

ax2 = nexttile(t);
hold(ax2,'on')

plot(ax2, x, r1, ...
    '-', ...
    'Color',baseline_rgb, ...
    'LineWidth',lw_main, ...
    'DisplayName',baseline_legend);

plot(ax2, x, r_selected, ...
    '-', ...
    'Color',selected_color, ...
    'LineWidth',lw_main, ...
    'DisplayName',selected_legend);

yline(ax2, 0, ...
    'k-', ...
    'LineWidth',lw_zero, ...
    'HandleVisibility','off');

box(ax2,'on')

set(ax2, ...
    'FontSize',font_size, ...
    'XTick',month_ticks, ...
    'XTickLabel',month_labels, ...
    'TickDir','in', ...
    'TickLabelInterpreter','latex');

xlim(ax2,[left_edge right_edge])
ylim(ax2,ylim_r)

ylabel(ax2, ...
    '$\mathrm{d}^{-1}$', ...
    'FontSize',label_size, ...
    'Interpreter','latex');

title(ax2, ...
    'Accumulation Rate ($\mu-l$)', ...
    'FontSize',title_size, ...
    'Interpreter','latex');

%% ========================================================================
% Bottom panel: growth and loss rates
% =========================================================================

ax1 = nexttile(t);
hold(ax1,'on')

% Baseline growth
plot(ax1, x, mu1, ...
    mu_line_style, ...
    'Color',baseline_rgb, ...
    'LineWidth',lw_main, ...
    'DisplayName','Case 1 growth ($\mu_1$)');

% Baseline loss
plot(ax1, x, loss1, ...
    loss_line_style, ...
    'Color',baseline_rgb, ...
    'LineWidth',lw_main, ...
    'DisplayName','Case 1 loss ($l_1$)');

% Selected-case growth
plot(ax1, x, mu_selected, ...
    mu_line_style, ...
    'Color',selected_color, ...
    'LineWidth',lw_main, ...
    'DisplayName',sprintf( ...
        'Case %d growth ($\\mu_%d$)', ...
        case_to_plot, case_to_plot));

% Selected-case loss
plot(ax1, x, loss_selected, ...
    loss_line_style, ...
    'Color',selected_color, ...
    'LineWidth',lw_main, ...
    'DisplayName',sprintf( ...
        'Case %d loss ($l_%d$)', ...
        case_to_plot, case_to_plot));

box(ax1,'on')

set(ax1, ...
    'FontSize',font_size, ...
    'XTick',month_ticks, ...
    'XTickLabel',month_labels, ...
    'TickDir','in', ...
    'TickLabelInterpreter','latex');

xlim(ax1,[left_edge right_edge])
ylim(ax1,ylim_rate)

ylabel(ax1, ...
    '$\mathrm{d}^{-1}$', ...
    'FontSize',label_size, ...
    'Interpreter','latex');

xlabel(ax1, ...
    'Month', ...
    'FontSize',label_size, ...
    'Interpreter','latex');

title(ax1, ...
    'Growth (--) and Loss Rates (- -)', ...
    'FontSize',title_size, ...
    'Interpreter','latex');


%% ========================================================================
% Link main-panel x-axes
% =========================================================================

linkaxes([ax1, ax2, ax3], 'x')

%% ========================================================================
% Inset bar chart: annual-mean changes relative to Case 1
% =========================================================================
%
% The inset is created after the tiled layout has been rendered so its
% location can be defined relative to the final position of the top panel.

drawnow

original_units = ax3.Units;
ax3.Units = 'normalized';
top_position = ax3.Position;
ax3.Units = original_units;

% Position within the upper-right region of the top panel
inset_position = [ ...
    top_position(1) + 0.6 .* top_position(3), ...
    top_position(2) + 0.655 .* top_position(4), ...
    0.4 .* top_position(3), ...
    0.34 .* top_position(4)];

ax_inset = axes( ...
    'Parent',fig, ...
    'Units','normalized', ...
    'Position',inset_position);

hold(ax_inset,'on')

bar(ax_inset, ...
    1:3, ...
    annual_change_pct, ...
    0.65, ...
    'FaceColor',selected_color, ...
    'EdgeColor','none');

% Prominent zero-change line
yline(ax_inset, 0, ...
    'k-', ...
    'LineWidth',1.4, ...
    'HandleVisibility','off');

% Use symmetric limits so positive and negative changes are clear
inset_limit = 1.25 .* max(abs(annual_change_pct));

if ~isfinite(inset_limit) || inset_limit < 5
    inset_limit = 5;
end

ylim(ax_inset,[-inset_limit inset_limit])
xlim(ax_inset,[0.4 3.6])

set(ax_inset, ...
    'Color','w', ...
    'Box','off', ...
    'Layer','top', ...
    'LineWidth',0.9, ...
    'FontSize',12, ...
    'TickDir','in', ...
    'XTick',1:3, ...
    'XTickLabel',{ ...
        '$\Delta_{\%} \mu$', ...
        '$\Delta_{\%} l$', ...
        '$\Delta_{\%} \Sigma \textrm{P}$'}, ...
    'TickLabelInterpreter','latex');

% White figure background
set(fig, ...
    'Color', 'w', ...
    'InvertHardcopy', 'off');

% Thin black border around the entire figure
annotation(fig, 'rectangle', ...
    [0.002 0.002 0.996 0.996], ...
    'Color', 'k', ...
    'LineWidth', 0.8, ...
    'FaceColor', 'none');

copygraphics(fig, ...
    'ContentType','vector', ...
    'BackgroundColor','white');

%% ========================================================================
% Diagnostics
% =========================================================================

fprintf('\n============================================================\n')
fprintf('Selected comparison\n')
fprintf('============================================================\n')
fprintf('Case 1 versus Case %d: %s\n', ...
    case_to_plot, selected_name)

fprintf('\nBiomass peak x-locations:\n')
fprintf('  Case 1 baseline : %.3f\n', x_B_peak(1))
fprintf('  Case %d         : %.3f\n', ...
    case_to_plot, x_B_peak(case_to_plot))

fprintf('\nMaximum accumulation-rate x-locations:\n')
fprintf('  Case 1 baseline : x = %.3f, mu = %.3f\n', ...
    x_r_max(1), mu_at_r_max(1))

fprintf('  Case %d         : x = %.3f, mu = %.3f\n', ...
    case_to_plot, ...
    x_r_max(case_to_plot), ...
    mu_at_r_max(case_to_plot))

fprintf('\nAnnual-mean changes relative to Case 1:\n')
fprintf('  Delta mu   : %+.2f %%\n', delta_mu_pct)
fprintf('  Delta loss : %+.2f %%\n', delta_loss_pct)
fprintf('  Delta P    : %+.2f %%\n', delta_B_pct)
fprintf('============================================================\n\n')

%% ========================================================================
% Optional export
% =========================================================================

% pdf_name = sprintf('schematic_change_case%d.pdf', case_to_plot);
% png_name = sprintf('schematic_change_case%d.png', case_to_plot);
%
% exportgraphics(fig, pdf_name, 'ContentType','vector')
% exportgraphics(fig, png_name, 'Resolution',300)
