1function [q, v, acc] = plate_modal_free(tau, omega, zeta, q0, v0)
2% PLATE_MODAL_FREE Free decay of one mode after the contact ends.
3%
4% [q, v, acc] = plate_modal_free(tau, omega, zeta, q0, v0) is the
5% displacement, velocity and acceleration of a mass-normalized mode of
6% angular frequency omega and damping ratio zeta at the times tau >= 0
7% after the hammer left, starting from q0 and v0:
8%
9% q(tau) = e^(-alpha tau) (q0 cos(wd tau) + (v0 + alpha q0) / wd sin(wd tau)),
10% alpha = zeta omega, wd = omega sqrt(1 - zeta^2).
11%
12% A rigid-body mode with no restoring force (omega = 0, the two free
13% rotations the finite-element method keeps) simply drifts: q = q0 + v0 tau.
15tau = tau(:)';
16if omega <= 0
17 q = q0 + v0 * tau;
18 v = v0 * ones(size(tau));
19 acc = zeros(size(tau));
20 return
21end
22alpha = zeta * omega;
23wd = omega * sqrt(max(1 - zeta^2, 1e-12));
24decay = exp(-alpha * tau);
25a = q0;
26b = (v0 + alpha * q0) / wd;
27ph = wd * tau;
28cs = cos(ph);
29sn = sin(ph);
30q = decay .* (a * cs + b * sn);
31v = decay .* ((-alpha * a + wd * b) * cs + (-alpha * b - wd * a) * sn);
32acc = -2 * alpha * v - omega^2 * q;
33end