function [fn, zeta, A, Hhat, report] = refine_modes_levenberg( ...
    f, Hc, fn0, z0, A0, varargin)

% ==================== OPTIONS ====================
p = inputParser;
p.addParameter('FRF','',@(s)ischar(s)||isstring(s));
p.addParameter('Model','auto',@(s)ischar(s)||isstring(s));
p.addParameter('DCGain',1+0i,@(x)isnumeric(x));
p.addParameter('Areal',true,@(b)islogical(b)&&isscalar(b));
p.addParameter('Delay',0,@(x)isnumeric(x)&&isscalar(x));
p.addParameter('EstimateDelay',false,@(b)islogical(b)&&isscalar(b));
p.addParameter('FitBaseline',false,@(b)islogical(b)&&isscalar(b));
p.addParameter('MaxIters',120,@(x)isnumeric(x)&&isscalar(x)&&x>=1);
p.addParameter('Lambda0',1e-2,@(x)isnumeric(x)&&isscalar(x)&&x>0);
p.addParameter('MinSepPct',3e-4,@(x)isnumeric(x)&&isscalar(x));
p.addParameter('Verbose',true,@(b)islogical(b)&&isscalar(b));
p.parse(varargin{:});
o = p.Results;

f=f(:); Hc=Hc(:);

% ==================== NORMALIZATION (NEW) ====================
scale = max(abs(Hc));
if scale == 0, scale = 1; end
Hc = Hc / scale;

w=2*pi*f;
frf   = char(o.FRF);
model = lower(char(o.Model));
if strcmp(model,'auto')
    model = iff(isempty(frf),'multiplicative','additive');
end
Gdc   = o.DCGain / scale;
Areal=o.Areal;

% ==================== START PARAMS ====================
fn   = fn0(:);
zeta = z0(:);
A    = A0(:) / scale;

if isempty(A) || any(~isfinite(A))
    A = seed_amplitudes_ls(f,Hc,fn,zeta,frf,Gdc,o.Delay,Areal,model);
end

tau  = o.Delay;
dGr  = 0; dGi = 0;

params = pack(fn,zeta,A,tau,dGr,dGi);

r0   = objective(params);
best = struct('params',params,'err',sqrt(mean(r0.^2)));

lambda = o.Lambda0;
accepts=0;

% ==================== LM LOOP ====================
for it=1:o.MaxIters

    [r, J] = residual_and_jac(params);

    Hn = J.'*J;
    g  = J.'*r;

    step = -(Hn + lambda*diag(diag(Hn))) \ g;

    trial = project(apply_step(params, step));

    r_try = objective(trial);
    err_try = sqrt(mean(r_try.^2));

    if err_try < best.err - 1e-12
        params = trial;
        best.params = trial;
        best.err = err_try;
        lambda = max(lambda*0.5, 1e-8);
        accepts = accepts + 1;
    else
        lambda = min(lambda*3, 1e8);
    end
end

[fn, zeta, A, tau, dGr, dGi] = unpack(best.params);

Hhat = compose(best.params) * scale;  % rescale output
A    = A * scale;

report = struct('rmse_after', complex_rmse(Hc*scale, Hhat), ...
                'iters', o.MaxIters, ...
                'accepts', accepts, ...
                'lambda', lambda, ...
                'tau', tau, ...
                'dG', dGr+1i*dGi, ...
                'fitBaseline', o.FitBaseline, ...
                'model', model, ...
                'frf', frf);

% ==============================================================
% ==================== NESTED FUNCTIONS ========================
% ==============================================================

    function r = objective(p)
        H = compose(p);
        rC = Hc - H;
        r  = [real(rC); imag(rC)];
    end

    function [r, J] = residual_and_jac(p)

        [fnP, zP, AP, tauP, dGrP, dGiP] = unpack(p);
        N = numel(f);
        M = numel(fnP);

        Hcols  = zeros(N, M);
        dH_dfn = zeros(N, M);
        dH_dze = zeros(N, M);

        for k=1:M
            wn = 2*pi*fnP(k);
            z  = zP(k);
            [Hk, dHk_dwn, dHk_dz] = kernel_and_derivs(w, wn, z, frf);
            Hcols(:,k) = Hk;
            dH_dfn(:,k) = dHk_dwn * (2*pi);
            dH_dze(:,k) = dHk_dz;
        end

        if strcmp(model,'multiplicative')
            H = Gdc * (1 + Hcols*AP);
        else
            H = Gdc + Hcols*AP + (o.FitBaseline)*(dGrP + 1i*dGiP);
        end

        rC = Hc - H;
        r  = [real(rC); imag(rC)];

        JA  = [-real(Hcols); -imag(Hcols)];
        Jfn = [-real(dH_dfn); -imag(dH_dfn)];
        Jze = [-real(dH_dze); -imag(dH_dze)];

        J = [Jfn, Jze, JA];
    end

    function p2 = apply_step(p, s)
        [fnP, zP, AP, tauP, dGrP, dGiP] = unpack(p);
        nf=numel(fnP); off=0;
        fnP = fnP + s(off+(1:nf)); off=off+nf;
        zP  = zP  + s(off+(1:nf)); off=off+nf;
        AP  = AP  + s(off+(1:nf));
        p2 = pack(fnP,zP,AP,tauP,dGrP,dGiP);
    end

    function pP = project(p)

        [fnP, zP, AP, tauP, dGrP, dGiP] = unpack(p);

        % ==================== STABILIZERS ====================

        % 1) Constrain frequency to measured band
        fmin = min(f);
        fmax = max(f);
        fnP  = max(fnP, fmin*0.98);
        fnP  = min(fnP, fmax*1.02);

        % 2) Enforce damping floor
        zP   = max(zP, 0.002);    % 0.2% minimum damping
        zP   = min(zP, 0.25);

        % 3) Enforce separation
        [fnP, ord] = sort(fnP(:));
        zP = zP(ord);
        AP = AP(ord);

        for k=2:numel(fnP)
            if fnP(k) <= fnP(k-1)*(1+o.MinSepPct)
                fnP(k) = fnP(k-1)*(1+o.MinSepPct);
            end
        end

        pP = pack(fnP,zP,AP,tauP,dGrP,dGiP);
    end

    function H = compose(p)
        [fnP, zP, AP, tauP, dGrP, dGiP] = unpack(p);
        N = numel(f);
        M = numel(fnP);
        Hcols = zeros(N, M);
        for k=1:M
            Hcols(:,k) = frf_kernel(w, 2*pi*fnP(k), zP(k), frf);
        end
        if strcmp(model,'multiplicative')
            H = Gdc * (1 + Hcols*AP);
        else
            H = Gdc + Hcols*AP + (o.FitBaseline)*(dGrP + 1i*dGiP);
        end
    end

    function p = pack(fnP,zP,AP,tauP,dGrP,dGiP)
        p = [fnP(:); zP(:); AP(:)];
    end

    function [fnP,zP,AP,tauP,dGrP,dGiP] = unpack(v)
        nf = numel(fn);
        fnP = v(1:nf);
        zP  = v(nf+1:2*nf);
        AP  = v(2*nf+1:3*nf);
        tauP=0; dGrP=0; dGiP=0;
    end
end

% ==============================================================
% ==================== KERNEL DERIVATIVES ======================
% ==============================================================

function [Hk, dHk_dwn, dHk_dz] = kernel_and_derivs(w, wn, zeta, frfType)

den = (wn.^2 - w.^2) + 1i*(2*zeta*wn).*w;

switch lower(char(frfType))
    case 'disp/force', N = 1;
    case 'vel/force',  N = 1i*w;
    case 'acc/force',  N = -(w.^2);
    otherwise,         N = (w.^2);
end

Hk = N ./ den;

dden_dwn = 2*wn + 1i*(2*zeta).*w;
dHk_dwn  = -(N .* dden_dwn) ./ (den.^2);

dden_dz  = 1i*(2*wn).*w;
dHk_dz   = -(N .* dden_dz) ./ (den.^2);
end

function rm = complex_rmse(a,b)
r=a(:)-b(:);
rm=sqrt(mean(abs(r).^2));
end

function y  = iff(c,a,b)
if c, y=a; else, y=b; end
end