function scatterCompare4Models( ...
    T_quad_rigid, T_quad_flex, T_simple_rigid, T_simple_flex, ...
    xVar, yVar, fixed_qr, fixed_qf, fixed_sr, fixed_sf, labels, varargin)
% scatterCompare4Models
% Confronta 4 dataset sullo stesso plot:
%  1) Quad rigida      (marker o)
%  2) Quad flessibile  (marker s)  <-- colorata con CVar
%  3) Semplice rigida  (marker ^)
%  4) Semplice fless.  (marker d)  <-- colorata con CVar
%
% La colorbar e la colormap si applicano SOLO ai due dataset colorati.
%
% INPUT:
%  - T_*          : tabelle
%  - xVar, yVar   : nomi colonne per X e Y
%  - fixed_*      : containers.Map con chiavi=nomi colonne e valori=fissi
%  - labels       : cell array 1x4 {"Quad rigida","Quad elastica","Cost rigida","Cost elastica"}
%
% OPZIONI (Name-Value):
%  'Trend'        "none"|"linear"|"quadratic"  (default "linear")
%  'CVar'         variabile per colore (default "Deflection of the wheel [m]")
%  'Tol'          containers.Map di tolleranze numeriche (default vuota)
%  'ShowData'     true/false (default true)
%  'PointSize'    scalar (default 45)

p = inputParser;
p.addParameter('Trend',"linear", @(s) any(strcmpi(string(s),["none","linear","quadratic"])));
p.addParameter('CVar',"Deflection of the wheel [m]", @(s) isstring(s) || ischar(s));
p.addParameter('Tol', containers.Map(), @(x) isa(x,'containers.Map'));
p.addParameter('ShowData', true, @(b) islogical(b) || isnumeric(b));
p.addParameter('PointSize', 130, @(n) isnumeric(n) && isscalar(n) && n>0);
p.parse(varargin{:});
opt = p.Results;

xVar = string(xVar);
yVar = string(yVar);
cVar = string(opt.CVar);

% ===== 1) Filtra i 4 dataset =====
Q1 = filterTable(T_quad_rigid,  fixed_qr, opt.Tol);
Q2 = filterTable(T_quad_flex,   fixed_qf, opt.Tol);
S1 = filterTable(T_simple_rigid,fixed_sr, opt.Tol);
S2 = filterTable(T_simple_flex, fixed_sf, opt.Tol);

disp("Rows after filtering:")
disp("Q1 = " + height(Q1))
disp("Q2 = " + height(Q2))
disp("S1 = " + height(S1))
disp("S2 = " + height(S2))

% ===== 2) Controlli colonne =====
needXY = [xVar yVar];
assert(all(ismember(needXY, string(T_quad_rigid.Properties.VariableNames))),  "xVar/yVar mancanti in quad rigida");
assert(all(ismember(needXY, string(T_quad_flex.Properties.VariableNames))),   "xVar/yVar mancanti in quad flessibile");
assert(all(ismember(needXY, string(T_simple_rigid.Properties.VariableNames))),"xVar/yVar mancanti in semplice rigida");
assert(all(ismember(needXY, string(T_simple_flex.Properties.VariableNames))), "xVar/yVar mancanti in semplice flessibile");

% ===== 3) Setup figura =====
figure('Name', yVar + " vs " + xVar); hold on; grid on; box on;

% Stili B/N (linee nere, marker diversi)
st(1) = struct('marker','o','line','-');
st(2) = struct('marker','o','line','--');
st(3) = struct('marker','^','line','-');
st(4) = struct('marker','^','line','--');

% ===== 4) Colormap SOLO per dataset 2 e 4 (Q2, S2) =====
hasC2 = ~isempty(Q2) && ismember(cVar, string(Q2.Properties.VariableNames));
hasC4 = ~isempty(S2) && ismember(cVar, string(S2.Properties.VariableNames));

% Range globale del colore solo dai dataset che lo hanno
cAll = [];
if hasC2, cAll = [cAll; Q2.(cVar)]; end
if hasC4, cAll = [cAll; S2.(cVar)]; end

useColor = ~isempty(cAll);
if useColor
    colormap(jet);
    cLim = [min(cAll) max(cAll)];
    if cLim(1) == cLim(2), cLim = cLim + [-1 1]*1e-6; end
    clim(cLim);
end

% ===== 5) Plot dei 4 dataset =====
hLeg = gobjects(1,4);

% (1) Quad rigida (NON colorata)
hLeg(1) = plotModel(Q1, xVar, yVar, st(1), opt.Trend, opt.ShowData, opt.PointSize, ...
                    false, cVar);

% (2) Quad flessibile (colorata se cVar esiste)
hLeg(2) = plotModel(Q2, xVar, yVar, st(2), opt.Trend, opt.ShowData, opt.PointSize, ...
                    true, cVar);

% (3) Semplice rigida (NON colorata)
hLeg(3) = plotModel(S1, xVar, yVar, st(3), opt.Trend, opt.ShowData, opt.PointSize, ...
                    false, cVar);

% (4) Semplice flessibile (colorata se cVar esiste)
hLeg(4) = plotModel(S2, xVar, yVar, st(4), opt.Trend, opt.ShowData, opt.PointSize, ...
                    true, cVar);

xlabel(xVar,'Interpreter','none', 'FontSize',30);
ylabel(yVar,'Interpreter','none', 'FontSize',30);
% title(yVar + " vs " + xVar, 'Interpreter','none');

lgd= legend(hLeg, labels{:}, 'Location','best');
lgd.FontSize = 20;

ax = gca;
ax.FontSize = 25;      % numeri sugli assi più grandi
ax.LineWidth = 2;    % assi un po' più marcati

% Colorbar SOLO se realmente usata (dataset 2 o 4 con cVar presente)
if useColor
    cb = colorbar;
    cb.FontSize = 25;
    ylabel(cb, cVar, 'Interpreter','none','FontSize',25);
end

hold off;

end

% ========================= HELPER FUNCTIONS =========================

function hLegend = plotModel(Ts, xVar, yVar, st, trend, showData, pSize, wantColor, cVar)
% Ritorna un handle "pulito" per legenda.
hLegend = plot(nan,nan,'k','LineStyle',st.line,'Marker',st.marker,'LineWidth',2);

if isempty(Ts) || height(Ts)==0
    return;
end

x = Ts.(xVar);
y = Ts.(yVar);

% punti
if showData
    if wantColor && ismember(cVar, string(Ts.Properties.VariableNames))
        c = Ts.(cVar);
        scatter(x, y, pSize, c, st.marker, ...
            'MarkerFaceColor','flat', 'MarkerEdgeColor','k', 'LineWidth',2.0);
    else
        scatter(x, y, pSize, st.marker, ...
            'MarkerFaceColor','w', 'MarkerEdgeColor','k', 'LineWidth',2.0);
    end
end

% trend (sempre nero)
if strcmpi(trend,"linear") && numel(x) >= 2
    pp = polyfit(x,y,1);
    xx = linspace(min(x),max(x),200);
    plot(xx, polyval(pp,xx), 'k', 'LineStyle',st.line, 'LineWidth',4.0);
elseif strcmpi(trend,"quadratic") && numel(x) >= 3
    pp = polyfit(x,y,2);
    xx = linspace(min(x),max(x),200);
    plot(xx, polyval(pp,xx), 'k', 'LineStyle',st.line, 'LineWidth',4.0);
end
end

function Ts = filterTable(T, fixed, tolMap)
% Filtro righe dove le colonne in fixed assumono il valore richiesto.
mask = true(height(T),1);
keysF = string(keys(fixed));

for k = 1:numel(keysF)
    v = keysF(k);

    if ~ismember(v, string(T.Properties.VariableNames))
        error("Colonna '%s' non presente nella tabella.", v);
    end

    col = T.(v);
    val = fixed(v);

    tol = 0;
    if isa(tolMap,'containers.Map') && isKey(tolMap, v)
        tol = tolMap(v);
    end

    if isnumeric(col)
        if tol == 0
            mask = mask & (col == val);
        else
            mask = mask & (abs(col - val) <= tol);
        end
    else
        mask = mask & (string(col) == string(val));
    end
end

Ts = T(mask,:);

end
