function [fnr, dampr, Ar] = trans_mag_extract_40(THM, ndof, iat)

    % Keep original copy ONLY for reference if needed later
    THM_o = THM; %#ok<NASGU>

    % If iat==2, input magnitude is POWER (|T|^2)
    % Convert to amplitude before fitting
    if iat==2 && size(THM,2)>=2
        THM(:,2) = sqrt(THM(:,2));
    end

    % ---------- Accept [f |T|] or [f Re Im] ----------
    if size(THM,2) == 2
        f    = THM(:,1);
        Tmag = THM(:,2);
    else
        f    = THM(:,1);
        Hc   = THM(:,2) + 1j*THM(:,3);
        Tmag = abs(Hc);
    end

    % ---------- Sanitize ----------
    good = isfinite(f) & isfinite(Tmag) & (Tmag > 0);
    f    = f(good);
    Tmag = Tmag(good);

    % ---------- Fit ----------
    [params, Tfit, ~] = trans_mag_pickfit([f, Tmag], ndof, ...
        'NumModes',           ndof, ...
        'DoParabolicPolish',  true, ...
        'DoHalfPowerTouchup', true, ...
        'DoLocalRefine',      false, ...
        'DoGlobalRefine',     false, ...
        'UseResidualMode',    true, ...
        'ShowProgress',       false, ...
        'Plot',               false, ...
        'MaxIters',           200, ...
        'Display',            'off', ...
        'DoSKPolish',         false);

    % ---------- Unpack + sort ----------
    fnr   = params(:,1);
    dampr = params(:,2);
    Ar    = params(:,3);

    [fnr, ord] = sort(fnr(:));
    dampr      = dampr(ord);
    Ar         = Ar(ord);

    % ---------- IMPORTANT FIX ----------
    % Preserve the MEASURED data (do NOT overwrite with fit)
    Tmag_input = Tmag;

    % ---------- Reconstruct fit from (fn, zeta) ----------
    [Tfit_plot, ~] = reconstruct_fit_from_fn_zeta( ...
        f, Tmag_input, fnr, dampr, 1);

    % ---------- Domain ----------
    fmin = min(f);
    fmax = max(f);

    % ---------- Amplitude vs Power handling ----------
    if iat == 1
        y_label  = 'Trans (G/G)';
        t_string = 'Transmissibility';
        md = 5;
    else
        y_label  = 'Trans (G^2/G^2)';
        t_string = 'Power Transmissibility';
        Tfit_plot  = Tfit_plot.^2;
        Tmag_input = Tmag_input.^2;
        md = 6;
    end

    % ---------- Plot data ----------
    ppp1 = [f, Tfit_plot];    % curve-fit
    ppp2 = [f, Tmag_input];   % MEASURED input

    x_label = 'Frequency (Hz)';
    fig_num = 10;

    leg1 = 'Curve-fit';
    leg2 = 'Input';

    plot_loglog_function_md_two_h2( ...
        fig_num, x_label, y_label, t_string, ...
        ppp1, ppp2, leg1, leg2, fmin, fmax, md);

    % ---------- Console summary ----------
    q = [Ar fnr dampr];
    q = sortrows(q, 2);

    fprintf('\n\n  A        fn(Hz)    damp\n');
    for i = 1:size(q,1)
        fprintf(' %8.5g  %8.3f  %8.5f\n', q(i,1), q(i,2), q(i,3));
    end
end
