function H = plate_hertz(params)
% PLATE_HERTZ  Hertz contact constants of the hammer tip on the plate.
%
%   H = plate_hertz(params) returns
%     estar           effective modulus, 1 / ((1 - nu_h^2) / E_h + (1 - nu^2) / E)  [Pa]
%     contact_k       Hertz stiffness 4/3 E* sqrt(R): force = k delta^(3/2)          [N/m^1.5]
%     duration        contact time of a sphere on a rigid half-space,
%                     2.94 (m^2 / (R E*^2 v))^(1/5)                                 [s]
%     impulse         momentum transferred, (1 + e) m v                              [N s]
%     shape_integral  mean of the sin^(3/2) pulse shape over its duration,
%                     Gamma(5/4) / (sqrt(pi) Gamma(7/4))
%     peak_force      impulse / (duration * shape_integral)                           [N]
%     hunt_crossley   the damping factor of the coupled model, 3 (1 - e) / (2 v)      [s/m]
%
%   The Rayleigh-Ritz method applies the pulse of this duration and impulse
%   directly; the finite-element method integrates the Hunt-Crossley force
%   k delta^(3/2) (1 + hunt_crossley delta_dot) against the moving plate and
%   uses the duration only to size its time step.

hm = params.hammer;
mat = params.material;
inv_estar = (1 - hm.nu^2) / hm.E + (1 - mat.nu^2) / mat.E;
H.estar = 1 / inv_estar;
H.contact_k = 4 / 3 * H.estar * sqrt(hm.tip_radius);
speed = max(hm.speed, 1e-9);
H.duration = 2.94 * (hm.mass^2 / (hm.tip_radius * H.estar^2 * speed))^0.2;
H.impulse = (1 + hm.restitution) * hm.mass * hm.speed;
H.shape_integral = gamma(1.25) / (sqrt(pi) * gamma(1.75));
H.peak_force = H.impulse / (H.duration * H.shape_integral);
H.hunt_crossley = 1.5 * (1 - hm.restitution) / speed;
end