function C = plate_rr_contact(params, P, freqs, zeta, phi_hit, t)
% PLATE_RR_CONTACT  The prescribed Hertz pulse and the modal states it leaves.
%
%   C = plate_rr_contact(params, P, freqs, zeta, phi_hit, t) applies the
%   force history
%
%     F(s) = I sin(pi s / T)^(3/2) / (T Gamma(5/4) / (sqrt(pi) Gamma(7/4))),  0 <= s <= T,
%
%   of duration T and impulse I from plate_hertz (the rounded pulse of a
%   Hertz contact on a rigid target) to every mode, whose response is the
%   Duhamel integral with the damped impulse response
%
%     h(u) = e^(-alpha u) sin(wd u) / wd,   alpha = zeta omega,  wd = omega sqrt(1 - zeta^2).
%
%   The state at separation (C.q0, C.v0) is a 64-point Gauss-Legendre
%   quadrature of the integral over [0, T]; the states at the audio samples
%   inside the contact (C.qs, C.vs at C.ts) and on a 600-point plotting grid
%   (C.qg, C.vg at C.tg) are 32-point quadratures over [0, s].  The plate
%   does not move the hammer back in this model: the force is what it would
%   be on a rigid target, which is why the peak force and the impulse differ
%   from the coupled finite-element contact (see plate_fem_contact).

H = P.hertz;
T = H.duration;
I = H.impulse;
nrm = T * H.shape_integral;
freqs = freqs(:)';
omega = 2 * pi * freqs;
zeta = zeta(:)';
alpha = zeta .* omega;
wd = omega .* sqrt(max(1 - zeta.^2, 1e-12));
phi = phi_hit(:)';
nm = numel(freqs);

% terminal state
[z, wz] = plate_gauss_legendre(64);
s = 0.5 * T * (z + 1);
ws = 0.5 * T * wz;
force = I * sin(pi * s / T).^1.5 / nrm;
u = T - s;
E = exp(-alpha' * u);
SW = sin(wd' * u);
CW = cos(wd' * u);
h = E .* SW ./ wd';
hd = E .* (CW - (alpha' ./ wd') .* SW);
q0 = phi .* (h * (ws .* force)')';
v0 = phi .* (hd * (ws .* force)')';

% states inside the contact: at the audio samples and on the plotting grid
ts = t(t < T);
[qs, vs, fs_] = rr_states(ts, T, I, nrm, alpha, wd, phi);
tg = linspace(0, T, 600);
[qg, vg, fg] = rr_states(tg, T, I, nrm, alpha, wd, phi);

C.model = 'prescribed Hertz pulse (rigid-target contact)';
C.tc = T;
C.q0 = q0;
C.v0 = v0;
C.ts = ts;
C.qs = qs;
C.vs = vs;
C.fs = fs_;
C.tg = tg;
C.qg = qg;
C.vg = vg;
C.fg = fg;
C.plate_v = phi * vg;
C.plate_w = phi * qg;
C.hammer_v = zeros(1, 0);
C.hammer_x = zeros(1, 0);
C.impulse = I;
C.peak_force = H.peak_force;
C.hammer_exit_v = NaN;
C.plate_exit_v = phi * v0';
C.relative_separation_v = NaN;
C.effective_restitution = params.hammer.restitution;
C.notes = {sprintf('prescribed Hertz pulse: %.4g us long, %.4g N peak, %.4g mN s', T * 1e6, H.peak_force, I * 1e3)};
end

function [q, v, f] = rr_states(times, T, I, nrm, alpha, wd, phi)
% modal states at the given times within the pulse, by quadrature over [0, time]
nm = numel(alpha);
nt = numel(times);
q = zeros(nm, nt);
v = zeros(nm, nt);
f = zeros(1, nt);
[z, wz] = plate_gauss_legendre(32);
for j = 1:nt
    tj = times(j);
    if tj <= 0
        continue
    end
    s = 0.5 * tj * (z + 1);
    ws = 0.5 * tj * wz;
    force = I * sin(pi * s / T).^1.5 / nrm;
    u = tj - s;
    E = exp(-alpha' * u);
    SW = sin(wd' * u);
    CW = cos(wd' * u);
    h = E .* SW ./ wd';
    hd = E .* (CW - (alpha' ./ wd') .* SW);
    q(:, j) = phi' .* (h * (ws .* force)');
    v(:, j) = phi' .* (hd * (ws .* force)');
    if tj < T
        f(j) = I * sin(pi * tj / T)^1.5 / nrm;
    end
end
end