function [fn,zeta,A,logrep] = collapse_duplicate_modes(fn,zeta,A,varargin)
% COLlAPSE_DUPLICATE_MODES  Merge modes that are effectively the same.
% Works with SISO (A is [Nm×1] or [1×Nm]) or Ny×Nm (MIMO residues).
p = inputParser;
p.addParameter('FreqTolPct', 1e-3);     % 0.1% of freq
p.addParameter('ZetaTolAbs', 0.01);     % |Δζ| <= 0.01 absolute
p.addParameter('PhaseTolDeg', 20);      % residue phase diff <= 20°
p.addParameter('MinAmpRatio', 0.05);    % keep if smaller amp not vanishing
p.parse(varargin{:});
opt = p.Results;

fn   = fn(:);
zeta = zeta(:);
if isvector(A), A = reshape(A,1,[]); end  % Ny×Nm
Ny = size(A,1);

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

i = 1; logrep.idx=[]; logrep.merged=false(0,1);
while i < numel(fn)
    f1=fn(i); f2=fn(i+1);
    rel = abs(f2 - f1) / max(f1,eps);
    if rel <= opt.FreqTolPct
        zOK = abs(zeta(i+1)-zeta(i)) <= opt.ZetaTolAbs;
        % residue checks (phase & non-trivial magnitude on both)
        Ai = A(:,i); Aj = A(:,i+1);
        ph = 0; rmag = 0;
        if Ny==1
            ph   = abs(angle(Ai) - angle(Aj));
            rmag = min(abs(Ai),abs(Aj)) / max(abs(Ai),abs(Aj)+eps);
        else
            % MAC as a robust phase/collinearity proxy
            mac = (abs(Ai' * Aj)^2) / (max(Ai'*Ai,eps)*max(Aj'*Aj,eps));
            ph  = acos( min(max(sqrt(mac),0),1) ); % 0..pi, small means aligned
            rmag= min(norm(Ai),norm(Aj)) / max(norm(Ai),norm(Aj)+eps);
        end
        phOK = ph <= deg2rad(opt.PhaseTolDeg);
        rOK  = rmag >= opt.MinAmpRatio;

        if zOK && phOK && rOK
            % merge: amplitude add, damping weighted by |A|
            w1 = norm(A(:,i)); w2 = norm(A(:,i+1)); w = w1+w2+eps;
            zeta(i) = (w1*zeta(i) + w2*zeta(i+1))/w;
            % keep frequency as weighted average too (keeps 1218)
            fn(i)   = (w1*f1 + w2*f2)/w;
            A(:,i)  = A(:,i) + A(:,i+1);
            % delete neighbor
            fn(i+1)=[]; zeta(i+1)=[]; A(:,i+1)=[];
            logrep.idx(end+1,:)=[i i+1]; logrep.merged(end+1,1)=true;
            % re-check new neighbor without i++
            continue
        end
    end
    i = i + 1;
end
end
