%%%%%%%%%%%%%%%%%%%%%%%%%
%% Fig. 5b — Ensemble Maps of loss rate Driver Terms
%%%%%%%%%%%%%%%%%%%%%%%%%
cd /Users/twrohr/Desktop/Matlab/Matlab/Grazing/CMIP6_Comp/NCC_Figs/Final_Figs
addpath(genpath('/Users/twrohr/Desktop/Matlab/Matlab/Utilities'))

for startup =1  
% Load regridded lat lon
cd /Users/twrohr/Desktop/data.nosync/SSP585
         load('Longhurst.mat');
          Lon =  LonGrid; Lat = LatGrid;

% Constants
Ra = 6.371e6;                  % Earth radius [m]
deg2rad = pi / 180;
dphi = deg2rad;               % 1 degree in radians
dlambda = deg2rad;

% Latitude vector (center of cells)
lat = linspace(-89.5, 89.5, 180);  % 180 latitudes from -89.5 to +89.5
lat_rad = lat * deg2rad;

% Area per 1x1 degree cell as a function of latitude
area_per_lat = Ra^2 * dphi * dlambda * cos(lat_rad);  % [m^2]

% Repeat across 360 longitudes to make full 180x360 grid
area_map = repmat(area_per_lat', 1, 360);  % [180x360]

% Convert lon for warppign --- Shift longitudes to [0, 360]
        lonvec = linspace(-179.5, 179.5, 360);
        latvec = linspace(-89.5, 89.5, 180);
        lonvec_wrapped = lonvec;
        lonvec_wrapped(lonvec_wrapped < 0) = lonvec_wrapped(lonvec_wrapped < 0) + 360;
        [lonvec_sorted, lon_order] = sort(lonvec_wrapped);

% Colormaps
    % Get full bam colormap (e.g., 256 colors)
    n_colors = 256;

    % --- Define BU Map  ---
    blue_rgb   = [0.35, 0.55, 0.85];  % muted slate blue
    white_rgb  = [1.0,  1.0,  1.0];   % white
    green_rgb  = [0.30, 0.65, 0.35];  % classy muted green
        % --- Interpolation ---
        n_half = floor(n_colors/2);
        blue_to_white = [linspace(blue_rgb(1), white_rgb(1), n_half)', ...
                         linspace(blue_rgb(2), white_rgb(2), n_half)', ...
                         linspace(blue_rgb(3), white_rgb(3), n_half)'];
        white_to_green = [linspace(white_rgb(1), green_rgb(1), n_colors - n_half)', ...
                          linspace(white_rgb(2), green_rgb(2), n_colors - n_half)', ...
                          linspace(white_rgb(3), green_rgb(3), n_colors - n_half)'];
        % --- Combine into full colormap ---
        cmapBU = [blue_to_white; white_to_green];
    
    % --- Define TD map  ---
    red_rgb     = [0.78, 0.25, 0.25];  % soft classy red
    mustard_rgb = [0.90, 0.7, 0.25];   % strong yellow
        % --- Interpolation ---
        red_to_white = [linspace(red_rgb(1), white_rgb(1), n_half)', ...
                        linspace(red_rgb(2), white_rgb(2), n_half)', ...
                        linspace(red_rgb(3), white_rgb(3), n_half)'];
        white_to_mustard = [linspace(white_rgb(1), mustard_rgb(1), n_colors - n_half)', ...
                            linspace(white_rgb(2), mustard_rgb(2), n_colors - n_half)', ...
                            linspace(white_rgb(3), mustard_rgb(3), n_colors - n_half)'];
        % --- Combine into full colormap ---
        cmapTD = [red_to_white; white_to_mustard];

    %%% Discrete
    green_start = [0.40, 0.60, 0.40];
    green_end   = [0.85, 0.95, 0.85];
    greens = [linspace(green_start(1), green_end(1), 6)', ...
              linspace(green_start(2), green_end(2), 6)', ...
              linspace(green_start(3), green_end(3), 6)'];
    red_start = [0.90, 0.70, 0.70];
    red_end   = [0.70, 0.00, 0.00];
    reds = [linspace(red_start(1), red_end(1), 6)', ...
            linspace(red_start(2), red_end(2), 6)', ...
            linspace(red_start(3), red_end(3), 6)'];
    custom_cmap = [greens; reds];

    % Model Colors
    model_names = {'iHAMMOCC','CMOC','WOMBAT','OECOv2','MARBL','CanOE','MEDUSA2.1',{'PISCESv2';'IPSL'},{'PISCESv2';'CNRM'},'BFM5.2','COBALTv2'};
    mod_order=[10,5,1,9,2,6,11,8,4,3,7,12];
    cmap1 = brewermap(13,'Oranges');
    cmap2 = brewermap(15,'Greens');
    cmap3 = brewermap(11,'Purples');
    cmap_mod = [cmap2([5,7,11,13,15],:); cmap1([5,7,11,13],:); cmap3([7,11],:); 0 0 0];

    CoreVar_Names = {'Phytoplankton Biomass',     'Zooplankton Biomass',         'NPP',                        'Grazing Losses',            'EP',                 ...
                     'Spec. Division Rate',       'Spec. Grazing Rate',          'GSP',                        'Z Loss rate',               'Non-Grazing Losses', ...
                     'Spec. Grazing Loss Rate',   'Fraction Small P',            'Fraction Small Z',           'SST' ,                      'MLD' ,               ...
                     'Spec. Loss Rate',           'Spec. Non-grazing Loss Rate'  'Export Efficeincy'           'Ecosystem Efficeincy'       'Surface [Phyto]' ,  ...
                     'DBM [Phyto]',               'DBM Depth',                   'Dilution Index',             'Diaz Frac',                 'Diatom Frac',           'Non- Diatom Frac'...
                     'Temperature (g)',           'Vertical Dist (g) ',          'Community Composition (g)',  'Phyto. Biomass (g)',        'Zooplankton Biomass (g)',...
                     'Temperature (n)',           'Vertical Dist (n) ',          'Community Composition (n)',  'Phyto. Biomass (n)',        'Zooplankton Biomass (n)',...
                     'Temperature (T)',           'Vertical Dist (T) ',          'Community Composition (T)',  'Phyto. Biomass (T)',        'Zooplankton Biomass (T)',...
                     'Dilution Index',};

     decomp_names = {{'$\Delta l_{T}$'; 'First-order Effect of' ; 'Temperture'}, {'$\Delta l_{VD}$'; 'First-order Effect of' ; 'Vertical Distribution'},...
                     {'$\Delta l_{CC}$';'First-order Effect of' ; 'Community Composition'}, {'$\Delta l_{\Sigma P}$';'First-order Effect of' ; 'Phytoplankton Biomass'},...
                     {'$\Delta l_{\Sigma Z}$';'First-order Effect of';'Zooplankton Biomass'}, {'$\Delta l_{I}$';'First-order Effect of';'Interaction Effects'}};
     decomp_symbols = {{'$\Delta l_{T}$ (d$^{-1}$)'}, {'$\Delta l_{VD}$ (d$^{-1}$)'},...
                     {'$\Delta l_{CC}$ (d$^{-1}$)'}, {'$\Delta l_{\Sigma P}$ (d$^{-1}$)'},...
                     {'$\Delta l_{\Sigma Z}$ (d$^{-1}$)'},{'$\Delta l_{Int}$ (d$^{-1}$)'}};
    
     decomp_names2 =  { '$\frac{\Delta \; l_{T}}{\Delta\;  \mathrm{SST}}$' , '$\frac{\Delta \; l_{VD} }{\Delta \; \mathrm{D}_L}$' ,...
                    '$\frac{\Delta \; l_{CC} }{\Delta\;  f_{\mathrm{SP}}}$' , '$\frac{\Delta \; l_{\Sigma P}}{\Delta \; \Sigma\,{\mathrm{P}}}$' ,...
                     '$\frac{\Delta \; l_{\Sigma Z} }{\Delta\;  \Sigma \mathrm{P} }$' };

     decomp_units = {'(d$^{-1}$ $^{\circ}$C warmer $^{-1}$)', '(d$^{-1}$ m $^{-1}$)', '(d$^{-1}$ $\%$ inc. small phy. $^{-1}$)', '(m$^{2}$ d$^{-1}$ mmol$^{-1}$)', '(m$^{2}$ d$^{-1}$ mmol$^{-1}$)'};

     Driver_label = {'$\Delta$ SST ($^\circ$C)',' $\Delta$ Dilution Length Scale ($D_{L}$; m)', ' $\Delta$ Fraction Sm. Phytoplankton ($f_{sp}$; $\%$)', '$\Delta$ $\Sigma$ Phytoplankton (mmolC m$^{-2}$)',...
                     '$\Delta$ $\Sigma$ Zooplankton (mmolC m$^{-2}$)','dummy'};

     % decomp_names =decomp_names([4,5,1,3,2]);
     % decomp_names2 =decomp_names2([4,5,1,3,2]);
     % Driver_label =Driver_label([4,5,1,3,2]);
     % decomp_units =decomp_units([4,5,1,3,2]);
     % decomp_symbols =decomp_symbols([4,5,1,3,2]); 
     CoreVar_Names{37} = '-';

     wedge_colors = [
    0.90 0.55 0.25;   % SST     → soft orange
    0.35 0.60 0.85;   % MLD     → clearer blue
    0.25 0.75 0.65;   % Diatoms → more greenish teal
    0.60 0.75 0.45;   % Phyto   → sage green
    0.85 0.55 0.75    % Zoo     → soft magenta
];
lighter_colors = wedge_colors + 0.6 .* (1 - wedge_colors);  % shift toward white, preserve hue

end

%% === Loss Rate =========================
% Note, this plots 5 variables - need to run twice to get full set of 6 terms plotted in 5b - changing the 'vord' in line 142 to include '6' then make by hand.

% =========================================
clc; close all
clearvars plot

% ---------------- CONFIG ----------------
% Loss-rate variables (primary)
vord        =  [4,5,1,3,2];  %  Note, this plots 5 variables - need to run twice to get full set of 6 terms plotted in 5b - changing the 'vord' in line 142 to include '6' then make by hand.        
%    SI: [4,5,1,3,2] | Fig 4: [4,2,1,3,2] L [5,3,6,3,2]] R. 
VARIDS      = [37,38,39,40,41,47]; 
ENVIDS      = [14,42,12,1,  2, 1];
var_ids     = VARIDS(vord);        
env_drivers = ENVIDS(vord);
vars2plot   = 1:5;
zonal_mean   = 1 ;% 1 for mean, 2 for median 

% Driver variables for new left column

%env_drivers = [1, 1, 1, 1, 1];

mods2plot  = 1:11;
model_id   = 1:11;
nbins      = 100;                  % vertical bins for violins
bgRed      = [0.85 0.30 0.30];     % muted red (y<0 band)
bgBlue     = [0.30 0.50 0.85];     % muted blue (y>0 band)
bgAlpha    = 0.06;                 % background transparency
lineThin   = 0.5;                  % per-model zonal curve width
lineMean   = 2.2;                  % ensemble mean width

% Dynamic latitude window — set once
minLat     = -65;        % <<< adjust as desired
maxLat     =  65;        % <<< adjust as desired
tick_step  = 20;         % for latitude labels

% Column width fractions (among remaining width after map)
zonal_frac  = 0.20;      % each zonal column fraction
% violin width will be: (1 - 2*zonal_frac) of the remaining width (so it’s reduced by one zonal)
map_aspect = 1.85;       % map_w ≈ map_aspect * map_h

% ---------------- DATA ----------------
basedir    = '/Users/twrohr/Desktop/data.nosync/SSP585';
regrid_dir = fullfile(basedir,'Abridged/regridded2');
load(fullfile(basedir,'All_Stats_Raw.mat'));   % expects: InterMod_mean_maps, Lat_t, Lon_t, latvec, zon_stats, cmap_mod, etc.

% Resolve names
NM           = resolve_names_from_stats(zon_stats);
modelNames   = NM.modelNames(model_id);
modelFields  = NM.modelFields(model_id);
cmapModels   = cmap_mod(model_id,:);

% Regridded *_rw files (for violins) in model order
mods0      = dir(fullfile(regrid_dir,'*.mat'));
mods0      = mods0(~[mods0.isdir]);
mod_order  = [10,5,1,9,2,6,11,8,4,3,7];  % provided order
filesOrdered = arrayfun(@(k) fullfile(regrid_dir, mods0(mod_order(k)).name), 1:numel(mod_order), 'uni',0);

% ---------------- LAYOUT ----------------
nrows    = numel(var_ids);
left_m   = 0.02; right_m = .02; top_m = 0.02; bot_m = 0.06;
col_gap  = 0.025;                        % gap between columns
row_gap  = 0.003;

row_h    = (1 - top_m - bot_m - (nrows-1)*row_gap) / nrows;   % target row height
map_h    = row_h;                          % fixed map height per spec
map_w    = map_aspect * map_h;             % width scaled from height

% Remaining width after placing the map + THREE gaps (col1-2, 2-3, 3-4)
rem_w    = 1 - left_m - right_m - map_w - 3*col_gap;
if rem_w <= 0
    error('Map width too large for margins/gaps. Reduce map_aspect or col_gap.');
end
zonal_w  = rem_w * zonal_frac;         % each zonal column
violin_w = rem_w - 1.3*zonal_w ;          % reduced by exactly one zonal width

figure('Units','normalized','Position',[0 0 1 1], 'Color','w');

% =========================================================
% ---------------- LOOP OVER VARIABLES (ROWS) -------------
% =========================================================
for r = vars2plot
    vidxs   = var_ids(r)-10;                 % loss-rate decomposition index
    vidx2  = vidxs + 5;                   % paired non-grazing index
    vidx3  = vidxs-27+1;                  % grazing sens index
    vidx4  = vidxs-27+5+1;                % non-grazing sens index
    pidx   = var_ids(r);
    % if any(vidxs ~= var_ids) ; 
    % else; pidx =vidxs +10; 
    % end
    vField = NM.varFields{vidxs};
    vLabel = NM.varNames{vidxs};

    vField2 = ''; vLabel2 = '';
    if vidx2 <= numel(NM.varFields)
        vField2 = NM.varFields{vidx2};
        vLabel2 = NM.varNames{vidx2};
    end

    % --- row vertical position ---
    topEdge = 1 - top_m - (r-1)*(row_h + row_gap);
    y0_row  = topEdge - row_h;

    % Precompute some x positions for the 4 columns
     x3_map     = left_m -.05;                          % Col 1: Map
     x4_lossZ   = x3_map     + map_w   + col_gap -.08;  % Col 2: Loss-rate Zonal
     x1_drivers = x4_lossZ   + zonal_w + col_gap  ; % Col 3: Drivers Zonal
     x2_violin  = x1_drivers + zonal_w + col_gap +.007;  % Col 4: Sesnitivity Violin
   

% =====================================================
% (1) MAP (First column) — create first to get projection & y-range
% =====================================================
    pos_map = [x3_map, y0_row, map_w, map_h];
    axm = axes('Units','normalized','Position',pos_map, ...
               'PositionConstraint','innerposition');
    axes(axm); axesm robinson; axis off; tightmap
    axesm('robinson','MapLonLimit',[0 360],'MapLatLimit',[minLat maxLat],'Origin',[0 180 0]);

    % Inter-model mean for this variable (TOTALS: vidx+10), mask to lat window
    plotvar = squeeze(InterMod_mean_maps(pidx,:,:)); % 
    lat_rows = Lat_t(:,1);
    mask_lat = (lat_rows < minLat) | (lat_rows > maxLat);
    plotvar(mask_lat,:) = NaN;

    Ra = georefcells([min(Lat_t(:,1)) max(Lat_t(:,1))], ...
                     [min(Lon_t(1,:)) max(Lon_t(1,:))], size(plotvar));
    geoshow(plotvar, Ra, 'DisplayType','texturemap');
    load coastlines; patchm(coastlat,coastlon,[.7 .7 .7]);

    % colors (reversed balance)
    % Use fixed or auto-centered limits as before
    climv = [ -.062, .062];
    colormap(axm, flipud(cmocean('balance')));
    clim(climv);

    % Annotations
    drawnow;                                
    map_pos = get(axm,'Position');          % [x y w h]
    lbl = decomp_names{vord(r)};
    annotation('textbox', [map_pos(1)+.065   map_pos(2)+.002 .1 .05], ...
        'String', lbl,'Interpreter', 'latex', ...
        'HorizontalAlignment', 'center', 'VerticalAlignment', 'middle', ...
        'Rotation', 90, 'FontSize', 14, 'EdgeColor', 'none');

    % Single colorbar on last row (match primary)
    cb_h    = 0.015; % .02
    cb_gap  = 0.023;
    if r==5
        cb = colorbar(axm,'southoutside');
        cb.Units = 'normalized';
        cb.TickLabelInterpreter = 'latex';
        cb.FontSize = 11;
        %cb.Position = [map_pos(1)+0.15*map_pos(3), map_pos(2)-cb_gap-cb_h, 0.6*map_pos(3), cb_h]; % [map_pos(1)+0.31*map_pos(3),  map_pos(2)-cb_gap-cb_h,   0.7*map_pos(3), cb_h]
        cb.Position = [map_pos(1)+0.2*map_pos(3),  map_pos(2)-cb_gap-cb_h,   0.6*map_pos(3), cb_h]
       
        pos = cb.Position; 
         annotation('textbox', [pos(1)-.0 pos(2)+pos(4)+.0 pos(3) 0.02], ...
               'String','(d$^{-1}$)', 'Interpreter','latex', ...
               'HorizontalAlignment','center', 'VerticalAlignment','bottom', ...
               'EdgeColor','none', 'FontSize', 15);
    end

    % Projection helpers for consistent y-range & ticks
    mstruct = gcm(axm);
    [~, y_minProj] = mfwdtran(mstruct, minLat, 0);
    [~, y_maxProj] = mfwdtran(mstruct, maxLat, 0);
    ybot = min(y_minProj, y_maxProj);
    ytop = max(y_minProj, y_maxProj);
    yl      = get(axm,'YLim');
    toNorm  = @(y) map_pos(2) + map_pos(4) * (y - yl(1)) / (yl(2)-yl(1));
    ybot_n  = toNorm(ybot);   
    ytop_n  = toNorm(ytop);

    % Shared tick-lats & projected positions
    tick_start = ceil(minLat / tick_step) * tick_step;
    tick_end   = floor(maxLat / tick_step) * tick_step;
    ticklats   = (tick_start:tick_step:tick_end)';
    [~, ytick_map] = mfwdtran(mstruct, ticklats, zeros(size(ticklats)));


end % rows

% =========================
% ----- Helpers ----------
% =========================
function out = ternary(cond, a, b)
if cond, out = a; else, out = b; end
end

function counts = histcounts_weighted(x, edges, w)
% Weighted histogram per bin edge
    x = x(:); w = w(:);
    idx   = discretize(x, edges);
    ok    = ~isnan(idx) & isfinite(w);
    counts_col = accumarray(idx(ok), w(ok), [numel(edges)-1, 1], @sum, 0);
    counts = counts_col.';   % row
end

function [q50,q25,q75] = local_wquant3(x,w)
% Robust weighted 50/25/75 percentiles
    ok = isfinite(x) & isfinite(w) & (w > 0);
    x  = x(ok); w = w(ok);
    if isempty(x) || sum(w)<=0, q50=NaN; q25=NaN; q75=NaN; return, end
    [xs,ix] = sort(x(:));
    ws = w(:); ws = ws(ix);
    cw = cumsum(ws); cw = cw./cw(end);
    [cu,iu] = unique(cw,'stable'); xsu = xs(iu);
    if numel(cu)==1, q50=xsu; q25=xsu; q75=xsu; return, end
    q25 = interp1(cu, xsu, 0.25, 'linear','extrap');
    q50 = interp1(cu, xsu, 0.50, 'linear','extrap');
    q75 = interp1(cu, xsu, 0.75, 'linear','extrap');
end

function out = smartfmt(v)
% Always scientific notation, 2 significant figures, 2 lines
    if ~isfinite(v) || v == 0
        mant = 0;  exp10 = 0;
    else
        exp10 = floor(log10(abs(v)));
        mant  = v / 10^exp10;
        mant  = round(mant, 2, 'significant');
        if abs(mant) >= 10
            mant  = mant / 10;
            exp10 = exp10 + 1;
        end
    end
    mant_str = sprintf('%+.1f', mant);
    if mant == 0 && exp10 ==0
        line1 = ''; line2 = '';
    else
        line1 = ['$' mant_str '$'];
        line2 = ['e ' num2str(exp10)];
    end
    out = {line1; line2};
end

function c = signcolor(v)
% Blue if +, red if -, gray if ~0/NaN
    if ~isfinite(v) || v==0
        c = [0.35 0.35 0.35];
    elseif v > 0
        c = [0.10 0.45 0.90];
    else
        c = [0.80 0.20 0.20];
    end
end

function draw_half_box(ax, c, side, med, q25, q75, maxW, col)
    if ~all(isfinite([med,q25,q75])) || q75<=q25, return; end
    w  = 0.55*maxW;
    x0 = c + strcmp(side,'left')*(-w) + strcmp(side,'right')*(0);
    r  = rectangle(ax,'Position',[x0, q25, w, (q75-q25)], ...
        'FaceColor',[0 0 0 0.12],'EdgeColor',col,'LineWidth',0.8);
    m  = plot(ax, [x0, x0+w], [med, med], 'Color', col, 'LineWidth', 1.6);
    uistack([r m],'top');
end

function f = weighted_frac(x, w, cond)
    valid = isfinite(x) & isfinite(w) & (w > 0);
    x    = x(valid);
    w    = w(valid);
    cond = cond(valid);

    den = sum(w);
    if den == 0
        f = NaN;
    else
        f = sum(w(cond)) / den;
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%% End — Combined Final Script %%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
