function [extra, err_after] = find_missing_overlapped_modes_splitrefine(THM, ndof, fn0, zeta0, varargin)
% Propose extra modes by splitting existing seeds; also returns merged-set RMSE.

p = inputParser;
p.addParameter('FRF','',@(s)ischar(s)||isstring(s));
p.addParameter('DCGain',1+0i,@(x)isnumeric(x));
p.addParameter('Delay',0,@(x)isnumeric(x)&&isscalar(x));
p.addParameter('PctWithin',[0.01 0.15],@isnumeric);
p.addParameter('WindowBW',1.75,@(x)isnumeric(x));
p.addParameter('ReplaceAllowed',true,@(b)islogical(b)&&isscalar(b)); %#ok<NASGU>
p.addParameter('MaxPasses',3,@(x)isnumeric(x)&&isscalar(x));
p.addParameter('Ntrials',300,@(x)isnumeric(x)&&isscalar(x));
p.addParameter('ImproveLocal',0.01,@(x)isnumeric(x));
p.addParameter('ImproveGlobal',0.003,@(x)isnumeric(x)); %#ok<NASGU>
p.addParameter('Verbose',true,@(b)islogical(b)&&isscalar(b));
p.parse(varargin{:});
o = p.Results;

[f,Hc] = thm_to_complex(THM); f=f(:); Hc=Hc(:);
frf = char(o.FRF); Gdc=o.DCGain; tau=o.Delay;

need = max(0, ndof - numel(fn0));
extra.fn = []; extra.zeta = []; err_after = NaN;
if need==0 || isempty(fn0)
    return;
end

rng(1);
accepted = 0;
for pass = 1:o.MaxPasses
    for k = 1:numel(fn0)
        fk=fn0(k); zk=zeta0(k);
        if ~isfinite(fk)||~isfinite(zk)||fk<=0||zk<=0, continue; end
        best = struct('fn',NaN,'zeta',NaN,'impr',-Inf);

        win = (f>=fk-o.WindowBW & f<=fk+o.WindowBW);

        for t=1:o.Ntrials
            r  = (o.PctWithin(1) + (o.PctWithin(2)-o.PctWithin(1))*rand())*fk;
            f1 = max(f(2), fk - r); f2 = min(f(end-1), fk + r);
            if f2 <= f1*(1+3e-4), f2 = f1*(1+3e-4)*1.001; end
            z1 = min(max(zk*(1+0.2*randn), 1e-4), 0.30);
            z2 = min(max(zk*(1+0.2*randn), 1e-4), 0.30);

            A1 = seed_amplitudes_ls(f(win), Hc(win), fk, zk, frf, Gdc, tau, true, 'auto');
            H1 = transfer_from_params_local(f, [A1], [fk], [zk], 'FRF',frf,'Model','auto','DCGain',Gdc,'Delay',tau);
            err1 = complex_rmse(Hc, H1);

            A2a = seed_amplitudes_ls(f(win), Hc(win), f1, z1, frf, Gdc, tau, true, 'auto');
            A2b = seed_amplitudes_ls(f(win), Hc(win), f2, z2, frf, Gdc, tau, true, 'auto');
            H2  = transfer_from_params_local(f, [A2a;A2b], [f1;f2], [z1;z2], 'FRF',frf,'Model','auto','DCGain',Gdc,'Delay',tau);
            err2 = complex_rmse(Hc, H2);

            impr = max(0,(err1 - err2)/max(err1,eps));
            if impr > best.impr
                if abs(f2-fk) >= abs(f1-fk), best.fn=f2; best.zeta=z2;
                else,                         best.fn=f1; best.zeta=z1;
                end
                best.impr = impr;
            end
        end

        if best.impr >= o.ImproveLocal
            extra.fn(end+1,1)   = best.fn;   %#ok<AGROW>
            extra.zeta(end+1,1) = best.zeta; %#ok<AGROW>
            accepted = accepted + 1;
            if accepted >= need, break; end
        end
    end
    if accepted >= need, break; end
end

% Compute RMSE for merged set suggested by this function
if ~isempty(extra.fn)
    fn_all   = [fn0(:); extra.fn(:)];
    zeta_all = [zeta0(:); extra.zeta(:)];
    A_all    = seed_amplitudes_ls(f, Hc, fn_all, zeta_all, frf, Gdc, tau, true, 'auto');
    H_all    = transfer_from_params_local(f, A_all, fn_all, zeta_all, 'FRF',frf,'Model','auto','DCGain',Gdc,'Delay',tau);
    err_after = complex_rmse(Hc, H_all);
end
end

% helpers
function [f,Hc] = thm_to_complex(T)
if size(T,2)>=3 && isreal(T(1,2)) && isreal(T(1,3))
    f=T(:,1); Hc=complex(T(:,2),T(:,3));
else
    f=T(:,1); Hc=T(:,2); if isreal(Hc), Hc=complex(Hc,0); end
end
f=f(:); Hc=Hc(:);
end
function rm=complex_rmse(a,b), r=a(:)-b(:); rm=sqrt(mean(abs(r).^2)); end
