function result = plate_fd(params)
% PLATE_FD  Explicit finite differences for the damped Kirchhoff-Love plate.
%
%   result = plate_fd(params) advances
%
%     w_tt + 2 sigma0 w_t - 2 sigma1 lap(w_t) + kappa^2 lap(lap(w)) = F(t) s(x,y) / mu
%
%   on a uniform nx x ny grid (nodes include the edges) with the scheme
%
%     (1 + sigma0 dt) w+ = 2 w - (1 - sigma0 dt) w- - dt^2 kappa^2 L(L w)
%                          + 2 sigma1 dt (L w - L w-) + dt^2 F s / mu
%
%   where w-, w, w+ are consecutive time levels and L is the five-point
%   Laplacian.  Applying L twice gives the 13-point biharmonic stencil.  The
%   scheme is second order in space and time (Bilbao, Numerical Sound
%   Synthesis, 2009, ch. 12).  On a rectangle, simply supported and clamped
%   edges enter through one layer of ghost nodes (plate_pad).  Free edges have
%   no ghost-layer form, so for them L and L(L .) are the sparse energy-method
%   operators of plate_matrices applied to every node, and the rigid-body
%   motions are projected out of the displacement after each step
%   (plate_rigid_modes).  The same path serves the centre-pinned plate, whose
%   pinned node is zeroed after each step, and the ellipse (plate_domain),
%   whose edge is a staircase on the grid with every condition assembled as a
%   sparse operator.
%
%   The step is dt = cfl * dt_max with dt_max from plate_fd_plan; the scheme
%   is unstable for cfl > 1 and the run stops as soon as the field diverges,
%   reporting the step at which it did (info.diverged_step).
%
%   The work is in three stages, each its own file, so that a host can run the
%   time loop elsewhere (the lab has a WebGPU port of plate_fd_loop) while the
%   setup and the output stay this code:
%     S      = plate_fd_setup(params)    grid, step, load, probe, listener, operators
%     R      = plate_fd_loop(S)          the time stepping; records probe, frames, audio bins
%     result = plate_fd_output(S, R)     stresses, audio, spectra, the result struct
%
%   Approximations.  Second-order truncation error in dx, dy, dt; the probe is
%   snapped to the nearest grid node; stress fields are evaluated only at the
%   saved frames, so the envelope is a sample of the true peak; the audio is
%   the Rayleigh integral box-averaged onto audio samples.
%
%   See plate_result_schema for the fields of result.

S = plate_fd_setup(params);
R = plate_fd_loop(S);
result = plate_fd_output(S, R);
end