clc; clear; close all;

%% ===================== LOAD TABLE =====================
load('TAB_tol_twin.mat');
T = TAB_tol_twin;

%% ===================== COLUMN NAMES =====================
col_t_ext   = 't ext [m]';
col_t_int   = 't int [m]';
col_t_rinf  = 't rinf [m]';

col_stress  = 'Maximum stress [Pa]';
col_stiff   = 'Wheel stiffness [kg/s^2]';
col_dur     = 'Durability index';

t_nom = 0.003;   % nominal thickness

%% ===================== FIGURE =====================
figure('Color','w','Position',[100 100 950 650]);
tiledlayout(3,1,'TileSpacing','compact','Padding','compact');

thicknessCols   = {col_t_ext, col_t_int, col_t_rinf};
thicknessLabels = {'Outer ring thickness','Inner ring thickness','Reinforcement thickness'};

for k = 1:3
    variedCol = thicknessCols{k};
    otherCols = thicknessCols;
    otherCols(k) = [];

    % === Trova tutte le combinazioni possibili delle altre due thickness ===
    G = findgroups(T.(otherCols{1}), T.(otherCols{2}));
    groups = unique(G);

    bestBlock = [];
    bestScore = -inf;

    for g = groups'
        Tb_try = T(G == g, :);

        % Deve esserci variazione solo nella colonna selezionata
        nVarMain = numel(unique(Tb_try.(variedCol)));
        nVarOth1 = numel(unique(Tb_try.(otherCols{1})));
        nVarOth2 = numel(unique(Tb_try.(otherCols{2})));

        if nVarMain > 1 && nVarOth1 == 1 && nVarOth2 == 1
            % punteggio: privilegia il blocco con più valori unici
            % e con gli altri due thickness più vicini al nominale
            fixed1 = unique(Tb_try.(otherCols{1}));
            fixed2 = unique(Tb_try.(otherCols{2}));

            score = 100*nVarMain - abs(fixed1 - t_nom) - abs(fixed2 - t_nom);

            if score > bestScore
                bestScore = score;
                bestBlock = Tb_try;
            end
        end
    end

    nexttile

    if isempty(bestBlock)
        axis off
        text(0.5,0.5,['No valid block found for ', thicknessLabels{k}], ...
            'HorizontalAlignment','center','FontSize',14)
        continue
    end

    %% ===================== CLEAN BLOCK =====================
    Tb = sortrows(bestBlock, variedCol);
    x_raw = Tb.(variedCol);

    % Se ci sono duplicati dello stesso spessore, media i risultati
    [x, ~, ic] = unique(x_raw);

    stress_v = accumarray(ic, Tb.(col_stress), [], @mean);
    stiff_v  = accumarray(ic, Tb.(col_stiff),  [], @mean);
    dur_v    = accumarray(ic, Tb.(col_dur),    [], @mean);

    % Trova il riferimento nominale
    idx0 = find(abs(x - t_nom) < 1e-12, 1);

    if isempty(idx0)
        % Se il nominale non esiste, usa il punto centrale
        [~, idx0] = min(abs(x - median(x)));
    end

    x0 = x(idx0);

    stress0 = stress_v(idx0);
    stiff0  = stiff_v(idx0);
    dur0    = dur_v(idx0);

    dx_pct = 100 * (x - x0) / x0;

    dStress_pct = 100 * (stress_v - stress0) / stress0;
    dStiff_pct  = 100 * (stiff_v  - stiff0 ) / stiff0;
    dDur_pct    = 100 * (dur_v    - dur0   ) / dur0;

    % ===================== ROBUSTNESS INDEX =====================
    R = ( max(abs(dStress_pct)) + max(abs(dStiff_pct)) + max(abs(dDur_pct)) ) / 3;

    % Conversione in score 1-5
    if R <= 5
        robustness_score = 5;
    elseif R <= 10
        robustness_score = 4;
    elseif R <= 20
        robustness_score = 3;
    elseif R <= 30
        robustness_score = 2;
    else
        robustness_score = 1;
    end

    fprintf('\n%s\n', thicknessLabels{k});
    fprintf('Robustness index R = %.2f %%\n', R);
    fprintf('Robustness score   = %d / 5\n', robustness_score);


    %% ===================== PLOT =====================
    hold on
    grid on
    box on

    plot(dx_pct, dStress_pct, '-o', 'LineWidth', 2.5, 'MarkerSize', 6)
    plot(dx_pct, dStiff_pct,  '-s', 'LineWidth', 2.5, 'MarkerSize', 6)
    plot(dx_pct, dDur_pct,    '-d', 'LineWidth', 2.5, 'MarkerSize', 6)

    xline(0,'--k','LineWidth',1.5);
    yline(0,'--k','LineWidth',1.5);

    ylabel('Response [\%]','Interpreter','latex','FontSize',15)
    title(thicknessLabels{k}, 'FontWeight','bold','Interpreter','latex','FontSize',18)

    if k == 3
        xlabel('Thickness variation [\%]','Interpreter','latex','FontSize',20)
    end

    legend('Stress','Stiffness','Durability', ...
        'Location','northwest','Interpreter','latex','FontSize',15)

    ax = gca;
    ax.FontSize = 14;
    ax.LineWidth = 1;
    ax.TickLabelInterpreter = 'latex';

    hold off
end