function AM = apparent_mass(omega, phi, zeta, omega_n, k2)
% APPARENT_MASS  Compute apparent mass FRF at the DOF-1 / DOF-2 interface
%
%   AM = apparent_mass(omega, phi, zeta, omega_n, k2)
%
%   Implements:
%
%              k2 * sum_n [ phi_n1*(phi_n1 - phi_n2) / Hn(omega) ]
%   AM(omega) = ---------------------------------------------------
%              -omega^2 * sum_n [ phi_n1^2 / Hn(omega) ]
%
%   where
%       Hn(omega) = omega_n^2 - omega^2 + 2i*zeta_n*omega_n*omega
%
% INPUTS
%   omega    [Nfreq x 1] or [1 x Nfreq]
%            frequency vector (rad/s)
%
%   phi      [NDOF x Nmodes]
%            mass-normalized mode shape matrix
%            rows = DOFs, columns = modes
%            row 1 = shaker/base drive DOF
%            row 2 = interface DOF
%
%   zeta     [Nmodes x 1] or [1 x Nmodes]
%            modal damping ratios
%
%   omega_n  [Nmodes x 1] or [1 x Nmodes]
%            undamped natural frequencies (rad/s)
%
%   k2       scalar
%            interface stiffness between DOF 1 and DOF 2
%
% OUTPUT
%   AM       [Nfreq x 1]
%            complex apparent mass (kg)
%
% NOTES
%   - phi must be mass-normalized: phi' * M * phi = I
%   - Elastic interface force only: F12 = k2*(X1 - X2)
%   - Interface damping force neglected
%   - At omega = 0, the formula has a singular denominator because of the
%     explicit -omega^2 term. DC/static limit should be handled separately.

    % -----------------------------
    % Input shaping
    % -----------------------------
    omega   = omega(:);      % [Nfreq x 1]
    omega_n = omega_n(:);    % [Nmodes x 1]
    zeta    = zeta(:);       % [Nmodes x 1]

    Nfreq  = length(omega);
    Nmodes = length(omega_n);

    % -----------------------------
    % Input checks
    % -----------------------------
    if nargin < 5
        error('apparent_mass requires 5 inputs: omega, phi, zeta, omega_n, k2');
    end

    if ~isnumeric(omega) || ~isnumeric(phi) || ~isnumeric(zeta) || ~isnumeric(omega_n) || ~isnumeric(k2)
        error('All inputs must be numeric.');
    end

    if ~isscalar(k2)
        error('k2 must be a scalar.');
    end

    if size(phi,1) < 2
        error('phi must have at least 2 rows (DOF 1 and DOF 2).');
    end

    if size(phi,2) ~= Nmodes
        error('Number of columns in phi (%d) must equal number of modes in omega_n (%d).', ...
            size(phi,2), Nmodes);
    end

    if length(zeta) ~= Nmodes
        error('Length of zeta (%d) must equal length of omega_n (%d).', ...
            length(zeta), Nmodes);
    end

    if any(~isfinite(omega))
        error('omega contains non-finite values.');
    end

    if any(~isfinite(omega_n)) || any(omega_n <= 0)
        error('omega_n must contain positive finite values.');
    end

    if any(~isfinite(zeta)) || any(zeta < 0)
        error('zeta must contain finite nonnegative values.');
    end

    if ~isfinite(k2)
        error('k2 must be finite.');
    end

    % -----------------------------
    % Mode shape components at DOFs 1 and 2
    % -----------------------------
    phi1 = phi(1,:).';   % [Nmodes x 1]
    phi2 = phi(2,:).';   % [Nmodes x 1]

    % -----------------------------
    % Modal denominator matrix H
    % H(i,n) = omega_n(n)^2 - omega(i)^2 + 2i*zeta(n)*omega_n(n)*omega(i)
    % size(H) = [Nfreq x Nmodes]
    % -----------------------------
    H = (omega_n.'.^2) ...
      - (omega.^2) ...
      + 2i * omega * (zeta .* omega_n).';

    % -----------------------------
    % Modal weights
    % -----------------------------
    w_num = phi1 .* (phi1 - phi2);   % [Nmodes x 1]
    w_den = phi1 .^ 2;               % [Nmodes x 1]

    % -----------------------------
    % Weighted modal sums over modes
    % For each frequency:
    %   sum_num(i) = sum_n w_num(n)/H(i,n)
    %   sum_den(i) = sum_n w_den(n)/H(i,n)
    % -----------------------------
    sum_num = (1 ./ H) * w_num;      % [Nfreq x 1]
    sum_den = (1 ./ H) * w_den;      % [Nfreq x 1]

    % -----------------------------
    % Apparent mass
    % -----------------------------
    AM = NaN(Nfreq,1);

    tol = 1e-14;
    idx = abs(omega) > tol & abs(sum_den) > tol;

    AM(idx) = (k2 .* sum_num(idx)) ./ (-omega(idx).^2 .* sum_den(idx));

    % Optional DC handling:
    % leave omega=0 as NaN because the dynamic formula is singular there.
    % If desired, user may replace AM(omega==0) with static limit externally.

end