function [fn, zeta, A, Hhat, info] = mdof_trans_main(THM, ndof, iat, imiss, ibh, varargin)

fprintf('\n\n mdof_trans_main \n\n');

% =====================================================================
% MDOF_TRANS_MAIN  (Transmissibility + Force-FRF pipeline)
%
% CLEAN ARCHITECTURE (as requested):
%
%   iat = 1  → Absolute acceleration transmissibility  u¨/y¨   (dimensionless)
%              - Identification pipeline unchanged (keeps your accurate fn, zeta)
%              - FINAL OUTPUT/PLOTTING uses ANALYTICAL multiplicative physics
%
%   iat = 2  → acc/force
%   iat = 3  → vel/force
%   iat = 4  → disp/force
%
% Force FRFs use residue-based additive model for BOTH fitting and plotting.
%
% Input THM (supports three forms):
%   1) Single FRF (matrix):     [f Re Im] or [f H]
%   2) Multi FRFs (cell):       {THM1, THM2, ...} all on same f
%   3) Multi FRFs (matrix):     [f, H1, H2, ...] with each Hk COMPLEX column
%
% Options (subset):
%   'PrimaryChan'              which FRF to use as the SISO fit channel (default: 1)
%   'Plot','Verbose','ForceAreal','EstimateDelay','DelayRange','DelayGrid',
%   'FitBaseline','MinSepPct','UseJitterBoost','JB_Rounds','JB_Starts',
%   'JB_FnPct','JB_ZetaPct','JB_LineGrid','ShowProgress','MaxSeconds'
%
%   --- Speed knobs added ---
%   'SkipBasinIfLMGainBelow'
%   'BasinTrials','BasinMinTrials','BasinAutoScale','BasinDecimate'
%   'BasinDecimTargetN','BasinFocusWinFactor','BasinUseParallel'
%   'BasinSobol','BasinEarlyStopPatience','BasinRefineMaxIters','BasinSeed'
%
% NOTE:
%   This file includes all LOCAL helper functions used here.
%   It still depends on your existing project functions:
%     estimate_modes_halfpower_40, preidentify_modes_complex_local,
%     seed_amplitudes_ls, find_missing_modes_by_ratio_addonly,
%     find_missing_overlapped_modes_splitrefine, refine_modes_trustregion,
%     refine_modes_jitterboost, resolve_close_modes_merge_if_needed,
%     refine_modes_levenberg, resolve_close_modes_bic_cmif,
%     collapse_duplicate_modes, refine_modes_basin_mc_fast/refine_modes_basin_mc
%   Those are unchanged by this plotting fix.
% =====================================================================

% -------------------- options --------------------
p = inputParser;
p.addParameter('Plot', true, @(b)islogical(b)&&isscalar(b));
p.addParameter('Verbose', true, @(b)islogical(b)&&isscalar(b));
p.addParameter('ForceAreal', true, @(b)islogical(b)&&isscalar(b));
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', [], @(x)islogical(x)||isempty(x));
p.addParameter('MinSepPct', 3e-4, @(x)isnumeric(x)&&isscalar(x));
p.addParameter('UseJitterBoost', true, @(b)islogical(b)&&isscalar(b));
p.addParameter('JB_Rounds', 3, @(x)isnumeric(x)&&isscalar(x));
p.addParameter('JB_Starts', 30, @(x)isnumeric(x)&&isscalar(x));
p.addParameter('JB_FnPct', [1.0 0.6 0.3], @isnumeric);
p.addParameter('JB_ZetaPct', [0.06 0.04 0.02], @isnumeric);
p.addParameter('JB_LineGrid', 7, @(x)isnumeric(x)&&isscalar(x)&&x>=3);
p.addParameter('ShowProgress', true, @(b)islogical(b)&&isscalar(b));
p.addParameter('MaxSeconds', inf, @(x)isnumeric(x)&&isscalar(x)&&x>0);
p.addParameter('PrimaryChan', 1, @(x)isnumeric(x)&&isscalar(x)&&x>=1);

% speed knobs
p.addParameter('SkipBasinIfLMGainBelow', 1.0e-3, @(x)isnumeric(x)&&isscalar(x)&&x>=0);
p.addParameter('BasinTrials', 1200, @(x)isnumeric(x)&&isscalar(x)&&x>=0);
p.addParameter('BasinMinTrials', 200, @(x)isnumeric(x)&&isscalar(x)&&x>=0);
p.addParameter('BasinAutoScale', true, @(b)islogical(b)&&isscalar(b));
p.addParameter('BasinDecimate', true, @(b)islogical(b)&&isscalar(b));
p.addParameter('BasinDecimTargetN', [], @(x)isempty(x)||(isnumeric(x)&&isscalar(x)&&x>=64));
p.addParameter('BasinFocusWinFactor', 1.75, @(x)isnumeric(x)&&isscalar(x)&&x>0);
p.addParameter('BasinUseParallel', false, @(b)islogical(b)&&isscalar(b)); % ignored by fast impl
p.addParameter('BasinSobol', true, @(b)islogical(b)&&isscalar(b));
p.addParameter('BasinEarlyStopPatience', [], @(x)isempty(x)||(isnumeric(x)&&isscalar(x)&&x>=10));
p.addParameter('BasinRefineMaxIters', 40, @(x)isnumeric(x)&&isscalar(x)&&x>=5);
p.addParameter('BasinSeed', []);
p.parse(varargin{:});
opt = p.Results;

t0 = tic; timedout = @() toc(t0) > opt.MaxSeconds;

% -------------------- progress bar --------------------
hwb = []; prog = 0;
weights = struct('init',0.12,'missing',0.15,'overlap',0.15,'trust',0.15, ...
                 'jitter',0.12,'resolver1',0.06,'levenberg',0.10,'bic',0.10,'collapse',0.05);
if opt.ShowProgress
    try
        hwb = waitbar(0,'Initializing…','Name','mdof_trans_main','WindowStyle','modal','Visible','on');
    catch, hwb = []; end %#ok<CTCH>
end
cleanupObj = onCleanup(@() try_close_waitbar(hwb)); %#ok<NASGU>
upd = @(frac,msg) update_waitbar_local(hwb, frac, msg);

% -------------------- data & FRF choice --------------------
Ny = 1; ich = 1; Hmat_all = []; % Ny x Nf

if iscell(THM)
    Ny = numel(THM);
    Hrows = cell(1,Ny);
    for j = 1:Ny
        if exist('thm_to_complex','file') == 2
            [fj, Hj] = thm_to_complex(THM{j});
        else
            [fj, Hj] = thm_to_complex_local(THM{j});
        end
        if j==1
            f = fj(:);
        else
            if numel(fj)~=numel(f) || max(abs(fj(:)-f))>eps(max(f))
                error('All FRFs must share the same frequency vector.');
            end
        end
        Hrows{j} = Hj(:).';
    end
    Hmat_all = cat(1, Hrows{:});              % Ny x Nf
    ich = max(1, min(opt.PrimaryChan, Ny));
    Hc = Hmat_all(ich,:).';                   % Nf x 1
elseif isnumeric(THM)
    % Numeric single or multi
    if size(THM,2) >= 3 && any(any(~isreal(THM(:,2:end))))
        % Multi-response numeric with complex columns
        f = double(THM(:,1));
        Hmulti = THM(:,2:end);
        Ny = size(Hmulti,2);
        ich = max(1, min(opt.PrimaryChan, Ny));
        Hc = Hmulti(:,ich);
        Hmat_all = Hmulti.';                  % Ny x Nf
    else
        % Single FRF numeric: [f H] or [f Re Im]
        if exist('thm_to_complex','file') == 2
            [f, Hc] = thm_to_complex(THM);
        else
            [f, Hc] = thm_to_complex_local(THM);
        end
        f=f(:); Hc=Hc(:); Hmat_all = Hc.'; ich=1; Ny=1;
    end
else
    error('Unsupported THM type. Provide matrix or cell array.');
end

% Build a NUMERIC SISO view for helpers that only accept [f Re Im]
THM_primary = [f, real(Hc), imag(Hc)];

switch iat
    case 1, frf='';           model='multiplicative'; % (transmissibility)
    case 2, frf='acc/force';  model='additive';
    case 3, frf='vel/force';  model='additive';
    case 4, frf='disp/force'; model='additive';
    otherwise, error('Unknown iat=%d. Valid: 1..4', iat);
end

Gdc = pick_DCGain(frf, f, Hc);
if isempty(opt.FitBaseline), fitBaseline = strcmp(frf,'disp/force'); else, fitBaseline = opt.FitBaseline; end

% Prefer project helpers if present
if exist('transfer_from_params','file') == 2
    xfer = @(A,fnn,zz,tau) transfer_from_params(f, A, fnn, zz, 'FRF', frf, 'Model','auto', 'DCGain', Gdc, 'Delay', tau);
else
    xfer = @(A,fnn,zz,tau) transfer_from_params_local(f, A, fnn, zz, 'FRF', frf, 'Model','auto', 'DCGain', Gdc, 'Delay', tau);
end
if exist('complex_rmse','file') == 2, crmse = @(H1,H2) complex_rmse(H1,H2);
else,                                 crmse = @(H1,H2) complex_rmse_local(H1,H2);
end

% Baseline “no modes” (fit model baseline)
H_base_fit = xfer(zeros(0,1), zeros(0,1), zeros(0,1), 0);
err_best = crmse(Hc, H_base_fit);
best = struct('fn',[],'zeta',[],'A',[],'tau',0,'Hhat_fit',H_base_fit,'err',err_best);

% Strict accept-if-better
function ok = improves(newErr, oldErr)
    tol = 1e-9 * max(1, oldErr);
    ok = isfinite(newErr) && (newErr < oldErr - tol);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% -------------------- pre-ID assist (seed candidates) --------------------
try
    kseed = ndof;
    [fn_pre, zeta_pre] = preidentify_modes_complex_local( ...
        f, Hc, kseed, 'MinSepPct', opt.MinSepPct, 'Verbose', opt.Verbose);

    [fn_pre, zeta_pre] = force_exact_ndof_seeds( ...
        f, Hc, fn_pre, zeta_pre, ndof, opt.MinSepPct);

    if ~isempty(fn_pre)
        A_pre  = seed_amplitudes_ls(f, Hc, fn_pre, zeta_pre, frf, Gdc, 0, opt.ForceAreal, 'auto');
        H_pre_fit  = xfer(A_pre, fn_pre, zeta_pre, 0);
        err_pre = crmse(Hc, H_pre_fit);
        if improves(err_pre, best.err)
            best = struct('fn', fn_pre, 'zeta', zeta_pre, 'A', A_pre, ...
                          'tau', 0, 'Hhat_fit', H_pre_fit, 'err', err_pre);
        end
        if opt.Verbose
            fprintf('[pre-ID enforced %d] ', ndof);
            fprintf('%.3f ', fn_pre); fprintf('\n');
        end
    end
catch
    if opt.Verbose
        warning('preidentify seed stage skipped.');
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% -------------------- 1) Initial extraction --------------------
upd(prog,'Initial extraction…');
nd_seed = max(ndof, ndof + (ndof>=3)*2);
quadOpts = struct('Enable',true, 'Weight',2.5, 'UseAbsImag',true, 'RefinePk',true);
jitOpts  = struct('Enable',true, 'Ntrials',900, 'ImproveThresh',0.008, ...
                  'MaxAddsPerMode',2, 'FnPct',0.10, 'MinSepPct',2e-4, 'WindowBW',1.5);
mcOpts   = struct('Enable',true, 'PctWithin',0.10, 'Ntrials',300, 'ImproveThresh',0.015);

[fn0, zeta0, info0, A0] = estimate_modes_halfpower_40(THM_primary, nd_seed, ...
    'FRF', frf, 'Model', model, 'plot', logical(opt.Plot), 'Areal', opt.ForceAreal, ...
    'zeta_guess', 0.010, 'MinRelSep', 1e-4, 'Quad', quadOpts, 'Jitter', jitOpts, 'MC', mcOpts);

if numel(fn0) > ndof, [fn0, zeta0, A0] = pick_top_separated(fn0, zeta0, A0, ndof); end
H0_fit = xfer(A0, fn0, zeta0, 0);

err0 = (isstruct(info0)&&isfield(info0,'rmse_after')&&isfinite(info0.rmse_after))*info0.rmse_after ...
     + ~(isstruct(info0)&&isfield(info0,'rmse_after'))*crmse(Hc,H0_fit);

if improves(err0,best.err)
    best = struct('fn',fn0,'zeta',zeta0,'A',A0,'tau',0,'Hhat_fit',H0_fit,'err',err0);
end
prog = prog + weights.init; upd(prog,'Initial extraction done'); if timedout(), goto_final(); return; end

% -------------------- 2) Missing-mode scan --------------------
added_tbl = []; cand_tbl = []; err_miss = NaN;
if imiss == 1 || numel(best.fn) < ndof
    upd(prog,'Residual-based missing-mode search…');
    Hfit_current = xfer(best.A, best.fn, best.zeta, best.tau);
    [fn_m, zeta_m, A_m, added_tbl, cand_tbl, err_miss] = find_missing_modes_by_ratio_addonly( ...
        THM_primary, best.fn, best.zeta, best.A, ndof, ...
        'Hfit', Hfit_current, 'FRF', frf, 'Model', model, ...
        'DCGain', Gdc, 'Delay', best.tau, 'Areal', opt.ForceAreal, ...
        'MinRelSep', opt.MinSepPct, 'MaxAdds', 2, 'ImproveThr', 0.010, ...
        'ZetaGrid', [0.004 0.007 0.01 0.015 0.02 0.03 0.05 0.08], ...
        'WinBW', 2.2, 'ForceOne', true, 'Verbose', opt.Verbose, 'Plot', logical(opt.Plot));
    if improves(err_miss, best.err)
        Hm_fit = xfer(A_m, fn_m, zeta_m, best.tau);
        best = struct('fn',fn_m,'zeta',zeta_m,'A',A_m,'tau',best.tau,'Hhat_fit',Hm_fit,'err',err_miss);
    end
    prog = prog + weights.missing; upd(prog,'Missing-mode scan done');
end
if timedout(), goto_final(); return; end

% -------------------- 3) Overlap split-refine --------------------
extra = struct('fn',[],'zeta',[]); err_overlap = NaN;
if numel(best.fn) < ndof
    upd(prog,'Overlap split-refine…');
    [extra, err_overlap] = find_missing_overlapped_modes_splitrefine( ...
        THM_primary, ndof, best.fn, best.zeta, ...
        'FRF', frf, 'DCGain', Gdc, 'Delay', best.tau, ...
        'PctWithin', [0.01 0.15], 'WindowBW', 1.75, ...
        'ReplaceAllowed', true, 'MaxPasses', 3, 'Ntrials', 450, ...
        'ImproveLocal', 0.01, 'ImproveGlobal', 0.003, 'Verbose', opt.Verbose);

    if ~isempty(extra.fn)
        fn_try   = [best.fn(:);   extra.fn(:)];
        zeta_try = [best.zeta(:); extra.zeta(:)];
        A_try    = seed_amplitudes_ls(f, Hc, fn_try, zeta_try, frf, Gdc, best.tau, opt.ForceAreal, 'auto');
        H_try_fit    = xfer(A_try, fn_try, zeta_try, best.tau);
        err_try  = crmse(Hc, H_try_fit);
        if isfinite(err_overlap), err_try = err_overlap; end
        if improves(err_try, best.err)
            best = struct('fn',fn_try,'zeta',zeta_try,'A',A_try,'tau',best.tau,'Hhat_fit',H_try_fit,'err',err_try);
        end
    end
    prog = prog + weights.overlap; upd(prog,'Overlap split-refine done');
end
if timedout(), goto_final(); return; end

% -------------------- 4) Trust-region refine --------------------
upd(prog,'Trust-region refinement…');
[fn_r, zeta_r, A_r, H_r_fit, report_r] = refine_modes_trustregion( ...
    f, Hc, best.fn, best.zeta, best.A, ...
    'FRF', frf, 'Model', model, 'DCGain', Gdc, 'Areal', opt.ForceAreal, ...
    'MinSepPct', opt.MinSepPct, 'Verbose', opt.Verbose, ...
    'Delay', best.tau, 'EstimateDelay', opt.EstimateDelay, ...
    'DelayRange', opt.DelayRange, 'DelayGrid', opt.DelayGrid, ...
    'FitBaseline', fitBaseline);
err_ref = report_r.rmse_after;
if improves(err_ref, best.err)
    best = struct('fn',fn_r,'zeta',zeta_r,'A',A_r,'tau',report_r.tau,'Hhat_fit',H_r_fit,'err',err_ref);
end
prog = prog + weights.trust; upd(prog,'Trust-region done'); if timedout(), goto_final(); return; end

% -------------------- 5) JitterBoost --------------------
err_jb = NaN;
if opt.UseJitterBoost
    upd(prog,'JitterBoost refinement…');
    [fn_j, zeta_j, A_j, H_j_fit, report_jb] = refine_modes_jitterboost( ...
        f, Hc, best.fn, best.zeta, best.A, ...
        'FRF', frf, 'Model', model, 'DCGain', Gdc, 'Delay', best.tau, ...
        'Areal', opt.ForceAreal, ...
        'Rounds', opt.JB_Rounds, 'StartsPerRound', opt.JB_Starts, ...
        'FnPct', opt.JB_FnPct, 'ZetaPct', opt.JB_ZetaPct, ...
        'LocalLine', true, 'LineGrid', opt.JB_LineGrid, ...
        'MinSepPct', opt.MinSepPct, 'Verbose', opt.Verbose);
    err_jb = report_jb.rmse_after;
    if improves(err_jb, best.err)
        best = struct('fn',fn_j,'zeta',zeta_j,'A',A_j,'tau',best.tau,'Hhat_fit',H_j_fit,'err',err_jb);
    end
    prog = prog + weights.jitter; upd(prog,'JitterBoost done');
end
if timedout(), goto_final(); return; end

% -------------------- 6) Close-mode resolver --------------------
upd(prog,'Close-mode resolver…');
[fn_res, zeta_res, A_res, H_res_fit, rep_res] = resolve_close_modes_merge_if_needed( ...
    f, Hc, best.fn, best.zeta, best.A, ...
    'FRF', frf, 'Model', model, 'DCGain', Gdc, 'Delay', best.tau, ...
    'Areal', opt.ForceAreal, ...
    'NearPct', 0.01, 'WinFactor', 2.0, 'Criterion', 'BIC', 'DeltaCrit', 2, ...
    'MaxIters', 60, 'Verbose', opt.Verbose);
err_res = crmse(Hc, H_res_fit);
did_merge1 = has_merge_flag(rep_res);
if err_res < best.err - 1e-9*max(1,best.err) || (did_merge1 && err_res <= best.err*(1+1e-12))
    best.fn = fn_res; best.zeta = zeta_res; best.A = A_res; best.Hhat_fit = H_res_fit; best.err = err_res;
end
prog = prog + weights.resolver1; upd(prog,'Close-mode resolver done'); if timedout(), goto_final(); return; end

% -------------------- 7) Levenberg–Marquardt --------------------
upd(prog,'Levenberg–Marquardt refinement…');
[fn_lm, zeta_lm, A_lm, H_lm_fit, rep_lm] = refine_modes_levenberg( ...
    f, Hc, best.fn, best.zeta, best.A, ...
    'FRF', frf, 'Model', model, 'DCGain', Gdc, ...
    'Areal', opt.ForceAreal, ...
    'Delay', best.tau, 'EstimateDelay', opt.EstimateDelay, ...
    'FitBaseline', fitBaseline, ...
    'MaxIters', 120, 'Lambda0', 1e-2, 'MinSepPct', opt.MinSepPct, ...
    'Verbose', opt.Verbose);
err_lm = rep_lm.rmse_after;
if err_lm < best.err - 1e-9*max(1,best.err)
    best = struct('fn',fn_lm,'zeta',zeta_lm,'A',A_lm,'tau',rep_lm.tau,'Hhat_fit',H_lm_fit,'err',err_lm);
end
prog = prog + weights.levenberg; upd(prog,'LM refine done'); if timedout(), goto_final(); return; end

% -------------------- 8) BIC/CMIF 1-vs-2 (multi-FRF aware) --------------------
upd(prog,'Post-hoc 1-vs-2 (BIC/CMIF)…');
if Ny > 1
    A_mat = ls_residues_multi(f, Hmat_all, best.fn, best.zeta, frf, model, Gdc, best.tau, opt.ForceAreal);
    Hmat_in = Hmat_all.';   % Nf x Ny
    useCMIF = true;
else
    A_mat = reshape(best.A, 1, []);
    Hmat_in = Hc;           % Nf x 1
    useCMIF = false;
end
[fn_b, zeta_b, A_mat_b, H_b_fit, rep_b] = resolve_close_modes_bic_cmif( ...
    f, Hmat_in, best.fn, best.zeta, A_mat, ...
    'FRF', frf, 'NearPct', 0.01, 'WinFactor', 2.0, ...
    'DeltaBIC', 2.0, 'UseCMIF', useCMIF, ...
    'Delay', best.tau, 'DCGain', Gdc, 'Verbose', opt.Verbose);

if Ny > 1
    err_b = mean(sqrt(mean(abs(H_b_fit - Hmat_in).^2, 1)));
else
    err_b = crmse(Hc, H_b_fit);
end
did_merge_b = has_merge_flag(rep_b);
if err_b < best.err - 1e-9*max(1,best.err) || (did_merge_b && err_b <= best.err*(1+1e-12))
    best.fn   = fn_b;
    best.zeta = zeta_b;
    if Ny > 1, best.A = A_mat_b(ich,:).'; else, best.A = A_mat_b.'; end
    best.Hhat_fit = xfer(best.A, best.fn, best.zeta, best.tau);
    if Ny > 1, best.err = mean(sqrt(mean(abs(H_b_fit - Hmat_in).^2, 1)));
    else,      best.err = crmse(Hc, best.Hhat_fit);
    end
end
prog = prog + weights.bic; upd(prog,'BIC/CMIF pass done'); if timedout(), goto_final(); return; end

% -------------------- 9) Duplicate collapse --------------------
upd(prog,'Collapsing near-duplicate modes…');
A_mat_c = reshape(best.A,1,[]);
[fn_c, zeta_c, A_mat_c, rep_c] = collapse_duplicate_modes( ...
    best.fn, best.zeta, A_mat_c, ...
    'FreqTolPct', 2.0e-3, 'ZetaTolAbs', 1.0, 'PhaseTolDeg', 200, 'MinAmpRatio', 0.0);
H_c_fit  = xfer(A_mat_c.', fn_c, zeta_c, best.tau);
err_c = crmse(Hc, H_c_fit);
did_collapse = has_merge_flag(rep_c);
if err_c <= best.err*(1+1e-12) || did_collapse
    best.fn = fn_c; best.zeta = zeta_c; best.A = A_mat_c.'; best.Hhat_fit = H_c_fit; best.err = err_c;
end
prog = min(1, prog + weights.collapse); upd(prog,'Collapse done');

% -------------------- Quick LM probe + conditional fast basin-hop --------------------
upd(prog,'Quick LM probe…');
err_before = best.err;
[fn_q, zeta_q, A_q, H_q_fit, rep_q] = refine_modes_levenberg( ...
    f, Hc, best.fn, best.zeta, best.A, ...
    'FRF', frf, 'Model', model, 'DCGain', Gdc, ...
    'Areal', false, 'Delay', best.tau, 'EstimateDelay', false, ...
    'FitBaseline', fitBaseline, 'MaxIters', 40, 'Lambda0', 1e-2, ...
    'MinSepPct', opt.MinSepPct, 'Verbose', opt.Verbose);
err_q = rep_q.rmse_after;
gain_q = max(0, (err_before - err_q) / max(err_before, eps));
if err_q < best.err * (1 - 1e-9)
    best.fn = fn_q; best.zeta = zeta_q; best.A = A_q; best.Hhat_fit = H_q_fit; best.err = err_q;
end

run_basin = gain_q >= opt.SkipBasinIfLMGainBelow;
if run_basin && ibh==1
    upd(prog,'Basin-hop refinement…');
    trials = opt.BasinTrials;
    if opt.BasinAutoScale
        scale = 0.3 + 0.7*min(1, 10*gain_q); % ∈ [0.3, 1.0]
        trials = max(opt.BasinMinTrials, round(scale * opt.BasinTrials));
    end
    decimN = opt.BasinDecimTargetN; if isempty(decimN), decimN = min(1200, numel(f)); end
    patience = opt.BasinEarlyStopPatience; if isempty(patience), patience = max(150, round(0.25*trials)); end

    if exist('refine_modes_basin_mc_fast','file')==2
        [fn_mc, zeta_mc, A_mc, H_mc_fit] = refine_modes_basin_mc_fast( ...
            f, Hc, best.fn, best.zeta, best.A, ...
            'FRF', frf, 'Model', model, 'DCGain', Gdc, 'Delay', best.tau, ...
            'Areal', false, 'MinSepPct', opt.MinSepPct, ...
            'Trials', trials, 'SigmaFnPct', [0.12 0.01], 'SigmaZetaPct', [0.25 0.03], ...
            'Decimate', opt.BasinDecimate, 'DecimTargetN', decimN, 'FocusWinFactor', opt.BasinFocusWinFactor, ...
            'Sobol', opt.BasinSobol, 'UseParallel', opt.BasinUseParallel, ...
            'EarlyStopPatience', patience, 'RefineMaxIters', opt.BasinRefineMaxIters, 'Seed', opt.BasinSeed);
    else
        [fn_mc, zeta_mc, A_mc, H_mc_fit] = refine_modes_basin_mc( ...
            f, Hc, best.fn, best.zeta, best.A, ...
            'FRF', frf, 'Model', model, 'DCGain', Gdc, 'Delay', best.tau, ...
            'Areal', false, 'MinSepPct', opt.MinSepPct, ...
            'Trials', trials, 'SigmaFnPct', [0.12 0.01], 'SigmaZetaPct', [0.25 0.03], ...
            'TemperatureScale', 0.01, 'WeightPeaks', false, 'FitBaseline', false);
    end
    err_mc = crmse(Hc, H_mc_fit);
    if err_mc < best.err * 0.999
        best.fn = fn_mc; best.zeta = zeta_mc; best.A = A_mc; best.Hhat_fit = H_mc_fit; best.err = err_mc;
    end
else
    if opt.Verbose
        fprintf('Skipping basin-hop (quick LM gain %.4g < %.4g or ibh==0).\n', gain_q, opt.SkipBasinIfLMGainBelow);
    end
end

% -------------------- Final LM polish (full) --------------------
upd(prog,'Final Levenberg–Marquardt…');
[fn_lm2, zeta_lm2, A_lm2, H_lm2_fit, rep_lm2] = refine_modes_levenberg( ...
    f, Hc, best.fn, best.zeta, best.A, ...
    'FRF', frf, 'Model', model, 'DCGain', Gdc, ...
    'Areal', false, 'Delay', best.tau, 'EstimateDelay', true, ...
    'FitBaseline', false, 'MaxIters', 200, 'Lambda0', 1e-2, ...
    'MinSepPct', opt.MinSepPct, 'Verbose', opt.Verbose);
if rep_lm2.rmse_after < best.err * 0.999
    best.fn = fn_lm2; best.zeta = zeta_lm2; best.A = A_lm2; best.Hhat_fit = H_lm2_fit; best.err = rep_lm2.rmse_after;
end

% -------------------- Final tidy & info --------------------
goto_final();

function goto_final()

    % ---- sort modes for output only ----
    fn   = best.fn(:);
    zeta = best.zeta(:);
    A    = best.A(:);

    [fn,ord] = sort(fn);
    zeta = zeta(ord);
    A    = A(ord);

    % enforce requested ndof only for output variables
    nkeep = min([ndof numel(fn) numel(zeta) numel(A)]);
    fn   = fn(1:nkeep);
    zeta = zeta(1:nkeep);
    A    = A(1:nkeep);

    % ---- DO NOT rebuild the FRF ----
    % use the model produced by the solver
    Hhat = best.Hhat_fit;

    % error relative to measured FRF
    err_final = crmse(Hc,Hhat);

    upd(1,'Done');

    % plotting
    if opt.Plot
        plot_fit_complex(f,Hc,Hhat,iat,ich);
    end

    % info structure
    info = struct();
    info.frf       = frf;
    info.model     = model;
    info.ndof_out  = numel(fn);
    info.err       = err_final;
    info.primary   = ich;

end

end

% ==================== locals ====================


function Gdc = pick_DCGain(frf, f, Hc)
switch char(frf)
    case '',                         Gdc = 1+0i;     % transmissibility baseline
    case {'acc/force','vel/force'},  Gdc = 0;
    case 'disp/force'
        nDC = max(5, round(0.06*numel(f)));
        Gdc = median(Hc(1:nDC)); if ~isfinite(Gdc), Gdc=0; end
    otherwise, Gdc = 1+0i;
end
end

function update_waitbar_local(hwb, frac, msg)
if isempty(hwb) || ~ishandle(hwb), return; end
try
    if strcmp(get(hwb,'Type'),'figure'), figure(hwb); end
    waitbar(max(0,min(1,frac)), hwb, msg);
    drawnow limitrate nocallbacks
catch
end
end

function try_close_waitbar(hwb)
try, if ~isempty(hwb) && ishandle(hwb), close(hwb); end, catch, end
end

function [fnS, zS, AS] = pick_top_separated(fn, z, A, k)
[~,ord] = sort(abs(A), 'descend'); ord = ord(1:min(k,numel(ord)));
fnS = fn(ord); zS = z(ord); AS = A(ord);
[fnS,ord2] = sort(fnS); zS = zS(ord2); AS = AS(ord2);
end

function [f, Hc] = thm_to_complex_local(THM1)
if size(THM1,2) < 2, error('THM must be [f, H] or [f, Re, Im].'); end
f = double(THM1(:,1));
if size(THM1,2) >= 3 && isreal(THM1(1,2)) && isreal(THM1(1,3))
    Hc = complex(double(THM1(:,2)), double(THM1(:,3)));
else
    Hc = THM1(:,2); if ~isa(Hc,'double'), Hc = double(Hc); end
    if isreal(Hc) && all(isfinite(Hc)) && all(Hc >= 0), Hc = complex(Hc, 0); end
end
f = f(:); Hc = Hc(:);
end

function Hhat = transfer_from_params_local(f, A, fn, zeta, varargin)
p = inputParser;
p.addParameter('DCGain',1+0i);
p.addParameter('Delay',0);
p.addParameter('FRF','');
p.addParameter('Model','auto');
p.parse(varargin{:});
G=p.Results.DCGain; tau=p.Results.Delay;
frfType = lower(char(p.Results.FRF));
modelForm = p.Results.Model;
if strcmpi(modelForm,'auto'), modelForm = iff(strcmp(frfType,''),'multiplicative','additive'); end

f = double(f(:)); A=A(:); fn=fn(:); zeta=zeta(:); w=2*pi*f;
Hcols = zeros(numel(w), numel(fn));
for k = 1:numel(fn)
    Hcols(:,k) = frf_kernel_local(w, 2*pi*fn(k), zeta(k), frfType);
end
if isempty(A), Hbody = zeros(size(f)); else, Hbody = Hcols*A; end

if strcmpi(modelForm,'multiplicative')
    Hhat = G * (1 + Hbody);
else
    Hhat = G + Hbody;
end
if tau~=0, Hhat = exp(-1i*w*tau).*Hhat; end
end

function v = iff(cond,a,b), if cond, v=a; else, v=b; end, end

function Hk = frf_kernel_local(w, wn, zeta, frfType)
den = (wn.^2 - w.^2) + 1i*(2*zeta*wn).*w;
switch lower(frfType)
    case 'disp/force', Hk = 1 ./ den;
    case 'vel/force',  Hk = (1i*w) ./ den;
    case 'acc/force',  Hk = (-w.^2) ./ den;
    otherwise
        % legacy transmissibility kernel used by the FITTING model:
        % H = 1 + sum A_k * (w^2/den).  (kept to preserve your fn/zeta results)
        Hk = (w.^2) ./ den;
end
end

% function rm = complex_rmse_local(Hc, Hhat)
% r = Hc - Hhat; rm = sqrt(mean(abs(r).^2));
% end

function A_mat = ls_residues_multi(f, Hmat, fn, zeta, frf, model, Gdc, tau, Areal)
% Solve residues for each channel with shared poles (fn,zeta)
w = 2*pi*f(:);
Nm = numel(fn); Ny = size(Hmat,1);
Hcols = zeros(numel(w), Nm);
for k=1:Nm, Hcols(:,k) = frf_kernel_local(w, 2*pi*fn(k), zeta(k), lower(char(frf))); end
A_mat = zeros(Ny, Nm);
useMult = strcmpi(model,'multiplicative') || isempty(frf);
for y = 1:Ny
    yvec = Hmat(y,:).';
    if tau~=0, yvec = yvec .* exp(1i*w*tau); end
    if useMult && Gdc~=0
        rhs = yvec./Gdc - 1;
    else
        rhs = yvec - Gdc;
    end
    Ay = Hcols \ rhs;
    if Areal, Ay = real(Ay); end
    A_mat(y,:) = Ay(:).';
end
end

function tf = has_merge_flag(R)
% Safely detect a "merge happened" flag in various report shapes.
tf = false;
try
    if iscell(R) && ~isempty(R), R = R{1}; end
    if ~isstruct(R), return; end
    if isfield(R,'log')
        L = R.log;
        if isstruct(L) && isfield(L,'merge')
            tf = any(L.merge(:));
        elseif isnumeric(L) || islogical(L)
            tf = any(L(:)~=0);
        end
    elseif isfield(R,'merged')
        tf = any(R.merged(:)~=0);
    elseif isfield(R,'idx')
        tf = ~isempty(R.idx);
    end
catch
    tf = false;
end
end

% ==================== helper to enforce exactly ndof seeds ====================
function [fn_out, z_out] = force_exact_ndof_seeds(f, Hc, fn_in, z_in, ndof, minSepPct)
% Ensure exactly ndof seeds by (1) reducing extras with separation-aware
% farthest-point sampling using |A| tie-break; (2) filling deficits from |H-1| peaks.

f = f(:); Hc = Hc(:);
fn = fn_in(:); z  = z_in(:);
ok = isfinite(fn) & isfinite(z) & fn>0 & z>0;
fn = fn(ok); z = z(ok);

if isempty(fn)
    [fn, z] = fill_from_R(f, Hc, [], [], ndof, minSepPct);
end

if numel(fn) > ndof
    try
        Aall = seed_amplitudes_ls(f, Hc, fn, z, '', 1+0i, 0, false, 'auto');
        s = abs(Aall(:));
    catch
        s = ones(numel(fn),1);
    end
    [~, j0] = max(s);
    pick = j0;
    while numel(pick) < ndof
        d = inf(numel(fn),1);
        for j=1:numel(fn)
            d(j) = min(abs(fn(j) - fn(pick)));
        end
        d(pick) = -inf;
        [~, jnext] = max(d + 0.01*s);
        if any(abs(fn(jnext) - fn(pick)) <= minSepPct*max(fn(jnext), eps))
            d(jnext) = -inf; [~, jnext] = max(d + 0.01*s);
        end
        pick(end+1) = jnext; %#ok<AGROW>
    end
    pick = unique(pick,'stable');
    fn = fn(pick); z = z(pick);
end

if numel(fn) < ndof
    need = ndof - numel(fn);
    [fn, z] = fill_from_R(f, Hc, fn, z, numel(fn)+need, minSepPct);
    if numel(fn) > ndof
        [fn, ord] = sort(fn); z = z(ord);
        fn = fn(1:ndof); z = z(1:ndof);
    end
end

[fn_out, ord] = sort(fn(:));
z_out = z(ord);
end

function [fn_new, z_new] = fill_from_R(f, Hc, fn_exist, z_exist, targetK, minSepPct)
% Add seeds from |H-1| peaks with half-height ζ; enforce min separation
R = abs(Hc - 1); R = R./max(R+eps);
N = numel(f);
win = max(3, round(0.0012*N));
R = conv(R(:), ones(win,1)/win, 'same');

y = R(:); n = numel(y);
locs = find([false; (y(2:n-1) > y(1:n-2)) & (y(2:n-1) >= y(3:n)); false]);
pks  = y(locs);
[~,ord] = sort(pks,'descend');
cand = f(locs(ord));
fn_acc = fn_exist(:); z_acc = z_exist(:);

for i = 1:numel(cand)
    if numel(fn_acc) >= targetK, break; end
    fi = cand(i);
    if ~isempty(fn_acc) && any(abs(fn_acc - fi) <= minSepPct*max(fi, eps)), continue; end

    [~, i0] = min(abs(f - fi));
    pk = R(i0); if ~(isfinite(pk) && pk>0), continue; end
    hh = 0.5*pk;
    iL = i0; while iL>1        && R(iL) >= hh, iL = iL-1; end
    iR = i0; while iR<numel(R) && R(iR) >= hh, iR = iR+1; end
    if iL>1 && iR<numel(R)
        fL = f(iL)   + (hh - R(iL))   * (f(iL+1)-f(iL)) / (R(iL+1)-R(iL) + eps);
        fR = f(iR-1) + (hh - R(iR-1)) * (f(iR)-f(iR-1)) / (R(iR)-R(iR-1) + eps);
        bw = max(median(diff(f)), fR - fL);
        zhat = max(5e-3, min(0.12, bw/(2*max(fi, eps))));
    else
        zhat = 0.02;
    end
    fn_acc(end+1,1) = fi; %#ok<AGROW>
    z_acc(end+1,1)  = zhat; %#ok<AGROW>
end

fn_new = fn_acc; z_new = z_acc;
end