function [Tfit_plot, Ahat] = reconstruct_fit_from_fn_zeta(f, Tmag_in, fn, zeta, iat)
% Rebuild |T(f)| using fixed fn, zeta; solve only for amplitudes A >= 0.
% Inputs:
%   f        : frequency vector (Hz)
%   Tmag_in  : measured magnitude (linear units if iat==1, power if iat==2)
%   fn,zeta  : column vectors (same length = nmodes)
%   iat      : 1 -> amplitude domain; 2 -> power domain
% Outputs:
%   Tfit_plot: fitted magnitude in the SAME domain as Tmag_in
%   Ahat     : fitted modal amplitudes (amplitude domain)

    f = f(:); Tmag_in = Tmag_in(:);
    nm = numel(fn);

    % Work internally in amplitude domain; convert if input is power
    if iat==2
        y = sqrt(max(Tmag_in, 0));   % amplitude-equivalent
    else
        y = Tmag_in;
    end

    % Build complex modal basis G(:,i) = H_i(f)
    G = zeros(numel(f), nm);
    r = f ./ fn(:).';                 % broadcast: (#f x #modes)
    den = (1 - r.^2) + 1j*(2*zeta(:).'.*r);
    G = 1 ./ den;                     % each column is mode i

    % Initial guess: amplitude from peak rule A ≈ 2ζ * peak(|T|)
    A0 = 2*zeta(:).*max(y); A0 = max(A0, 1e-12);

    % Solve A >= 0 to match |G*A| to y (nonlinear least squares on logA)
    u0 = log(max(A0,1e-12));
    function err = obj(u)
        A = exp(u);
        yh = abs(G * A);
        w  = 1./max(y,1e-6);          % de-emphasize large peaks dominating
        err = sum(((yh - y).*w).^2);
    end
    opts = optimset('fminsearch'); opts.Display='off'; opts.MaxIter=2000; opts.MaxFunEvals=8000;
    uhat = fminsearch(@obj, u0, opts);
    Ahat = exp(uhat);

    % Final fit (amplitude domain), then convert back for plotting if needed
    Tfit_amp = abs(G * Ahat);
    if iat==2
        Tfit_plot = Tfit_amp.^2;      % back to power to match your plot
    else
        Tfit_plot = Tfit_amp;
    end
end
