%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Fig. S4 Plot 2015 Loss Rates
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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 -- 0.45, 0.60, 0.80
    white_rgb  = [1.0,  1.0,  1.0];   % white
    green_rgb  = [0.30, 0.65, 0.35];  % classy muted green -- 0.40, 0.60, 0.40
        % --- 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 -- 0.75, 0.30, 0.30
    mustard_rgb = [0.90, 0.7, 0.25];  % strong end member yellow -- 0.80, 0.65, 0.20
        % --- 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];


   
    %%% Discete
    % Interpolate 6 steps from strong to weak BU
    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)'];
    % Interpolate 6 steps from weak to strong TD
    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)'];
    % Combine
    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)',};


end



%% === Plot One Variable + Right Zonal Panel ===
var_id   = 3;                   % single variable index
model_id = 1:11;                % which models to include

cd /Users/twrohr/Desktop/data.nosync/SSP585
load('All_Stats_2015.mat')      

figure('Units','normalized','Position',[0 0 1 1.5]); set(gcf,'color','w');

% Resolve names from your stats struct that holds zonal results
 NM = resolve_names_from_stats(zon_stats);   % same helper as before
 vidx   = var_id;
 vField = NM.varFields{vidx};
 vLabel = NM.varNames{vidx};
        
map_left  = 0.06;      % map left
map_width = 0.48;      % map width
map_height= 0.40;      % map height
map_top   = 0.95;      % top edge of the map (so y = top - height)
gap_h     = 0.035;     % gap between map and zonal panel
zonal_w   = 0.22;      % <- YOU control the zonal panel width


%%%%%%%%%%%%%%%%%%%%%
%%% --- MAP ----
%%%%%%%%%%%%%%%%%%%%%
for Plot_Map = 1
    map_pos0  = [map_left, map_top-map_height, map_width, map_height];
    axm = axes('Units','normalized','Position',map_pos0, ...
               'PositionConstraint','innerposition');  % colorbar won't resize the map
    axes(axm); axesm robinson; axis off; tightmap
    axesm('robinson','MapLonLimit',[0 360],'MapLatLimit',[-75 75],'Origin',[0 180 0]);

    % Create Map and Raster 
    plotvar = squeeze(InterMod_mean_maps(16,:,:));
    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 (reverse the colormap)
    pr = prctile(plotvar(:),[1 99]);
    cmins = [pr];
    colormap(axm, flipud(cmocean('thermal')));  % <<< reversed
    clim(cmins);

    %title('Change Phytoplankton Specific Loss Rate', 'Interpreter','latex','FontSize',25);

    % -------- Colorbar BELOW the map on a separate strip --------
    drawnow;                       % let tightmap settle the final position
    map_pos = get(axm,'Position'); % [x y w h] normalized

    cb_h  = 0.025;   % colorbar height
    cb_gap= 0.018;   % gap between map and colorbar

    cb = colorbar(axm,'southoutside');   
    cb.Units = 'normalized'; cb.FontSize = 14; cb.TickLabelInterpreter = 'latex';
    cb.Position = [map_pos(1)+.11, map_pos(2)-cb_gap-cb_h, map_pos(3)-.22, cb_h];  % separate axis
    pos = cb.Position;  % [x y w h], normalized
    annotation('textbox', [pos(1) pos(2)+pos(4) pos(3) 0.02], ...
               'String',' $l$ (d$^{-1}$)', 'Interpreter','latex', 'HorizontalAlignment','center', 'VerticalAlignment','bottom', 'EdgeColor','none', 'FontSize', 16);
end


%%%%%%%%%%%%%%%%%%%
%%% Zonal Averages 
%%%%%%%%%%%%%%%%%%%
for Zonal_plot = 1

% --- Frame → exact globe vertical span (map data units) ---
setm(axm,'Frame','on');
hFrame = handlem('Frame');  if isempty(hFrame), hFrame = framem; end
if numel(hFrame) > 1, hFrame = hFrame(1); end
Yf = get(hFrame,'YData'); if iscell(Yf), Yf = Yf{1}; end
ybot = min(Yf);  ytop = max(Yf);

% Map axes position -> zonal axes position (normalized units)
map_pos = get(axm,'Position');     % [x y w h]
yl      = get(axm,'YLim');         % map-data Y
toNorm  = @(y) map_pos(2) + map_pos(4) * (y - yl(1)) / (yl(2)-yl(1));
ybot_n  = toNorm(ybot);  ytop_n = toNorm(ytop);

gap_h   = 0.015;                                  % horizontal gap to the right panel
zonal_x = map_pos(1) + map_pos(3) + gap_h;
zon_pos = [zonal_x, ybot_n, zonal_w, max(ytop_n - ybot_n, 1e-9)];

axz = axes('Units','normalized','Position',zon_pos,'Clipping','on'); 
hold(axz,'on');

% --- Project latitude centers to same Y as the map (for vertical alignment) ---
mstruct  = gcm(axm);
lat_cent = latvec(:);
[~, y_cent] = mfwdtran(mstruct, lat_cent, zeros(size(lat_cent)));

% --- Collect model curves (Nlat x Nmod) from zonal stats (medians) ---
Nlat = numel(lat_cent);
mods = numel(model_id);
Z = nan(Nlat, mods);
Zstd = nan(Nlat, mods);
for ii = 1:mods
    mi = model_id(ii); mField = NM.modelFields{mi};
    if isfield(zon_stats.(vField), mField) && isfield(zon_stats.(vField).(mField),'mean')
        Z(:,ii) = zon_stats.(vField).(mField).mean(:);
        Zstd(:,ii) = zon_stats.(vField).(mField).std(:);
        % Div Rates
        Z2(:,ii) = zon_stats.Phytotoplankton_Specific_Division_Rate.(mField).mean(:); 
    end
end

if all(isnan(Z),'all')
    set(axz,'YLim',[ybot ytop]);
    ylabel(axz,'Latitude','Interpreter','latex');
    title(axz,'No zonal data','Interpreter','latex');
    continue
end

% --- Ensemble mean ± std across models (per latitude) ---
ens_mean  = mean(Z, 2, 'omitnan');
%ens_std   = std(Z, 0, 2, 'omitnan');
ens_std = sqrt( std(Z, 0, 2, 'omitnan').^2 + (sum(Zstd.^2,2,'omitnan')/11^2) ); % Ensemble total uncertinaty (inter and intra model std)
ens_mean2 = mean(Z2, 2, 'omitnan');

% --- X limits (tight, symmetric) ---
xmin_d = 0;%min([Z(:); ens_mean-ens_std],[],'omitnan');
xmax_d = 0.9;%max([Z(:); ens_mean+ens_std],[],'omitnan');
rngx   = max([xmin_d,xmax_d]); if ~isfinite(rngx) || rngx==0, rngx = 1; end
xlim(axz, [xmin_d,xmax_d]);
set(axz, 'YLim', [ybot ytop]);      % perfect vertical alignment with map

% ===== Background shading =====
xl = xlim(axz);  
patch('XData', [xl(1) 0 0 xl(1)], 'YData', [ybot ybot ytop ytop], ...
      'FaceColor',[0.85 0.30 0.30],'FaceAlpha',0.04,'EdgeColor','none','Parent',axz);
patch('XData', [0 xl(2) xl(2) 0], 'YData', [ybot ybot ytop ytop], ...
      'FaceColor',[0.30 0.50 0.85],'FaceAlpha',0.04,'EdgeColor','none','Parent',axz);


% ===== Individual models =====
for ii = 1:mods
    mi = model_id(ii);
    if ~all(isnan(Z(:,ii)))
        plot(axz, Z(:,ii), y_cent, 'LineWidth', 1, 'Color', cmap_mod(mi,:));
    end
end

% ===== ±1&2 σ band + mean =====
good = isfinite(ens_mean) & isfinite(ens_std) & isfinite(y_cent);
if any(good)
    % 1 STD
    X1 = [ens_mean(good)-ens_std(good); flipud(ens_mean(good)+ens_std(good))];
    Y1 = [y_cent(good);                  flipud(y_cent(good))];
    patch('XData',X1,'YData',Y1,'FaceColor',[0 0 0], ...
          'FaceAlpha',0.2,'EdgeColor','none','Parent',axz);
   % Means
    plot(axz, ens_mean, y_cent, 'k-', 'LineWidth', 6);
    plot(axz, ens_mean2, y_cent, 'b-', 'LineWidth', 2);
    
end



% ===== Axis cosmetics =====
xticks(0:.2:1)
ticklats = (-80:20:80)';
[~, ytick_map] = mfwdtran(mstruct, ticklats, zeros(size(ticklats)));
set(axz,'YTick', ytick_map, ...
        'YTickLabel', arrayfun(@(v)sprintf('%g',v), ticklats, 'uni',0), ...
        'Layer','top','FontSize',16,'TickLabelInterpreter','latex', ...
        'Box','on','YGrid','on','XGrid','on', ...
        'XAxisLocation','top');                  % (1) ticks at TOP
xline(axz, 0, '--', 'Color', 'k', 'LineWidth', 1.2);
%xlabel(axz, 'Relative Change in Phytoplankton Specific Rates (\%)', ...
%       'Interpreter','latex','FontSize',16);
%ylabel(axz, 'Latitude', 'Interpreter','latex','FontSize',20);

% =======================================================================
% (2) Small axis directly BELOW the zonal plot for global medians “ticks”
% =======================================================================
glob_med = nan(mods,1);
for ii = 1:mods
    mi    = model_id(ii);
    mFld  = NM.modelFields{mi};
    if isfield(mod_stats.(vField), mFld) && isfield(mod_stats.(vField).(mFld),'median')
        glob_med(ii) = mod_stats.(vField).(mFld).median;
    end
end

% Narrow tick axis (touching, same width)
posZ    = get(axz,'Position');
hTickAx = 0.022;         % height of tick strip (small)
gapTick = 0.00;        % tiny gap between zonal and tick strip
axTick  = axes('Units','normalized', ...
               'Position',[posZ(1), posZ(2)-gapTick-hTickAx, posZ(3), hTickAx], ...
               'Color','none','Box','off');
set(axTick,'XLim',xlim(axz),'YLim',[0 1], ...
           'XColor','none','YColor','none','HitTest','off');  % invisible frame
hold(axTick,'on');

% Draw the short colored ticks
for ii = 1:mods
    mi = model_id(ii);
    if ~isfinite(glob_med(ii)), continue; end
    plot(axTick, [glob_med(ii) glob_med(ii)], [0 0.8], ...
         '-', 'Color', cmap_mod(mi,:), 'LineWidth', 1.5);
end

% =======================================================================
% (3) Label strip below ticks: names evenly spaced, sorted by median,
%     with diagonal connectors from tick base to the label (vertical text)
% =======================================================================
valid   = isfinite(glob_med);
idxKeep = find(valid);
[~, ord] = sort(glob_med(valid), 'ascend');
idxSrt  = idxKeep(ord);                    % model indices in ascending median
nLab    = numel(idxSrt);
if nLab > 0
    % --- make the label strip a little lower & a hair taller ---
    hLabAx   = 0.034;     % label strip height (slightly taller)
    gapLab   = 0.002;     % tiny gap between tick strip and label strip
    extraDrop = 0.030;    % << additional drop to push labels farther down

    axLab  = axes('Units','normalized', ...
                  'Position',[posZ(1), ...
                              posZ(2) - gapTick - hTickAx - gapLab - hLabAx - extraDrop, ...
                              posZ(3), hLabAx], ...
                  'Color','none','Box','off');
    set(axLab,'XLim',xlim(axz),'YLim',[0 1], 'XColor','none','YColor','none', 'HitTest','off');
    hold(axLab,'on');

    % label x-positions: evenly spaced across the panel width
    xl   = xlim(axz);
    pad  = 0.03*range(xl);
    xlab = linspace(xl(1)+pad, xl(2)-pad, nLab);
%
    for j = 1:nLab
        mi   = model_id(idxSrt(j));
        name = model_names{mi};
        col  = cmap_mod(mi,:);

        % ----- vertical label text (rotated 90°) -----
        % place it low in the strip so lines are longer/clearer
        text(axLab, xlab(j)+.01, 0.01, name, ...
             'Color', col, 'FontSize', 12, 'Interpreter','latex', ...
             'HorizontalAlignment','center', 'VerticalAlignment','bottom', ...
             'Rotation', 90);

        % ----- connector from tick base (axTick, y=0) to near top of label axis -----
        [x1n, y1n] = axData2norm(axTick, glob_med(idxSrt(j)), 0.0);
        [x2n, y2n] = axData2norm(axLab,  xlab(j), 0.98);   % go to top of label strip
        annotation('line',[x1n x2n],[y1n y2n+.007], 'Color', col, 'LineWidth', 0.8);
    end
end

end % Zonal_plot (for)




%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Add Violin Distribution (biomass-weighted) + ENSEMBLE violin
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for Violin_dist = 1

    % ---- axis placement (kept from your code) ----
    dist_pos = [map_pos(1), ...
                ybot_n - max(ytop_n - ybot_n, 1e-9) - 0.05, ...
                zonal_w + map_width + gap_h, ...
                max(ytop_n - ybot_n, 1e-9)-.1];
    axd = axes('Position', dist_pos); hold(axd,'on');
    set(axd, 'TickDir','in', ...          % ticks go inward (cleaner for violins)
         'TickLength',[0.007 0.007], ...
         'LineWidth',1.0);

    % ---- axis / bins ----
    ymin = 0; ymax = 1.1;      % you hard-coded this in Code 1
    nbins = 250;

    modsN   = numel(model_id);
    centers = 1:(modsN+1);      % +1 for ENSEMBLE
    maxBarW = 0.50;

    y_edges   = linspace(ymin, ymax, nbins+1);
    y_centers = 0.5*(y_edges(1:end-1) + y_edges(2:end));
    y = y_centers(:);

    % ---- background halves ----
    xlim(axd, [.25, modsN+1.75]);
    ylim(axd, [ymin, ymax]);
    
    set(axd,'Layer','top','FontSize',15,'TickLabelInterpreter','latex');
    set(axd,'XTick',centers,'XTickLabel',[]);

    % ---- y label (you had a couple versions; keep one) ----
    ylabel(axd, 'Grazing Pressure, $l_g$ (d$^{-1}$)', 'Interpreter','latex','fontsize',20);

    % ---- x labels (models + ensemble) ----
    for ii = 1:modsN
        text(axd, centers(ii), ymin - 0.04*(ymax-ymin), model_names{ii}, ...
            'HorizontalAlignment','center', 'VerticalAlignment','top', ...
            'FontSize',16, 'Color', cmap_mod(ii,:), 'Interpreter','latex');
    end
    text(axd, centers(modsN+1), ymin - 0.04*(ymax-ymin), '\textbf{Ensemble}', ...
        'HorizontalAlignment','center', 'VerticalAlignment','top', ...
        'FontSize',16, 'Color','k', 'Interpreter','latex');

    xline(axd, centers, ':', 'Color', [0 0 0 0.25]);

    % -------- per-model half-violins (your original loop) --------
    cd /Users/twrohr/Desktop/data.nosync/SSP585/Abridged/regridded2
    mods = dir('*.mat');

    % ACCUMULATORS FOR ENSEMBLE violin (pooled weighted samples)
    g_all  = [];  w_all  = [];
    ng_all = [];  wn_all = [];

    for ii = 1:modsN

        % pull values + biomass weights
        S = load(mods(mod_order(model_id(ii))).name, 'CoreVar_Inv_adj', 'Bio_Weights');
        V = S.CoreVar_Inv_adj;
        A = S.Bio_Weights;

        % Select variables to plot on left (g) and Right (ng) side of violin plots
        g   = squeeze(mean(V(11,1:60,:,:),2,'omitnan'));   % 11 for Fig4d; 16 For Sup Fig 5.
        ng  = squeeze(mean(V(17,1:60,:,:),2,'omitnan'));   % 17 for Fig4d;  6 For Sup Fig 5.

        % valid masks
        maskG   = isfinite(g)  & isfinite(A) & (A > 0);
        maskNG  = isfinite(ng) & isfinite(A) & (A > 0);

        % “real data” masks: ignore pure zeros (matches Code 2)
        maskG_data  = maskG  & (g  ~= 0);
        maskNG_data = maskNG & (ng ~= 0);

        % accumulate ensemble samples
        if any(maskG_data(:))
            g_all = [g_all; g(maskG_data)];   %#ok<AGROW>
            w_all = [w_all; A(maskG_data)];   %#ok<AGROW>
        end
        if any(maskNG_data(:))
            ng_all = [ng_all; ng(maskNG_data)]; %#ok<AGROW>
            wn_all = [wn_all; A(maskNG_data)];  %#ok<AGROW>
        end

        % nothing to draw?
        if ~any(maskG_data(:)) && ~any(maskNG_data(:)), continue; end

        % KDE grid along y
        fG   = zeros(size(y));
        fNG  = zeros(size(y));
        if any(maskG_data(:))
            fG = ksdensity(g(maskG_data),  y, 'Weights', A(maskG_data), ...
                           'Support','unbounded', 'Function','pdf');
        end
        if any(maskNG_data(:))
            fNG = ksdensity(ng(maskNG_data), y, 'Weights', A(maskNG_data), ...
                            'Support','unbounded', 'Function','pdf');
        end

        % independent width scaling
        col = cmap_mod(model_id(ii),:);
        mG   = max(fG);  if ~(isfinite(mG) && mG>0),  mG  = 1; end
        mNG  = max(fNG); if ~(isfinite(mNG) && mNG>0), mNG = 1; end
        wG   = -(fG  / mG)  * maxBarW;   % left
        wNG  =  (fNG / mNG) * maxBarW;   % right

        % LEFT (g): filled
        if any(maskG_data(:)) && any(isfinite(wG))
            xL = centers(ii) + wG;
            xC = centers(ii) + zeros(size(y));
            patch('Parent',axd, ...
                  'XData',[xL; flipud(xC)], 'YData',[y;  flipud(y)], ...
                  'FaceColor',col, 'FaceAlpha',0.7, 'EdgeColor','none');
        end

        % RIGHT (ng): outline-ish
        if any(maskNG_data(:)) && any(isfinite(wNG))
            xC = centers(ii) + zeros(size(y));
            xR = centers(ii) + wNG;
            patch('Parent',axd, ...
                  'XData',[xC; flipud(xR)], 'YData',[y;  flipud(y)], ...
                  'FaceColor',col, 'FaceAlpha',0.1, ...
                  'EdgeColor',cmap_mod(ii,:),'LineWidth',0.5);
        end

        % medians & IQR (your original, but use *_data masks)
        g_med = NaN; g_25 = NaN; g_75 = NaN;
        ng_med = NaN; ng_25 = NaN; ng_75 = NaN;
        if any(maskG_data(:))
            [g_med, g_25, g_75] = local_wquant3(g(maskG_data), A(maskG_data));
        end
        if any(maskNG_data(:))
            [ng_med, ng_25, ng_75] = local_wquant3(ng(maskNG_data), A(maskNG_data));
        end

        cap  = 0.08; lw = 1.2;

        if isfinite(g_med)
            wGm = interp1(y, wG, g_med, 'linear','extrap');
            if isfinite(wGm)
                xIQR = centers(ii) + 0.5*wGm;  colk = 'k';
                h1 = plot(axd, [xIQR xIQR], [g_25 g_75],   'Color',colk,'LineWidth', lw);
                h2 = plot(axd, [xIQR-cap/2, xIQR+cap/2], [g_25 g_25], 'Color',colk,'LineWidth', 0.7);
                h3 = plot(axd, [xIQR-cap/2, xIQR+cap/2], [g_75 g_75], 'Color',colk,'LineWidth', 0.7);
                h4 = scatter(axd, xIQR, g_med, 50, 'w', 'filled', 'Marker','d', ...
                             'MarkerEdgeColor','k', 'LineWidth',1);
                uistack([h1 h2 h3 h4],'top');
            end
        end

        if isfinite(ng_med)
            wNGm = interp1(y, wNG, ng_med, 'linear','extrap');
            if isfinite(wNGm)
                xIQR = centers(ii) + 0.5*wNGm;  colk = 'k';
                h1 = plot(axd, [xIQR xIQR], [ng_25 ng_75], 'Color',colk,'LineWidth', lw);
                h2 = plot(axd, [xIQR-cap/2, xIQR+cap/2], [ng_25 ng_25], 'Color',colk,'LineWidth', 0.7);
                h3 = plot(axd, [xIQR-cap/2, xIQR+cap/2], [ng_75 ng_75], 'Color',colk,'LineWidth', 0.7);
                h4 = scatter(axd, xIQR, ng_med, 50, 'k', 'filled', 'Marker','d', ...
                             'MarkerEdgeColor','k', 'LineWidth',1);
                uistack([h1 h2 h3 h4],'top');
            end
        end

    end % per-model loop

    % ------------------------------------------------------------------
    % ENSEMBLE violin at far right (centers(modsN+1)), like Code 2
    % ------------------------------------------------------------------
    cx          = centers(modsN+1);
    maxBarW_ens = maxBarW;

    % Left (g_all): solid black
    if ~isempty(g_all) && any(w_all>0)
        fg = ksdensity(g_all, y, 'Weights', w_all, ...
                       'Support','unbounded','Function','pdf');
        if any(fg>0), fg = fg ./ max(fg); else, fg = zeros(size(fg)); end
        wG = -maxBarW_ens * fg;

        patch(axd, [cx + wG; cx*ones(size(y))], [y; flipud(y)], ...
              'k', 'EdgeColor','none', 'FaceAlpha', 0.70);

        [g_medE, g_q25E, g_q75E] = local_wquant3(g_all, w_all);
        cap = 0.08; lw = 1.1; xIQR = cx - 0.5*maxBarW_ens;
        plot(axd, [xIQR xIQR],[g_q25E min(ymax,g_q75E)],'k-','LineWidth',lw);
        plot(axd, [xIQR-cap/2 xIQR+cap/2],[g_q25E g_q25E],'k-','LineWidth',lw);
        if ymax > g_q75E
            plot(axd, [xIQR-cap/2 xIQR+cap/2],[g_q75E g_q75E],'k-','LineWidth',lw);
        end
        scatter(axd, xIQR, g_medE, 55, 'w', 'filled', 'Marker','d', ...
                'MarkerEdgeColor','k', 'LineWidth',1.1);
    end

    % Right (ng_all): translucent black (MISSING IN YOUR CURRENT FIG)
    if ~isempty(ng_all) && any(wn_all>0)
        fn = ksdensity(ng_all, y, 'Weights', wn_all, ...
                       'Support','unbounded','Function','pdf');
        if any(fn>0), fn = fn ./ max(fn); else, fn = zeros(size(fn)); end
        wN =  maxBarW_ens * fn;
    
        % half-violin patch (right side)
        patch(axd, [cx + wN; cx*ones(size(y))], [y; flipud(y)], ...
              'k', 'EdgeColor','k', 'FaceAlpha', 0.10, 'LineWidth', 0.6);
    
        % stats markers (right side)
        [ng_medE, ng_q25E, ng_q75E] = local_wquant3(ng_all, wn_all);
        cap = 0.08; lw = 1.1; xIQR = cx + 0.5*maxBarW_ens;
    
        plot(axd, [xIQR xIQR], [ng_q25E min(ymax,ng_q75E)], 'k-', 'LineWidth', lw);
        plot(axd, [xIQR-cap/2 xIQR+cap/2], [ng_q25E ng_q25E], 'k-', 'LineWidth', lw);
        if ymax > ng_q75E
            plot(axd, [xIQR-cap/2 xIQR+cap/2], [ng_q75E ng_q75E], 'k-', 'LineWidth', lw);
        end
        scatter(axd, xIQR, ng_medE, 55, 'k', 'filled', ...
                'Marker','d', 'MarkerEdgeColor','k', 'LineWidth',1.1);
    end

        % Make a full box and inward ticks
        set(axd, 'Box','on', ...
                 'TickDir','in', ...
                 'TickLength',[0.010 0.010], ...
                 'LineWidth',1.1);
        
        % Choose y-ticks explicitly (nice, avoids weird auto ticks)
        yt = 0:0.2:1.0;   % <-- adjust if you want (must sit within [ymin ymax])
        set(axd,'YTick',yt);
        
        % LEFT axis (grazing)
        yyaxis(axd,'left');
        set(axd,'YColor','k');   % make sure left spine stays black
        ylabel(axd,'Grazing Pressure, $l_g$ (d$^{-1}$)','Interpreter','latex','fontsize',20);  % Figure 4b,c
        %ylabel(axd,'Loss Rates, $l$ (d$^{-1}$)','Interpreter','latex','fontsize',20);           % Supplementary Figure 5a

        % RIGHT axis (non-grazing): same ticks/limits, ticks inward
        yyaxis(axd,'right');
        set(axd,'YColor','k', 'YTick',yt, 'YLim',[ymin ymax]);   % match left exactly
        ylabel(axd,'Non-grazing Pressure, $l_{ng}$ (d$^{-1}$)','Interpreter','latex','fontsize',20); % Figure 4b,c
        %ylabel(axd,'Division Rate, $\mu$ (d$^{-1}$)','Interpreter','latex','fontsize',20);           % Supplementary Figure 5a
        
        % Go back to left so any later plotting uses left axis
        yyaxis(axd,'left');

end % Violin_dist

    




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%  Helper Functions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xn, yn] = axData2norm(ax, x, y)
% Convert data coords (x,y) in axes AX to figure normalized coords
    pos = get(ax,'Position');     % [x y w h] in normalized figure units
    xl  = get(ax,'XLim');  yl = get(ax,'YLim');
    xn  = pos(1) + (x - xl(1)) / (xl(2)-xl(1)) * pos(3);
    yn  = pos(2) + (y - yl(1)) / (yl(2)-yl(1)) * pos(4);
end

function counts = histcounts_weighted(x, edges, w)
    % x: values (vector), edges: bin edges, w: weights (same size as x)
    % returns 1 x (numBins) row vector of weighted sums per bin
    idx   = discretize(x, edges);                 % bin index per sample
    valid = ~isnan(idx) & isfinite(w);
    counts_col = accumarray(idx(valid), w(valid), [numel(edges)-1, 1], @sum, 0);
    counts = counts_col.';                        % make it a row vector
end

function q = local_wquant(x,w,p)
    % returns NaN if not enough information
    ok = isfinite(x) & isfinite(w) & (w > 0);
    x = x(ok); w = w(ok);
    if numel(x) == 0 || sum(w) <= 0
        q = NaN; return
    end
    % sort by x
    [xs, ix] = sort(x(:));
    ws = w(:); ws = ws(ix);

    % cumulative weights in [0,1]
    cw = cumsum(ws);
    cw = cw ./ cw(end);

    % collapse duplicate cumulative weights before interp1
    [cu, iu] = unique(cw, 'stable');
    xsu = xs(iu);

    % if still only one point, return it
    if numel(cu) == 1
        q = xsu(1); return
    end

    % interpolate the inverse CDF
    q = interp1(cu, xsu, p, 'linear', 'extrap');
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
    % IQR
    q25 = interp1(cu, xsu, 0.25, 'linear','extrap');
    q50 = interp1(cu, xsu, 0.50, 'linear','extrap');
    q75 = interp1(cu, xsu, 0.75, 'linear','extrap');
    % pm 10%
    % q25 = interp1(cu, xsu, 0.40, 'linear','extrap');
    % q50 = interp1(cu, xsu, 0.50, 'linear','extrap');
    % q75 = interp1(cu, xsu, 0.60, 'linear','extrap');
end
