function [fn_out, zeta_out, A_out, added_tbl, cand_tbl, err_after] = find_missing_modes_by_ratio_addonly( ...
    THM, fn_in, zeta_in, A_in, ndof, varargin)

p = inputParser;
p.addParameter('Hfit', [], @(x)isnumeric(x));
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('MinRelSep',4e-4,@(x)isnumeric(x)&&isscalar(x));
p.addParameter('MaxAdds',2,@(x)isnumeric(x)&&isscalar(x));
p.addParameter('ImproveThr',0.01,@(x)isnumeric(x)&&isscalar(x));
p.addParameter('ZetaGrid',[0.004 0.007 0.01 0.015 0.02 0.03 0.05 0.08],@isnumeric);
p.addParameter('WinBW',2.2,@(x)isnumeric(x)&&isscalar(x));
p.addParameter('ForceOne',true,@(b)islogical(b)&&isscalar(b));
p.addParameter('Verbose',true,@(b)islogical(b)&&isscalar(b));
p.addParameter('Plot',false,@(b)islogical(b)&&isscalar(b));
p.parse(varargin{:});
opt = p.Results;

[f, Hc] = thm_to_complex(THM);
f = f(:); Hc = Hc(:);
frf = char(opt.FRF);
model = iff(strcmpi(opt.Model,'auto'), iff(isempty(frf),'multiplicative','additive'), lower(char(opt.Model)));

% Baseline current fit (if not provided)
if isempty(opt.Hfit)
    Hfit = transfer_from_params_local(f, A_in, fn_in, zeta_in, ...
        'FRF',frf,'Model',model,'DCGain',opt.DCGain,'Delay',opt.Delay);
else
    Hfit = opt.Hfit(:);
end

% Residual magnitude
if strcmp(model,'multiplicative')
    R = (Hc./max(opt.DCGain,eps)) - (Hfit./max(opt.DCGain,eps));
else
    R = Hc - Hfit;
end
rmag = abs(R);

% Candidate peaks
sm = movmedian(rmag, max(5, round(numel(f)/400)));
if exist('sgolayfilt','file')==2 && numel(sm)>=11, sm = sgolayfilt(sm,3,11); end
[pk,loc] = findpeaks(sm, 'NPeaks', max(8, 2*opt.MaxAdds), 'SortStr','descend', ...
                          'MinPeakDistance', max(5, round(numel(f)/400)));
cand_tbl = table(f(loc), pk, 'VariableNames', {'f_Hz','res_peak'});

fn_added=[]; zeta_added=[]; A_added=[];
for ci = 1:numel(loc)
    if numel(fn_in) + numel(fn_added) >= ndof, break; end
    fpk = f(loc(ci));
    best_impr = -Inf; best_z = NaN; best_A = 0;

    for z = opt.ZetaGrid(:).'
        fn_try = [fn_in; fn_added; fpk];
        z_try  = [zeta_in; zeta_added; z];
        A_try  = [A_in; A_added; 0];

        A_last = seed_amplitudes_ls(f, Hc, fpk, z, frf, opt.DCGain, opt.Delay, opt.Areal, model);
        A_try(end) = A_last;

        Htry = transfer_from_params_local(f, A_try, fn_try, z_try, ...
            'FRF',frf,'Model',model,'DCGain',opt.DCGain,'Delay',opt.Delay);

        err0 = complex_rmse(Hc, Hfit);
        err1 = complex_rmse(Hc, Htry);
        impr = max(0, (err0 - err1)/max(err0,eps));

        if impr > best_impr
            best_impr = impr; best_z = z; best_A = A_last;
        end
    end

    if best_impr >= opt.ImproveThr || (opt.ForceOne && isempty(fn_added))
        if any( abs(fpk - [fn_in; fn_added]) < opt.MinRelSep*max(fpk,1) )
            continue;
        end
        fn_added(end+1,1)=fpk; zeta_added(end+1,1)=best_z; A_added(end+1,1)=best_A; %#ok<AGROW>
    end
    if numel(fn_added) >= opt.MaxAdds, break; end
end

% Merge & cap
fn_out   = [fn_in(:);   fn_added(:)];
zeta_out = [zeta_in(:); zeta_added(:)];
A_out    = [A_in(:);    A_added(:)];
if numel(fn_out)>ndof, fn_out=fn_out(1:ndof); zeta_out=zeta_out(1:ndof); A_out=A_out(1:ndof); end

% Compute error of the proposed set
H_out    = transfer_from_params_local(f, A_out, fn_out, zeta_out, ...
               'FRF',frf,'Model',model,'DCGain',opt.DCGain,'Delay',opt.Delay);
err_after = complex_rmse(Hc, H_out);

added_tbl = table(fn_added(:), zeta_added(:), A_added(:), ...
    'VariableNames', {'fn_Hz','zeta','A'});

% Optional plot
if opt.Plot
    figure('Color','w','Name','Missing-mode scan (residual)');
    subplot(2,1,1); plot(f, rmag, '-', f(loc), rmag(loc),'ro'); grid on
    xlabel('Hz'); ylabel('|Residual|'); title('Residual magnitude & candidate peaks')
    subplot(2,1,2); plot(f, real(Hc), f, real(Hfit), '--'); grid on
    xlabel('Hz'); ylabel('Real'); legend('meas','fit')
end
end

% ---- helpers ----
function [f,Hc] = thm_to_complex(THM)
if size(THM,2)>=3 && isreal(THM(1,2)) && isreal(THM(1,3))
    f = THM(:,1); Hc = complex(THM(:,2), THM(:,3));
else
    f = THM(:,1); Hc = THM(:,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
function y  = iff(c,a,b), if c, y=a; else, y=b; end, end
