%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Supplemntary Figure 10 - Sensitvity and Driver Violions %%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plots ensemble distribution of raw change, percent driver change, and sensitivty
% (normalizing by percent change in driver)
% Created Jan 1 2026
%
% ADJUSTMENTS (this edit):
%  1) Add ENSEMBLE box/whisker inside the same loop as models, as "ii = nmod+1"
%     and place it at jitter position = 0 so spacing is even and centered.
%     Also make ensemble marker size + line width bigger.
%  2) var=47: mirror ensemble violin left/right again.
%     But for individual model box/whiskers: DO NOT repeat left+right;
%     instead compute ONE set of per-model quantiles from var47 (biomass-weighted),
%     and plot them centered across both sides at cx + jitter, with jitterWidth doubled.
%  3) Replace 'g'/'ng' with 'grazing'/'non-grazing' and move just above axis
%     (outside plot area). No labels for var=47.
%
% Notes:
%  - For i != 47: left = i-10 (grazing), right = i-5 (non-grazing)
%  - For i == 47: left = 47, right = 47 (mirrored ensemble)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear; close all; clc;

cd /Users/twrohr/Desktop/Matlab/Matlab/Grazing/CMIP6_Comp/NCC_Figs
addpath(genpath('/Users/twrohr/Desktop/Matlab/Matlab/Utilities'))
dataDir = '/Users/twrohr/Desktop/data.nosync/SSP585/Abridged/regridded2';

for startup = 1 

    % ---------------- Model names / ordering / 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];
    model_id    = 1:11;
    nmod        = numel(model_id);

    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];
    cmap_mod = cmap_mod(1:11,:);

    wedge_colors = [
        0.90 0.55 0.25;   % SST
        0.35 0.60 0.85;   % MLD
        0.25 0.75 0.65;   % Diatoms
        0.60 0.75 0.45;   % Phyto
        0.85 0.55 0.75;   % Zoo
        0.20 0.20 0.20];  % Interactions
    lighter_colors = wedge_colors + 0.6 .* (1 - wedge_colors);

    % Reorder to match your plotting order (unchanged from your snippet)
    wedge_colors   = wedge_colors([4,5,1,3,2,6],:);
    lighter_colors = lighter_colors([4,5,1,3,2,6],:);

    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 Length Scalee',             '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_symbols = {'$\Delta l_{g,\Sigma P}$', '$\Delta l_{ng,\Sigma P}$', ...
                    '$\Delta l_{g,\Sigma Z}$','$\Delta l_{ng,\Sigma Z}$', ...
                    '$\Delta l_{g,\ T}$',  '$\Delta l_{ng,T}$', ...
                    '$\Delta l_{g,CC}$', '$\Delta l_{ng, CC}$', ...
                    '$\Delta l_{g, VD}$','$\Delta l_{ng, VD}$', ...
                    '$\Delta l_{g, Int}$'};

    decomp_symbols = {'$\Delta l_{T}$','$\Delta l_{VD}$',...
                     '$\Delta l_{CC}$','$\Delta l_{\Sigma P}$',...
                     '$\Delta l_{\Sigma Z}$','$\Delta l_{Int}$'};
    decomp_symbols = decomp_symbols([4,5,1,3,2,6]);

    decomp_symbols2 = {'$l_{g,\Sigma P}$', '$l_{ng,\Sigma P}$', ...
                    '$l_{g,\Sigma Z}$','$l_{ng,\Sigma Z}$', ...
                    '$l_{g,T}$',  '$l_{ng,T}$', ...
                    '$l_{g,CC}$', '$l_{ng,CC}$', ...
                    '$l_{g,VD}$','$l_{mg,VD}$', ...
                    '$l_{Int}$'};

    Sens_names2 =  { '$\frac{\Delta \; l_{P}}{\Delta\;\Sigma \textrm{P}}$' , '$\frac{\Delta \; l_{\Sigma \textrm{Z}} }{\Delta \; \Sigma Z}$' ,...
                    '$\frac{\Delta \; l_{T} }{\Delta\;  \textrm{SST}}$' , '$\frac{\Delta \; l_{\textrm{CC}}}{\Delta \; f_{\mathrm{SP}}}$' ,...
                     '$\frac{\Delta \;  l_{\mathrm{VD}}}{\Delta\;  \mathrm{D}_\mathrm{L} }$' };
    Sens_names2  =  { '$\frac{\Delta \; l_{T}}{\Delta_{\%} \;  \mathrm{SST}}$' , '$\frac{\Delta \; l_{VD} }{\Delta_{\%}  \; \mathrm{D}_\mathrm{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{Z} }$' };

    Sens_names2 = Sens_names2([4,5,1,3,2]);

end


% ---------------- Panel variable lists ----------------
mvar_list   = [40,41,37,39,38,47];
mirror_id   = [false,false,false,false,false,true];

mvar_list2  = [1,2,14,12,42];
mirror_id2  = [true,true,true,true,true];

mvar_list3  = [4,5,1,3,2]+10;
mirror_id3  = [false,false,false,false,false];


% ---------------- Plot settings (shared) ----------------
nbins = 250; maxBarW = 0.5; dx = 1.0; x0 = 1.0; halfOffset = 0.25; jitterWidth = 0.2;
markerSize = 22; lw_iqr = 1.5; capSize = 6;
markerSizeE = 34; lw_iqrE = 2.1; capSizeE = 8;
fsPct = 16; showEnsembleStats = true; doTDmaskRight_default = false;

% y-lims (keep your current)
ymin1 = -0.06; ymax1 = 0.08; 
ymin2 = -50;   ymax2 = 50;
ymin3 = -.0042;   ymax3 = .0055;

% y grids
y1 = linspace(ymin1,ymax1,nbins+1); y1 = 0.5*(y1(1:end-1)+y1(2:end)); y1 = y1(:);
y2 = linspace(ymin2,ymax2,nbins+1); y2 = 0.5*(y2(1:end-1)+y2(2:end)); y2 = y2(:);
y3 = linspace(ymin3,ymax3,nbins+1); y3 = 0.5*(y3(1:end-1)+y3(2:end)); y3 = y3(:);

% ---------------- Files ----------------
cd(dataDir); mods = dir('*.mat'); if isempty(mods), error('No .mat files found in %s', dataDir); end

% ---------------- CoreVar_Names must exist ----------------
if ~exist('CoreVar_Names','var') || isempty(CoreVar_Names)
    error('CoreVar_Names must be defined in the workspace (used for x-label annotations).');
end

% ---------------- Figure with 2 panels ----------------
figure('Color','w','Units','normalized','Position',[0.05 0.06 0.92 0.88]);
W1 = 0.75;          % ax1 width
W2 = .9;      % scale by number of variables
L1 = 0.06;

ax1 = axes('Position',[L1 0.69 W1 0.26]); hold(ax1,'on');
ax2 = axes('Position',[L1 0.36 W1 0.26]); hold(ax2,'on');
ax3 = axes('Position',[L1 0.03 W2 0.26]); hold(ax3,'on');



plot_panel(ax1, mods, dataDir, mvar_list3, mirror_id3, 'Sens', y3, ymin3, ymax3, ...
    model_id, nmod, mod_order, cmap_mod, wedge_colors, lighter_colors, dx, x0, maxBarW, halfOffset, jitterWidth, ...
    markerSize, lw_iqr, capSize, markerSizeE, lw_iqrE, capSizeE, fsPct, showEnsembleStats, doTDmaskRight_default, Sens_names2);

        % ----  xtick labels
        set(gca, 'FontSize', 18)
        xticks(ax1, 1:numel(Sens_names2))
        xticklabels(ax1, [])
        fig   = ancestor(ax1,'figure');
        axpos = ax1.Position;
        xt    = ax1.XTick;
        xl    = ax1.XLim;
        
        x_norm = axpos(1) + axpos(3) * (xt - xl(1)) ./ diff(xl);
        y_norm = axpos(2) - 0.045;   % vertical offset (tune if needed)
        
        for i = 1:numel(Sens_names2)
            annotation(fig,'textbox', ...
                [x_norm(i)-0.02, y_norm, 0.04, 0.04], ...
                'String', Sens_names2{i}, ...
                'Interpreter','latex', ...
                'FontSize', 25, ...
                'HorizontalAlignment','center', ...
                'VerticalAlignment','top', ...
                'EdgeColor','none');
        end 
        ylabel(ax1, {'Sensitivity of Loss Rates to ';'Relative Change in Driver';'($\textrm{d}^{-1} \; \% ^{-1} )$'},'FontSize',20, 'Interpreter','latex')
   
plot_panel(ax2, mods, dataDir, mvar_list2, mirror_id2, 'CoreVar_Inv_pc', y2, ymin2, ymax2, ...
    model_id, nmod, mod_order, cmap_mod, lighter_colors, lighter_colors, dx, x0, maxBarW, halfOffset, jitterWidth, ...
    markerSize, lw_iqr, capSize, markerSizeE, lw_iqrE, capSizeE, fsPct, showEnsembleStats, doTDmaskRight_default, CoreVar_Names);

     xticklabels(ax2, {'$\Delta_{\%} \Sigma \textrm{P}$','$\Delta_{\%} \Sigma \textrm{Z}$','$\Delta_{\%} \textrm{SST}$','$\Delta_{\%}$ f$_{\textrm{sp}}$', '$\Delta_{\%}$ \textrm{D}$_{\textrm{L}}$'})
     set(gca, 'FontSize', 18)
     ylabel(ax2,{'Relative Change in Driver'; '($\%$)'},'FontSize',20, 'Interpreter','latex')
   

plot_panel(ax3, mods, dataDir, mvar_list, mirror_id, 'CoreVar_Inv_rw', y1, ymin1, ymax1, ...
    model_id, nmod, mod_order, cmap_mod, wedge_colors, lighter_colors, dx, x0, maxBarW, halfOffset, jitterWidth, ...
    markerSize, lw_iqr, capSize, markerSizeE, lw_iqrE, capSizeE, fsPct, showEnsembleStats, doTDmaskRight_default, decomp_symbols2);

    xticklabels(ax3, decomp_symbols); set(gca, 'FontSize', 18)
    ylabel(ax3, {'First-order effect on'; 'loss rates; $\Delta l_i$ (d$^{-1} $)'},'FontSize',20, 'Interpreter','latex')
   
    

% ---------------- Legend -----------------
for write_legend =1

    % ===== Legend in the blank right-hand space next to ax2+ax3 (spanning both) =====
pos2 = get(ax1,'Position'); pos3 = get(ax2,'Position');  % [x y w h] normalized
gap  = 0.010;  padR = 0.012;                              % tweak if needed
xLeft = pos2(1) + pos2(3) + gap;                          % to the right of panels
yBot  = pos3(2);
yTop  = pos2(2) + pos2(4);
wLeg  = max(0.001, 1 - xLeft - padR);                     % whatever space remains
hLeg  = yTop - yBot;

axLeg = axes('Position',[xLeft yBot wLeg hLeg],'Color','none','XLim',[0 1],'YLim',[0 1], ...
             'XTick',[],'YTick',[],'Box','off','XColor','none','YColor','none'); axis(axLeg,'off');

ys = linspace(0.95,0.05,nmod);  % even vertical spread across BOTH panels
xText = 0.02; txtFS = 16; padAfter = 0.03; targetX = 0.78; targetY = 0.52; alpha = 0.92;

for j=1:nmod
    c = cmap_mod(j,:);
    ht = text(axLeg,xText,ys(j),model_names{j},'Color',c,'Interpreter','latex', ...
              'HorizontalAlignment','left','VerticalAlignment','middle','FontSize',txtFS,'Clipping','off');
    drawnow;
    ext = get(ht,'Extent'); x1 = ext(1)+ext(3)+padAfter; y1 = ext(2)+0.5*ext(4);
    x2 = x1 + alpha*(targetX-x1); y2 = y1 + alpha*(targetY-y1);
    line(axLeg,[x1 x2],[y1 y2],'Color',c,'LineWidth',1.3,'Clipping','off');
end


end


% ---------------- Helpers ----------------
function plot_panel(ax, mods, dataDir, mvar_list, mirror_id, Vfield, y, ymin, ymax, model_id, nmod, mod_order, cmap_mod, L_colors, R_colors, dx, x0, maxBarW, halfOffset, jitterWidth, markerSize, lw_iqr, capSize, markerSizeE, lw_iqrE, capSizeE, fsPct, showEnsembleStats, doTDmaskRight_default, Var_Names)

keep = ~isnan(mvar_list);
mvar_list = mvar_list(keep);
mirror_id = mirror_id(keep);
nVar = numel(mvar_list);
xCtr = x0 + (0:nVar-1)*dx;

xlim(ax,[xCtr(1)-0.6*dx, xCtr(end)+0.6*dx]); ylim(ax,[ymin ymax]);
patch(ax,[xlim(ax) fliplr(xlim(ax))],[0 0 ymax ymax],[0.3 0.5 0.85],'EdgeColor','none','FaceAlpha',0.04);
patch(ax,[xlim(ax) fliplr(xlim(ax))],[ymin ymin 0 0],[0.85 0.3 0.3],'EdgeColor','none','FaceAlpha',0.04);
yline(ax,0,'k-','LineWidth',0.9);
for k=1:nVar, xline(ax,xCtr(k),':','Color',[0 0 0 0.20],'LineWidth',0.9); end
set(ax,'XTick',xCtr,'XTickLabel',[],'TickLabelInterpreter','latex','FontSize',15);
ylabel(ax,'Rate contribution (d$^{-1}$)','Interpreter','latex','FontSize',18);
grid(ax,'on'); ax.GridAlpha=0.12; set(ax,'Box','on','Layer','top');

slotInts = -6:5; slotVals = (slotInts/6)*jitterWidth; i0 = find(slotInts==0,1); modelSlotVals = slotVals; modelSlotVals(i0) = [];

for k=1:nVar
    mvar_id = mvar_list(k); 
    cx = xCtr(k); isMirrored = logical(mirror_id(k));
    if isMirrored, left_id=mvar_id; right_id=mvar_id;      doTDmaskRight=false; aL = .5; aR = .5;   
    else,          left_id=mvar_id-10; right_id=mvar_id-5; doTDmaskRight=false; aL = .8; aR = .33;  end

    L_all=[]; wL_all=[]; R_all=[]; wR_all=[]; L_med=nan(nmod,1); L_q25=nan(nmod,1); L_q75=nan(nmod,1); R_med=nan(nmod,1); R_q25=nan(nmod,1); R_q75=nan(nmod,1);
    C_med=nan(nmod,1); C_q25=nan(nmod,1); C_q75=nan(nmod,1);
    
    % loop models
    for ii=1:nmod
        S = load(fullfile(dataDir,mods(mod_order(model_id(ii))).name), Vfield,'Bio_Weights');
        V = S.(Vfield); A = S.Bio_Weights;

        L = squeeze(V(left_id,:,:)); maskL = isfinite(L)&isfinite(A)&(A>0)&(L~=0);
        if any(maskL(:))
            wL = A(maskL); sL = sum(wL,'omitnan'); if ~(isfinite(sL)&&sL>0), sL=1; end; wL_eq = wL./sL;
            L_all=[L_all;L(maskL)]; wL_all=[wL_all;wL_eq];
            if isMirrored, [C_med(ii),C_q25(ii),C_q75(ii)] = local_wquant3(L(maskL),wL); else, [L_med(ii),L_q25(ii),L_q75(ii)] = local_wquant3(L(maskL),wL); end
        end

        R = squeeze(V(right_id,:,:));
        if doTDmaskRight
            BU=squeeze(V(1,:,:)); TD=squeeze(V(6,:,:)); TDmask = sign(BU).*sign(TD)==-1;
            maskR = isfinite(R)&isfinite(A)&(A>0)&TDmask&(R~=0);
        else
            maskR = isfinite(R)&isfinite(A)&(A>0)&(R~=0);
        end
        if any(maskR(:))
            wR = A(maskR); sR = sum(wR,'omitnan'); if ~(isfinite(sR)&&sR>0), sR=1; end; wR_eq = wR./sR;
            R_all=[R_all;R(maskR)]; wR_all=[wR_all;wR_eq];
            if ~isMirrored, [R_med(ii),R_q25(ii),R_q75(ii)] = local_wquant3(R(maskR),wR); end
        end
    end

    % Add Violins
    colL = L_colors(k,:); colR = R_colors(k,:);

    if ~isempty(L_all) && any(wL_all>0)
        fL = ksdensity(L_all,y,'Weights',wL_all,'Support','unbounded','Function','pdf'); if any(fL>0), fL=fL./max(fL); else, fL=zeros(size(fL)); end; wL = -maxBarW*fL;
        patch(ax,[cx+wL; cx*ones(size(y))],[y; flipud(y)],colL,'EdgeColor','none','FaceAlpha',aL); plot(ax,cx+wL,y,'Color',[colL 0.35],'LineWidth',0.50);
        f_pos_L = weighted_frac(L_all,wL_all,L_all>0); f_neg_L = weighted_frac(L_all,wL_all,L_all<0); dy=0.005*(ymax-ymin); xTxt=cx-0.17*maxBarW;
        if isfinite(f_pos_L)&&round(100*f_pos_L)>0, text(ax,xTxt,ymax-dy,sprintf('%d\\%%',round(100*f_pos_L)),'HorizontalAlignment','center','VerticalAlignment','top','Interpreter','latex','Color','k','FontSize',fsPct,'Clipping','on'); end
        if isfinite(f_neg_L)&&round(100*f_neg_L)>0, text(ax,xTxt,ymin+dy,sprintf('%d\\%%',round(100*f_neg_L)),'HorizontalAlignment','center','VerticalAlignment','bottom','Interpreter','latex','Color','k','FontSize',fsPct,'Clipping','on'); end
    else
        wL = zeros(size(y));
    end

    if ~isempty(R_all) && any(wR_all>0)
        fR = ksdensity(R_all,y,'Weights',wR_all,'Support','unbounded','Function','pdf'); if any(fR>0), fR=fR./max(fR); else, fR=zeros(size(fR)); end; wR = maxBarW*fR;
        patch(ax,[cx+wR; cx*ones(size(y))],[y; flipud(y)],colR,'EdgeColor',colR,'FaceAlpha',aR,'LineWidth',0.9); plot(ax,cx+wR,y,'Color',[colR 0.45],'LineWidth',0.60);
        f_pos_R = weighted_frac(R_all,wR_all,R_all>0); f_neg_R = weighted_frac(R_all,wR_all,R_all<0); dy=0.005*(ymax-ymin); xTxt=cx+0.17*maxBarW;
        if isfinite(f_pos_R)&&round(100*f_pos_R)>0, text(ax,xTxt,ymax-dy,sprintf('%d\\%%',round(100*f_pos_R)),'HorizontalAlignment','center','VerticalAlignment','top','Interpreter','latex','Color','k','FontSize',fsPct,'Clipping','on'); end
        if isfinite(f_neg_R)&&round(100*f_neg_R)>0, text(ax,xTxt,ymin+dy,sprintf('%d\\%%',round(100*f_neg_R)),'HorizontalAlignment','center','VerticalAlignment','bottom','Interpreter','latex','Color','k','FontSize',fsPct,'Clipping','on'); end
    else
        wR = zeros(size(y));
    end

    if showEnsembleStats
        if ~isempty(L_all) && any(wL_all>0), [LmE,L25E,L75E] = local_wquant3(L_all,wL_all); else, LmE=NaN; L25E=NaN; L75E=NaN; end
        if ~isempty(R_all) && any(wR_all>0), [RmE,R25E,R75E] = local_wquant3(R_all,wR_all); else, RmE=NaN; R25E=NaN; R75E=NaN; end
        if isMirrored, CmE=LmE; C25E=L25E; C75E=L75E; end
    end

    slotVals47 = (slotInts/6)*(2*jitterWidth); modelSlotVals47 = slotVals47; modelSlotVals47(i0) = [];
    for ii=1:(nmod+1)
        isEns = (ii==nmod+1);
        if isEns, c=[0 0 0]; ms=markerSizeE; lw=lw_iqrE; cs=capSizeE; jitThis=0; else, c=cmap_mod(ii,:); ms=markerSize; lw=lw_iqr; cs=capSize; jitThis = (isMirrored)*modelSlotVals47(ii) + (~isMirrored)*modelSlotVals(ii); end

        if ~isMirrored
            if isEns, med=LmE; q25=L25E; q75=L75E; else, med=L_med(ii); q25=L_q25(ii); q75=L_q75(ii); end
            if isfinite(med)
                xL = cx-halfOffset+jitThis; errorbar(ax,xL,med,med-q25,q75-med,'LineStyle','none','Color',c,'LineWidth',lw,'CapSize',cs);
                if isEns, sc=scatter(ax,xL,med,ms*1.35,'w','filled','Marker','d','MarkerEdgeColor','k','LineWidth',1.2); else, sc=scatter(ax,xL,med,ms,c,'filled','MarkerEdgeColor','k','LineWidth',0.40); end
                uistack(sc,'top');
            end
            if isEns, med=RmE; q25=R25E; q75=R75E; else, med=R_med(ii); q25=R_q25(ii); q75=R_q75(ii); end
            if isfinite(med)
                xR = cx+halfOffset+jitThis; errorbar(ax,xR,med,med-q25,q75-med,'LineStyle','none','Color',c,'LineWidth',lw,'CapSize',cs);
                if isEns, sc=scatter(ax,xR,med,ms*1.35,'w','filled','Marker','d','MarkerEdgeColor','k','LineWidth',1.2); else, sc=scatter(ax,xR,med,ms,c,'filled','MarkerEdgeColor','k','LineWidth',0.40); end
                uistack(sc,'top');
            end
        else
            if isEns, med=CmE; q25=C25E; q75=C75E; else, med=C_med(ii); q25=C_q25(ii); q75=C_q75(ii); end
            if isfinite(med)
                xC = cx+jitThis; errorbar(ax,xC,med,med-q25,q75-med,'LineStyle','none','Color',c,'LineWidth',lw,'CapSize',cs);
                if isEns, sc=scatter(ax,xC,med,ms*1.35,'w','filled','Marker','d','MarkerEdgeColor','k','LineWidth',1.2); else, sc=scatter(ax,xC,med,ms,c,'filled','MarkerEdgeColor','k','LineWidth',0.30); end
                uistack(sc,'top');
            end
        end
    end

    if ~isMirrored
        axPos=ax.Position; xl=xlim(ax); dx_edge=0.05; xL_edge=cx-dx_edge; xR_edge=cx+dx_edge;
        xL_norm=axPos(1)+(xL_edge-xl(1))/(xl(2)-xl(1))*axPos(3); xR_norm=axPos(1)+(xR_edge-xl(1))/(xl(2)-xl(1))*axPos(3); y_norm=axPos(2)+axPos(4)+0.003;
        annotation('textbox',[xL_norm-0.06 y_norm 0.06 0.03],'String','\bf grazing','Interpreter','latex','HorizontalAlignment','right','VerticalAlignment','bottom','EdgeColor','none','Color',colL,'FontSize',14);
        annotation('textbox',[xR_norm y_norm 0.08 0.03],'String','\bf non-grazing','Interpreter','latex','HorizontalAlignment','left','VerticalAlignment','bottom','EdgeColor','none','Color',colR,'FontSize',14);
    end

end

% % ---- “x labels” from CoreVar_Names in same band logic as your prior annotations ----
% axPos=ax.Position; xl=xlim(ax); x2n=@(x) axPos(1)+(x-xl(1))/(xl(2)-xl(1))*axPos(3); y_norm=axPos(2)-0.05;
% for k=1:nVar
%     mvar_id=mvar_list(k); if isnan(mvar_id), continue; end
%     cx=xCtr(k);
%     if mirror_id(k)
%         lab = Var_Names{mvar_id};
%         annotation('textbox',[x2n(cx)-0.10 y_norm 0.20 0.04],'String',lab,'Interpreter','latex','HorizontalAlignment','center','VerticalAlignment','bottom','EdgeColor','none','FontSize',18);
%     else
%         lid = mvar_id-10; rid = mvar_id-5; labL = Var_Names{lid}; labR = Var_Names{rid};
%         annotation('textbox',[x2n(cx-0.10)-0.12 y_norm 0.12 0.04],'String',labL,'Interpreter','latex','HorizontalAlignment','right','VerticalAlignment','bottom','EdgeColor','none','FontSize',18);
%         annotation('textbox',[x2n(cx+0.10) y_norm 0.14 0.04],'String',labR,'Interpreter','latex','HorizontalAlignment','left','VerticalAlignment','bottom','EdgeColor','none','FontSize',18);
%     end
% end

end

function [q50,q25,q75] = local_wquant3(x,w)
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, q25=xsu; q50=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 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
