function [zeta, Q] = plate_damping(freqs, support_part, params, P)
% PLATE_DAMPING  Modal damping ratios from three losses.
%
%   [zeta, Q] = plate_damping(freqs, support_part, params, P) combines, for
%   each mode of frequency f,
%     material loss     1 / Q_material                       (params.material.Q)
%     air loading       1 / Q_air,  Q_air = 1100 / sqrt(f / 1 kHz), floored at 50 Hz
%     thread dashpot    c phi_hole^2 / omega, the dashpot c at the hole acting on
%                       the mass-normalized mode's displacement phi_hole there
%   into 1/Q = sum of the three; zeta = min(0.08, 1 / (2 Q)) is the damping
%   ratio and Q the total quality factor.  Rigid-body modes (below 20 Hz, the
%   finite-element method keeps three) get zeta = 0.01 so that momentum can
%   flow into them during the collision without their being clamped.
%
%   These are the constants of the simulators this lab follows: order-of-
%   magnitude figures, not measurements of a particular plate.

freqs = freqs(:)';
support_part = support_part(:)';
omega = 2 * pi * freqs;
q_air = 1100 ./ sqrt(max(freqs / 1000, 0.05));
inv_q = 1 / params.material.Q + 1 ./ q_air;
inv_q = inv_q + params.thread.dashpot * support_part.^2 ./ max(omega, 1);
zeta = min(0.08, 0.5 * inv_q);
zeta(freqs < 20) = 0.01;
Q = 1 ./ inv_q;
end