function [q, v, acc] = plate_modal_free(tau, omega, zeta, q0, v0)
% PLATE_MODAL_FREE  Free decay of one mode after the contact ends.
%
%   [q, v, acc] = plate_modal_free(tau, omega, zeta, q0, v0) is the
%   displacement, velocity and acceleration of a mass-normalized mode of
%   angular frequency omega and damping ratio zeta at the times tau >= 0
%   after the hammer left, starting from q0 and v0:
%
%     q(tau) = e^(-alpha tau) (q0 cos(wd tau) + (v0 + alpha q0) / wd sin(wd tau)),
%     alpha = zeta omega,  wd = omega sqrt(1 - zeta^2).
%
%   A rigid-body mode with no restoring force (omega = 0, the two free
%   rotations the finite-element method keeps) simply drifts: q = q0 + v0 tau.

tau = tau(:)';
if omega <= 0
    q = q0 + v0 * tau;
    v = v0 * ones(size(tau));
    acc = zeros(size(tau));
    return
end
alpha = zeta * omega;
wd = omega * sqrt(max(1 - zeta^2, 1e-12));
decay = exp(-alpha * tau);
a = q0;
b = (v0 + alpha * q0) / wd;
ph = wd * tau;
cs = cos(ph);
sn = sin(ph);
q = decay .* (a * cs + b * sn);
v = decay .* ((-alpha * a + wd * b) * cs + (-alpha * b - wd * a) * sn);
acc = -2 * alpha * v - omega^2 * q;
end