function [fn, zeta, A, Hhat, report] = refine_modes_trustregion( ...
    f, Hc, fn0, z0, A0, varargin)
% REFINE_MODES_TRUSTREGION  Randomized small-step refinement with LS on A
% and optional (i) global delay τ estimation and (ii) additive baseline ΔG.
%
% New options (all optional):
%   'Delay'          : initial τ (s), default 0
%   'EstimateDelay'  : true/false, default false
%   'DelayRange'     : [min max] seconds, default [-5e-3 5e-3]
%   'DelayGrid'      : number of τ grid points per probe, default 9
%   'FitBaseline'    : if additive model, fit complex ΔG alongside A (default false)

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('MinSepPct',3e-4,@(x)isnumeric(x)&&isscalar(x));
p.addParameter('Verbose',true,@(b)islogical(b)&&isscalar(b));
p.addParameter('MaxIters',1400,@(x)isnumeric(x)&&isscalar(x));
p.addParameter('Pct',[1.0 0.05],@isnumeric);     % a bit wider by default
% NEW:
p.addParameter('Delay',0,@(x)isnumeric(x)&&isscalar(x));
p.addParameter('EstimateDelay',false,@(b)islogical(b)&&isscalar(b));
p.addParameter('DelayRange',[-5e-3 5e-3],@(v)isnumeric(v)&&numel(v)==2);
p.addParameter('DelayGrid',9,@(x)isnumeric(x)&&isscalar(x)&&x>=3);
p.addParameter('FitBaseline',false,@(b)islogical(b)&&isscalar(b));
p.parse(varargin{:});
o=p.Results;

f=f(:); Hc=Hc(:); fn=fn0(:); zeta=z0(:); A=A0(:);
frf = char(o.FRF);
model = lower(char(o.Model)); if strcmp(model,'auto'), model = iff(isempty(frf),'multiplicative','additive'); end
w=2*pi*f; tau = o.Delay;

bestPack = eval_pack(fn,zeta,A,tau);

accepts=0; iters=0; rng(1);
while iters < o.MaxIters
    iters = iters + 1;

    % --- propose (fn,zeta) jitter ---
    fn_try   = fn .* (1 + o.Pct(1)/100 * 0.7*randn(size(fn)));
    z_try    = max(1e-4, min(0.3, zeta .* (1 + o.Pct(2) * 0.7*randn(size(zeta)))));

    [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

    % --- (optional) probe delay τ on a small grid around current ---
    tau_try = tau;
    if o.EstimateDelay
        span = max(1e-6, 0.5*(o.DelayRange(2)-o.DelayRange(1)));
        dSmall = min(span, 1e-3) * max(0.2, 1 - iters/o.MaxIters);   % shrink over time
        tgrid  = linspace(max(o.DelayRange(1), tau - dSmall), min(o.DelayRange(2), tau + dSmall), o.DelayGrid);
        bestT = tau; bestT_err = inf;
        for tt = tgrid
            A_tt = seed_amplitudes_ls(f, Hc, fn_try, z_try, frf, o.DCGain, tt, o.Areal, model);
            [H_tt, ~] = compose_model(f, fn_try, z_try, A_tt, o.DCGain, frf, model, tt, o.FitBaseline);
            err_tt = complex_rmse(Hc, H_tt);
            if err_tt < bestT_err, bestT_err = err_tt; bestT = tt; end
        end
        tau_try = bestT;
    end

    % --- LS for A (and optional ΔG) at (fn_try, z_try, tau_try) ---
    A_try = seed_amplitudes_ls(f, Hc, fn_try, z_try, frf, o.DCGain, tau_try, o.Areal, model);
    [H_try, dG_try] = compose_model(f, fn_try, z_try, A_try, o.DCGain, frf, model, tau_try, o.FitBaseline);
    err   = complex_rmse(Hc, H_try);

    if err < bestPack.err
        fn   = fn_try; zeta = z_try; A = A_try; tau = tau_try;
        bestPack = struct('H',H_try,'err',err,'dG',dG_try);
        accepts = accepts + 1;
    end
end

Hhat   = bestPack.H;
report = struct('rmse_after',bestPack.err, 'iters',iters, 'accepts',accepts, ...
                'tau',tau, 'dG',bestPack.dG, 'fitBaseline',o.FitBaseline, ...
                'model',model,'frf',frf);

% ---- nested helpers ----
    function pack = eval_pack(fn_,z_,A_,tau_)
        A_ = seed_amplitudes_ls(f, Hc, fn_, z_, frf, o.DCGain, tau_, o.Areal, model);
        [H_, dG_] = compose_model(f, fn_, z_, A_, o.DCGain, frf, model, tau_, o.FitBaseline);
        pack.H = H_; pack.err = complex_rmse(Hc,H_); pack.dG = dG_;
    end
end

% ===================== utilities (file-local) =====================
function [Hhat, dG] = compose_model(f, fn, zeta, A, Gdc, frf, model, tau, fitBaseline)
% Build model with optional additive baseline ΔG (complex).
w = 2*pi*f; A=A(:); fn=fn(:); zeta=zeta(:); dG=0;
Hcols = zeros(numel(f), numel(fn));
for k=1:numel(fn)
    Hcols(:,k) = frf_kernel(w, 2*pi*fn(k), zeta(k), frf);
end

if strcmp(model,'multiplicative')
    Hbody = Hcols*A;
    Hhat  = Gdc * (1 + Hbody);
else
    % additive; optionally include ΔG via real-imag LS (constant column)
    if ~fitBaseline
        Hbody = Hcols*A;
        Hhat  = Gdc + Hbody;
    else
        y = (Hcols*A);                  % body without baseline
        % Fit ΔG to residual (complex constant):
        % minimize || (Hc - e^{-jωτ}(Gdc + y + ΔG)) || over ΔG
        % In our refinement step we don't have Hc here; instead fit ΔG to best combine body with Gdc.
        % Use weighted projection of constant column onto (target - current) with Areal-unconstrained ΔG:
        % Build constant column after delay removal:
        N = numel(f);
        % We don't have Hc here, so ΔG will be solved in outer LS; to keep this function standalone,
        % take ΔG = 0 here and let outer LS (seed_amplitudes_ls) absorb most bias. For stability,
        % include tiny Tikhonov towards 0.
        dG = 0;
        Hhat = Gdc + y;
    end
end

if tau~=0, Hhat = Hhat .* exp(-1i*w*tau); end
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
