function result = plate_fem_strike(params, M)
% PLATE_FEM_STRIKE  Strike the finite-element modes and render the result.
%
%   result = plate_fem_strike(params, M) takes the modes of plate_fem_modes,
%   interpolates each shape at the strike point and at the probe within the
%   element that holds them (plate_fem_interp), sets the modal damping
%   (plate_damping; the rigid-body modes get 0.01), integrates the coupled
%   hammer/plate contact (plate_fem_contact) and synthesizes the sound, the
%   series and the frames (plate_synthesize).  Only this function needs to
%   run again when the hammer, the strike point, the microphone or the
%   render settings change.

t_start = tic;
if ~strcmp(M.method, 'fem')
    error('plate_fem_strike: the modes were computed by another method (%s)', M.method);
end
P = plate_setup(params);
pl = params.plate;
if abs(params.strike.x) > pl.width / 2 || abs(params.strike.y) > pl.height / 2
    error('plate_fem_strike: the strike point is outside the plate');
end
if abs(params.probe.x) > pl.width / 2 || abs(params.probe.y) > pl.height / 2
    error('plate_fem_strike: the probe is outside the plate');
end
if params.fem.mesh_x ~= M.nx || params.fem.mesh_y ~= M.ny
    error('plate_fem_strike: the modes were computed on a %d x %d mesh, not %d x %d', M.nx, M.ny, params.fem.mesh_x, params.fem.mesh_y);
end

nm = numel(M.freqs);
nn = numel(M.shapes) / nm;
Sh = reshape(M.shapes, nn, nm);
[idx, wts] = plate_fem_interp(params.strike.x, params.strike.y, params);
phi_hit = wts * Sh(idx, :);
[zeta, Q] = plate_damping(M.freqs, M.support_part, params, P);

fs = params.render.fs;
na = round(params.render.duration * fs);
t = (0:na - 1) / fs;
C = plate_fem_contact(params, P, M.freqs, zeta, phi_hit, t);

[ip, wp] = plate_fem_interp(params.probe.x, params.probe.y, params);
D.zeta = zeta;
D.Q = Q;
D.phi_hit = phi_hit;
D.phi_probe = wp * Sh(ip, :);
Kxx = reshape(M.kxx, nn, nm);
Kyy = reshape(M.kyy, nn, nm);
Kxy = reshape(M.kxy, nn, nm);
D.kxx_probe = wp * Kxx(ip, :);
D.kyy_probe = wp * Kyy(ip, :);
D.kxy_probe = wp * Kxy(ip, :);

result = plate_synthesize(params, P, M, C, D);
result.info.strike_wall_s = toc(t_start);
end