function AM = apparent_mass_semidefinite(omega, phi, zeta, omega_n, k1)
% APPARENT_MASS_SEMIDEFINITE  Interface apparent mass M_A(omega) = F_12 / a_1
% for a multi-DOF system driven at DOF 1 (mass m1) through stiffness k1.
% Works for free-free / semidefinite systems (omega_n(1) = 0).
%
% Inputs
%   omega   : excitation frequencies (rad/s)
%   phi     : mode shape matrix (rows = DOFs, cols = modes), mass-normalized
%   zeta    : modal damping ratios, one per mode
%   omega_n : natural frequencies (rad/s), >= 0; 0 allowed for the rigid-body mode
%   k1      : interface stiffness between m1 and m2
%
% Output
%   AM      : complex apparent mass, same shape as omega
%
% Formulation.  With mass-normalized modes and unit force at DOF 1,
%
%   X_j(omega) = sum_n  phi_{j,n} phi_{1,n} / H_n(omega)
%   H_n(omega) = (omega_n^2 - omega^2) + 2j*zeta_n*omega_n*omega
%
% so that, defining F_12 = k1 (X_1 - X_2) and a_1 = -omega^2 X_1,
%
%   M_A(omega) = k1 * sum_n  phi_{1,n}(phi_{1,n}-phi_{2,n}) / H_n
%              --------------------------------------------------------
%                       -omega^2 * sum_n  phi_{1,n}^2 / H_n
%
% Static limit:  M_A(0) = (mass above the k1 interface) = m2 + m3 + ...
% which is obtained from the residue of the rigid-body pole; see DC block.

    omega   = omega(:);
    omega_n = omega_n(:);
    zeta    = zeta(:);
    Nfreq   = numel(omega);
    Nmodes  = numel(omega_n);

    if any(omega_n < 0)
        error('Natural frequencies omega_n must be non-negative.');
    end
    if size(phi,2) ~= Nmodes
        error('phi must have one column per mode (got %d cols, %d modes).', ...
              size(phi,2), Nmodes);
    end

    phi1 = phi(1,:).';
    phi2 = phi(2,:).';

    w_num = phi1 .* (phi1 - phi2);   % rigid-body contribution is zero
    w_den = phi1 .^ 2;

    S_num = zeros(Nfreq, 1);
    S_den = zeros(Nfreq, 1);

    for n = 1:Nmodes
        Hn = (omega_n(n)^2 - omega.^2) + 2i * zeta(n) * omega_n(n) .* omega;
        S_num = S_num + w_num(n) ./ Hn;
        S_den = S_den + w_den(n) ./ Hn;
    end

    AM = (k1 .* S_num) ./ (-omega.^2 .* S_den);

    % --- DC limit -----------------------------------------------------------
    % At omega = 0 both S_num and -omega^2*S_den are singular through the
    % rigid-body pole.  The limit is finite:
    %   F_12(0) = k1 * sum_{elastic} phi_{1,n}(phi_{1,n}-phi_{2,n}) / omega_n^2
    %   a_1(0)  = sum_{rb} phi_{1,n}^2   (rigid-body inertance residue)
    %   M_A(0)  = F_12(0) / a_1(0)
    tol = max(eps, 1e-12*max(abs(omega)));
    zero_idx = abs(omega) <= tol;
    if any(zero_idx)
        elastic = omega_n > 1e-6;
        rb      = ~elastic;
        if any(rb)
            F12_0 = k1 * sum( phi1(elastic) .* (phi1(elastic)-phi2(elastic)) ...
                              ./ omega_n(elastic).^2 );
            a1_0  = sum( phi1(rb).^2 );
            AM(zero_idx) = F12_0 / a1_0;
        else
            AM(zero_idx) = NaN;   % no rigid-body mode
        end
    end
end
