function dy = plate_fem_rhs(y, phi, nm, omega, twozw, contact_k, hunt_crossley, hammer_mass)
% PLATE_FEM_RHS  Time derivative of the coupled hammer/plate state.
%
%   dy = plate_fem_rhs(y, phi, nm, omega, twozw, contact_k, hunt_crossley, m)
%   for the state y = [q; v; x_h; v_h]:
%
%     q_tt = phi F - 2 zeta omega q_t - omega^2 q     (each mode, mass-normalized)
%     x_h_tt = -F / m                                  (the hammer)
%
%   with F the contact force of plate_fem_force and twozw = 2 zeta omega.

q = y(1:nm);
v = y(nm + 1:2 * nm);
F = plate_fem_force(y, phi, nm, contact_k, hunt_crossley);
qdd = phi' * F - twozw .* v - omega.^2 .* q;
dy = [v; qdd; y(2 * nm + 2); -F / hammer_mass];
end