function [fn_out, zeta_out, A_out, Hhat_out, report] = ...
    resolve_close_modes_merge_if_needed(f, Hc, fn, zeta, A, varargin)
% Resolve pairs of modes closer than NearPct (default 1%) by
% testing TWO-mode vs ONE-merged-mode using global RMSE + AIC/BIC.
%
% Usage:
%   [fn,zeta,A,Hhat,rep] = resolve_close_modes_merge_if_needed( ...
%       f,Hc,fn,zeta,A,'FRF',frf,'Model','auto','DCGain',Gdc,'Delay',tau, ...
%       'Areal',true,'NearPct',0.01,'WinFactor',2.0,'Criterion','BIC', ...
%       'DeltaCrit',2,'MaxIters',60,'Verbose',true)
%
% Output 'report' lists each resolved pair with metrics.

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('NearPct',0.01,@(x)isnumeric(x)&&isscalar(x)&&x>0);    % 1% default
p.addParameter('WinFactor',2.0,@(x)isnumeric(x)&&isscalar(x)&&x>0);   % window ~ factor×sum half-power BW
p.addParameter('Criterion','BIC',@(s)ischar(s)||isstring(s));         % 'BIC'|'AIC'
p.addParameter('DeltaCrit',2,@(x)isnumeric(x)&&isscalar(x));          % prefer simpler if Δ≥2
p.addParameter('MaxIters',60,@(x)isnumeric(x)&&isscalar(x)&&x>=1);    % local LM its
p.addParameter('Verbose',true,@(b)islogical(b)&&isscalar(b));
p.parse(varargin{:});
o=p.Results;

f=f(:); Hc=Hc(:); w=2*pi*f;
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;

% convenience
composeFull = @(fn_,z_,A_) transfer_from_params_local(f, A_, fn_, z_, ...
                        'FRF',frf,'Model',model,'DCGain',Gdc,'Delay',tau);
rmse = @(H) sqrt(mean(abs(Hc(:)-H(:)).^2));

% find close pairs (sorted by closeness)
[fn, ord] = sort(fn(:)); zeta = zeta(:); zeta=zeta(ord); A = A(:); A=A(ord);
pairs = [];
for k=1:numel(fn)-1
    if abs(fn(k+1)-fn(k))/max(1e-9,0.5*(fn(k+1)+fn(k))) <= o.NearPct
        pairs(end+1,:) = [k k+1]; %#ok<AGROW>
    end
end
fn_out=fn; zeta_out=zeta; A_out=A;
Hhat_curr = composeFull(fn_out,zeta_out,A_out);
err_curr  = rmse(Hhat_curr);

logs = [];

for ip = size(pairs,1):-1:1   % work from high indices downward (safe deletion)
    i = pairs(ip,1); j = pairs(ip,2);

    % Build window around the pair (for stability of local fits)
    BW_i = 2*zeta_out(i)*fn_out(i);  BW_j = 2*zeta_out(j)*fn_out(j);
    mid  = 0.5*(fn_out(i)+fn_out(j));
    pad  = o.WinFactor * 0.5*(BW_i + BW_j);
    win  = (f >= max(f(2),mid-pad) & f <= min(f(end-1),mid+pad));
    if nnz(win) < 12, win = true(size(f)); end  % fallback to global

    % Hypothesis T (two-mode): locally re-polish the pair only
    [fn2, z2, A2] = local_refit_pair(fn_out, zeta_out, A_out, [i j], ...
                         f, Hc, frf, model, Gdc, tau, Areal, o.MaxIters, win);
    H2 = composeFull(fn2,z2,A2);
    err2 = rmse(H2);

    % Hypothesis S (single merged): replace (i,j) with one mode then local fit
    [fn1, z1, A1] = local_merge_and_refit(fn_out, zeta_out, A_out, [i j], ...
                         f, Hc, frf, model, Gdc, tau, Areal, o.MaxIters, win);
    H1 = composeFull(fn1,z1,A1);
    err1 = rmse(H1);

    % Information criterion on GLOBAL residuals (real+imag stacked)
    n   = 2*numel(f);           % stacked residual length
    SSE1 = sum(abs(Hc(:)-H1(:)).^2);
    SSE2 = sum(abs(Hc(:)-H2(:)).^2);
    k1  = 3; k2 = 6;            % params for the pair (fn,zeta,A)*{1 or 2}
    switch upper(o.Criterion)
        case 'AIC'
            AIC1 = n*log(max(SSE1/n, eps)) + 2*k1;
            AIC2 = n*log(max(SSE2/n, eps)) + 2*k2;
            dC = AIC1 - AIC2;   % negative favors single
        otherwise % 'BIC'
            BIC1 = n*log(max(SSE1/n, eps)) + k1*log(n);
            BIC2 = n*log(max(SSE2/n, eps)) + k2*log(n);
            dC = BIC1 - BIC2;   % negative favors single
    end

    % Decision: merge if single is not worse in RMSE and wins by ΔCrit in IC
    doMerge = (err1 <= err2*(1+1e-12)) && (dC <= -abs(o.DeltaCrit));

    logs = [logs; table(i,j,mid,err2,err1,dC,doMerge, ...
        'VariableNames',{'idx_i','idx_j','mid_Hz','rmse_two','rmse_one','dCrit','merge'})]; %#ok<AGROW>

    if doMerge
        fn_out = fn1; zeta_out = z1; A_out = A1;
    else
        fn_out = fn2; zeta_out = z2; A_out = A2;
    end

    % Recompute pairs map indices after structural change:
    Hhat_curr = composeFull(fn_out,zeta_out,A_out);
    err_curr  = rmse(Hhat_curr);

    % Rebuild near list in the neighborhood only (simple approach: full rebuild)
    [fn_out, ord2] = sort(fn_out(:)); zeta_out=zeta_out(ord2); A_out=A_out(ord2);
end

Hhat_out = composeFull(fn_out,zeta_out,A_out);
report = struct('log',logs,'rmse_after',rmse(Hhat_out),'criterion',upper(o.Criterion));
end

% ====== local helpers (file-local) =====================================

function [fn2,z2,A2] = local_refit_pair(fn,zeta,A,pairIdx, ...
        f,Hc,frf,model,Gdc,tau,Areal,maxIters,win)
i=pairIdx(1); j=pairIdx(2);
keep = true(size(fn)); keep([i j]) = false;

% fixed part (others)
composeFull = @(fn_,z_,A_) transfer_from_params_local(f,A_,fn_,z_, ...
    'FRF',frf,'Model',model,'DCGain',Gdc,'Delay',tau);

% start from current values; small LM on just the pair
fn2=fn; z2=zeta; A2=A;

% small LM loop on (fn_i,z_i,A_i, fn_j,z_j,A_j)
params = [fn2([i j]); z2([i j]); A2([i j])];  % 6×1
for it=1:maxIters
    [r,J] = residual_J_pair(params, fn2,z2,A2,[i j], f,Hc,frf,model,Gdc,tau,Areal,win);
    Hn = J.'*J; g = J.'*r;
    step = - (Hn + 1e-2*diag(diag(Hn))) \ g;
    params = params + step;

    % project: bounds & small separation
    params(1:2) = min(max(params(1:2), f(2)), f(end-1));
    params(3:4) = min(max(params(3:4), 1e-4), 0.30);
    % enforce order & min separation
    [fpair,ord] = sort(params(1:2)); zpair=params(3:4); zpair=zpair(ord); Apair=params(5:6); Apair=Apair(ord);
    if fpair(2) <= fpair(1)*(1+3e-4), fpair(2)=fpair(1)*(1+3e-4); end
    params = [fpair; zpair; Apair];
end
fn2([i j]) = params(1:2); z2([i j]) = params(3:4); A2([i j]) = params(5:6);

% tidy
[fn2,ord2]=sort(fn2(:)); z2=z2(ord2); A2=A2(ord2);
end

function [fn1,z1,A1] = local_merge_and_refit(fn,zeta,A,pairIdx, ...
        f,Hc,frf,model,Gdc,tau,Areal,maxIters,win)
i=pairIdx(1); j=pairIdx(2);

% initial merged seeds
wA = abs([A(i) A(j)]) + eps;
f0 = (wA(1)*fn(i) + wA(2)*fn(j))/sum(wA);
z0 = (wA(1)*zeta(i) + wA(2)*zeta(j))/sum(wA);
A0 = (A(i) + A(j));

% replace pair by one entry
mask = true(size(fn)); mask([i j]) = false;
fn1 = fn(mask); z1 = zeta(mask); A1 = A(mask);
fn1(end+1,1)=f0; z1(end+1,1)=z0; A1(end+1,1)=A0;

% LM on the single mode only
idx = numel(fn1);
params = [fn1(idx); z1(idx); A1(idx)];  % 3×1
for it=1:maxIters
    [r,J] = residual_J_single(params, fn1,z1,A1, idx, f,Hc,frf,model,Gdc,tau,Areal,win);
    Hn = J.'*J; g = J.'*r;
    step = - (Hn + 1e-2*diag(diag(Hn))) \ g;
    params = params + step;

    % project
    params(1) = min(max(params(1), f(2)), f(end-1));
    params(2) = min(max(params(2), 1e-4), 0.30);
end
fn1(idx)=params(1); z1(idx)=params(2); A1(idx)=params(3);

% tidy
[fn1,ord2]=sort(fn1(:)); z1=z1(ord2); A1=A1(ord2);
end

function [r,J] = residual_J_pair(params, fn,zeta,A,pairIdx, ...
        f,Hc,frf,model,Gdc,tau,Areal,win)
i=pairIdx(1); j=pairIdx(2);
fi=params(1); fj=params(2); zi=params(3); zj=params(4); Ai=params(5); Aj=params(6);

% others fixed
mask = true(size(fn)); mask([i j])=false;
H_fixed = transfer_from_params_local(f, A(mask), fn(mask), zeta(mask), ...
            'FRF',frf,'Model',model,'DCGain',Gdc,'Delay',tau);

% var part (pair)
w = 2*pi*f;
[Hi, dHi_dwn, dHi_dz] = kernel_and_derivs(w, 2*pi*fi, zi, frf);
[Hj, dHj_dwn, dHj_dz] = kernel_and_derivs(w, 2*pi*fj, zj, frf);

if strcmp(model,'multiplicative')
    H = Gdc * (1 + (H_fixed - Gdc)/Gdc + Ai*Hi + Aj*Hj);
else
    H = H_fixed + Ai*Hi + Aj*Hj;
end
if tau~=0, H = H .* exp(-1i*w*tau); end
rC = Hc - H;

% restrict to window for Jacobian weighting (residual still global)
Wloc = true(size(f)); Wloc(win) = true; %#ok<NASGU> % keep option open
W = 1 + (abs(Hc)./max(abs(Hc)+eps)).^2;
Ws = sqrt(W(:));

% ∂H/∂θ (global), then stack −∂H/∂θ
dH_dfi = (strcmp(model,'multiplicative')*Gdc + strcmp(model,'additive')*1) * (Ai*(dHi_dwn*2*pi));
dH_dfj = (strcmp(model,'multiplicative')*Gdc + strcmp(model,'additive')*1) * (Aj*(dHj_dwn*2*pi));
dH_dzi = (strcmp(model,'multiplicative')*Gdc + strcmp(model,'additive')*1) * (Ai*dHi_dz);
dH_dzj = (strcmp(model,'multiplicative')*Gdc + strcmp(model,'additive')*1) * (Aj*dHj_dz);
dH_dAi = (strcmp(model,'multiplicative')*Gdc + strcmp(model,'additive')*1) * (Hi);
dH_dAj = (strcmp(model,'multiplicative')*Gdc + strcmp(model,'additive')*1) * (Hj);
if tau~=0
    E = exp(-1i*w*tau);
    dH_dfi = dH_dfi.*E; dH_dfj=dH_dfj.*E; dH_dzi=dH_dzi.*E; dH_dzj=dH_dzj.*E;
    dH_dAi = dH_dAi.*E; dH_dAj=dH_dAj.*E;
end

r = [real(rC).*Ws; imag(rC).*Ws];
J = [ ...
    -real(dH_dfi).*Ws, -real(dH_dfj).*Ws, -real(dH_dzi).*Ws, -real(dH_dzj).*Ws, -real(dH_dAi).*Ws, -real(dH_dAj).*Ws; ...
    -imag(dH_dfi).*Ws, -imag(dH_dfj).*Ws, -imag(dH_dzi).*Ws, -imag(dH_dzj).*Ws, -imag(dH_dAi).*Ws, -imag(dH_dAj).*Ws  ...
    ];
end

function [r,J] = residual_J_single(params, fn,zeta,A,idx, ...
        f,Hc,frf,model,Gdc,tau,Areal,win)
fk=params(1); zk=params(2); Ak=params(3);

mask = true(size(fn)); mask(idx)=false;
H_fixed = transfer_from_params_local(f, A(mask), fn(mask), zeta(mask), ...
            'FRF',frf,'Model',model,'DCGain',Gdc,'Delay',tau);

w = 2*pi*f;
[Hk, dHk_dwn, dHk_dz] = kernel_and_derivs(w, 2*pi*fk, zk, frf);
if strcmp(model,'multiplicative')
    H = Gdc * (1 + (H_fixed - Gdc)/Gdc + Ak*Hk);
else
    H = H_fixed + Ak*Hk;
end
if tau~=0, H = H .* exp(-1i*w*tau); end
rC = Hc - H;

W = 1 + (abs(Hc)./max(abs(Hc)+eps)).^2; Ws = sqrt(W(:));
dH_dfk = (strcmp(model,'multiplicative')*Gdc + strcmp(model,'additive')*1) * (Ak*(dHk_dwn*2*pi));
dH_dzk = (strcmp(model,'multiplicative')*Gdc + strcmp(model,'additive')*1) * (Ak*dHk_dz);
dH_dAk = (strcmp(model,'multiplicative')*Gdc + strcmp(model,'additive')*1) * (Hk);
if tau~=0
    E = exp(-1i*w*tau);
    dH_dfk=dH_dfk.*E; dH_dzk=dH_dzk.*E; dH_dAk=dH_dAk.*E;
end

r = [real(rC).*Ws; imag(rC).*Ws];
J = [ -real(dH_dfk).*Ws, -real(dH_dzk).*Ws, -real(dH_dAk).*Ws; ...
      -imag(dH_dfk).*Ws, -imag(dH_dzk).*Ws, -imag(dH_dAk).*Ws ];
end

function [Hk, dHk_dwn, dHk_dz] = kernel_and_derivs(w, wn, zeta, frfType)
den = (wn.^2 - w.^2) + 1i*(2*zeta*wn).*w;
switch lower(char(frfType))
    case 'disp/force', N = 1;
    case 'vel/force',  N = 1i*w;
    case 'acc/force',  N = -(w.^2);
    otherwise,         N =  (w.^2);   % transmissibility
end
Hk = N ./ den;
dden_dwn = 2*wn + 1i*(2*zeta).*w;
dden_dz  = 1i*(2*wn).*w;
dHk_dwn  = -(N .* dden_dwn) ./ (den.^2);
dHk_dz   = -(N .* dden_dz)  ./ (den.^2);
end

function y=iff(c,a,b), if c, y=a; else, y=b; end, end
