function [fn, zeta, A, Hhat, report] = refine_modes_jitterboost( ...
    f, Hc, fn0, z0, A0, varargin)
% REFINE_MODES_JITTERBOOST
% Multi-start jitter + per-mode micro line-search to squeeze extra accuracy.
% Amplitudes are re-fit (real by default) after every perturbation.
%
% Usage:
%   [fn,zeta,A,Hhat,report] = refine_modes_jitterboost(f,Hc,fn0,z0,A0, ...
%       'FRF',frf,'Model','auto','DCGain',Gdc,'Delay',0,'Areal',true, ...
%       'Rounds',3,'StartsPerRound',30,'FnPct',[1.0 0.6 0.3], ...
%       'ZetaPct',[0.06 0.04 0.02],'LocalLine',true,'LineGrid',7, ...
%       'MinSepPct',3e-4,'Verbose',true)
%
% Output report fields:
%   rmse_after, accepts, evals, rounds, best_round

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('Delay',0,@(x)isnumeric(x)&&isscalar(x));
p.addParameter('Areal',true,@(b)islogical(b)&&isscalar(b));

p.addParameter('Rounds',3,@(x)isnumeric(x)&&isscalar(x)&&x>=1);
p.addParameter('StartsPerRound',30,@(x)isnumeric(x)&&isscalar(x)&&x>=1);
p.addParameter('FnPct',[1.0 0.6 0.3],@isnumeric);     % % jitter of fn
p.addParameter('ZetaPct',[0.06 0.04 0.02],@isnumeric);% absolute fraction jitter of zeta
p.addParameter('LocalLine',true,@(b)islogical(b)&&isscalar(b));
p.addParameter('LineGrid',7,@(x)isnumeric(x)&&isscalar(x)&&x>=3);
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(:);
frf = char(o.FRF);
model = lower(char(o.Model)); if strcmp(model,'auto'), model = iff(isempty(frf),'multiplicative','additive'); end
Gdc = o.DCGain; tau = o.Delay; Areal = o.Areal;
w = 2*pi*f;

% Ensure vectors
fn0 = fn0(:); z0 = z0(:); A0 = A0(:);

% Best-so-far from the provided starting point
if isempty(A0) || any(~isfinite(A0))
    A0 = seed_amplitudes_ls(f, Hc, fn0, z0, frf, Gdc, tau, Areal, model);
end
Hbest = transfer_from_params_local(f, A0, fn0, z0, 'FRF',frf,'Model',model,'DCGain',Gdc,'Delay',tau);
bestErr = complex_rmse(Hc, Hbest);
fn = fn0; zeta = z0; A = A0;

rounds = o.Rounds;
if isscalar(o.FnPct), fnPcts = repmat(o.FnPct,1,rounds); else, fnPcts = o.FnPct; end
if isscalar(o.ZetaPct), zePcts = repmat(o.ZetaPct,1,rounds); else, zePcts = o.ZetaPct; end
fnPcts = fnPcts(1:rounds); zePcts = zePcts(1:rounds);

accepts = 0; evals = 0; best_round = 0;
rng(2); % deterministic

for r = 1:rounds
    fn_pct   = fnPcts(r);          % percent for frequencies
    zeta_pct = zePcts(r);          % fraction for zeta
    roundBestErr = bestErr; roundBest = struct('fn',fn,'zeta',zeta,'A',A,'H',Hbest);

    for s = 1:o.StartsPerRound
        % --- random jitter proposal ---
        fn_try   = fn .* (1 + (fn_pct/100) * 0.8*randn(size(fn)));
        z_try    = zeta .* (1 + zeta_pct * 0.8*randn(size(zeta)));
        z_try    = min(max(z_try, 1e-4), 0.30);

        % enforce bounds and separation
        fn_try = min(max(fn_try, max(f(2), 0.5*min(f))), min(f(end-1), 2.0*max(f)));
        [fn_try,ord] = sort(fn_try(:)); z_try = z_try(ord);
        for k=2:numel(fn_try)
            if fn_try(k) <= fn_try(k-1)*(1+o.MinSepPct)
                fn_try(k) = fn_try(k-1)*(1+o.MinSepPct);
            end
        end

        % re-fit A
        A_try = seed_amplitudes_ls(f, Hc, fn_try, z_try, frf, Gdc, tau, Areal, model);

        % optional per-mode micro line search on frequency (keeps zeta fixed)
        if o.LocalLine
            for k = 1:numel(fn_try)
                half = 0.5*(fn_pct/100)*fn_try(k);
                if half <= 0, continue; end
                fgrid = linspace(max(f(2), fn_try(k)-half), min(f(end-1), fn_try(k)+half), o.LineGrid);
                best_k_err = inf; best_fk = fn_try(k); best_Ak = A_try;
                for fk = fgrid
                    fn_tmp      = fn_try; fn_tmp(k) = fk;
                    A_tmp       = seed_amplitudes_ls(f, Hc, fn_tmp, z_try, frf, Gdc, tau, Areal, model);
                    H_tmp       = transfer_from_params_local(f, A_tmp, fn_tmp, z_try, 'FRF',frf,'Model',model,'DCGain',Gdc,'Delay',tau);
                    err_tmp     = complex_rmse(Hc, H_tmp);
                    evals = evals + 1;
                    if err_tmp < best_k_err
                        best_k_err = err_tmp; best_fk = fk; best_Ak = A_tmp;
                    end
                end
                fn_try = fn_try; fn_try(k) = best_fk; A_try = best_Ak;
            end
        end

        H_try   = transfer_from_params_local(f, A_try, fn_try, z_try, 'FRF',frf,'Model',model,'DCGain',Gdc,'Delay',tau);
        err_try = complex_rmse(Hc, H_try);
        evals   = evals + 1;

        if err_try < roundBestErr - 1e-12
            roundBestErr = err_try;
            roundBest = struct('fn',fn_try,'zeta',z_try,'A',A_try,'H',H_try);
        end
    end

    % accept round winner if improved over global best
    if roundBestErr < bestErr - 1e-12
        fn = roundBest.fn; zeta = roundBest.zeta; A = roundBest.A;
        Hbest = roundBest.H; bestErr = roundBestErr;
        accepts = accepts + 1; best_round = r;
    end

    if o.Verbose
        fprintf('[jitterboost] round %d/%d: best rmse=%.6g, global=%.6g\n', ...
            r, rounds, roundBestErr, bestErr);
    end
end

% light coordinate polish (greedy, very small steps)
fn_step_pct = max(0.25, fnPcts(end)*0.5);
ze_step_pct = max(0.015, zePcts(end)*0.5);
improved = true; sweep = 0;
while improved && sweep < 3
    improved = false; sweep = sweep + 1;
    for k=1:numel(fn)
        % try ± small fn step
        for sgn = [-1 1]
            fk = fn(k) * (1 + sgn*fn_step_pct/100);
            fk = min(max(fk, f(2)), f(end-1));
            fn_tmp = fn; fn_tmp(k) = fk;
            for j=2:numel(fn_tmp)
                if fn_tmp(j) <= fn_tmp(j-1)*(1+o.MinSepPct), fn_tmp(j)=fn_tmp(j-1)*(1+o.MinSepPct); end
            end
            A_tmp = seed_amplitudes_ls(f, Hc, fn_tmp, zeta, frf, Gdc, tau, Areal, model);
            H_tmp = transfer_from_params_local(f, A_tmp, fn_tmp, zeta, 'FRF',frf,'Model',model,'DCGain',Gdc,'Delay',tau);
            err   = complex_rmse(Hc, H_tmp);
            evals = evals + 1;
            if err < bestErr - 1e-12
                fn = fn_tmp; A = A_tmp; Hbest = H_tmp; bestErr = err; improved = true;
            end
        end
        % try ± small zeta step
        for sgn = [-1 1]
            zk = min(max(zeta(k) * (1 + sgn*ze_step_pct), 1e-4), 0.30);
            z_tmp = zeta; z_tmp(k) = zk;
            A_tmp = seed_amplitudes_ls(f, Hc, fn, z_tmp, frf, Gdc, tau, Areal, model);
            H_tmp = transfer_from_params_local(f, A_tmp, fn, z_tmp, 'FRF',frf,'Model',model,'DCGain',Gdc,'Delay',tau);
            err   = complex_rmse(Hc, H_tmp);
            evals = evals + 1;
            if err < bestErr - 1e-12
                zeta = z_tmp; A = A_tmp; Hbest = H_tmp; bestErr = err; improved = true;
            end
        end
    end
end

Hhat = Hbest;
report = struct('rmse_after', bestErr, 'accepts', accepts, ...
                'evals', evals, 'rounds', rounds, 'best_round', best_round);
end

% --- tiny helper (local) ---
function y = iff(c,a,b), if c, y=a; else, y=b; end, end
