function S = plate_fd_setup(params)
% PLATE_FD_SETUP  Everything the explicit finite-difference time loop needs, computed once.
%
%   S = plate_fd_setup(params) does the part of plate_fd that comes before the
%   time stepping: the derived constants (plate_setup), the grid and time step
%   (plate_fd_plan), the plate's footprint on the grid (plate_domain), the load
%   footprint and force history, the probe node, the frame schedule, the
%   listener geometry and, whenever the loop needs them, the sparse operators
%   and the rigid-body projector (plate_matrices, plate_rigid_modes).
%   plate_fd_loop advances the plate with S, and plate_fd_output turns what the
%   loop recorded into the result struct; plate_fd calls the three in turn.
%
%   The split lets a host run the loop somewhere else while the setup and the
%   output stay this code: the lab has a WebGPU port of plate_fd_loop that reads
%   exactly the fields below and hands its record to plate_fd_output.
%
%   Two forms of the spatial operators.  A rectangle with supported or clamped
%   edges uses the dense stencils on the whole ny x nx array (plate_pad,
%   plate_laplacian).  Everything else, free edges on any shape and every
%   condition on an ellipse, uses the sparse operators of plate_matrices on the
%   unknowns, embedded in vectors over all N = nx ny nodes so that the loop
%   keeps one array layout: rows and columns of nodes that are not unknowns are
%   zero, their node weight is one, and S.keep zeroes them after every step.
%
%   Fields.
%     dom                           the footprint (plate_domain)
%     nx, ny, dx, dy, x, y, X, Y    grid; w(i, j) lives at (x(j), y(i))
%     dt, nsteps, tt, F, fs         time step, step count, tt = (0:nsteps) dt,
%                                   force F(tt), audio sample rate
%     s                             load footprint, ny x nx [1/m^2], zero off the plate
%     ip, jp                        row and column of the probe node
%     frame_steps, nf               steps at which the field is kept
%     Lis, bins, wgt, na            plate_listener geometry; delay bin and
%                                   weight of every node in column order; number
%                                   of audio samples to accumulate
%     free                          the edges are free (with or without the pin)
%     nunknowns                     nodes in the system (the pinned node included)
%     sparse                        the loop uses Lm, Km, mvec, keep, Rb, Qb
%     Lm, Km, mvec                  N x N Laplacian and stiffness, N node weights
%     keep                          ny x nx, 1 at the unknowns and 0 elsewhere
%                                   (off the plate, on a supported ring, at the pin)
%     nr, Rb, Qb                    rigid-body motions to remove (0, 2 or 3):
%                                   N x nr basis and projector, zero off the plate
%     fixed                         node index of the pinned node, or empty
%     c_old, c_new, k2dt2, s1dt2,   constants of the update
%     fdt2, blowup                    wn = c_new (2 w - c_old wm - k2dt2 B
%                                          + s1dt2 (Lw - Lw_m) + fdt2 F(n) s)
%                                   with c_old = 1 - sigma0 dt, c_new =
%                                   1 / (1 + sigma0 dt), k2dt2 = kappa^2 dt^2,
%                                   s1dt2 = 2 sigma1 dt, fdt2 = dt^2 / mu; blowup
%                                   is the deflection that declares the run
%                                   unstable
%     notes                         remarks about this particular run
%     t_start                       tic at entry, for the wall time

S.t_start = tic;
P = plate_setup(params);
plan = plate_fd_plan(params);
nx = plan.nx; ny = plan.ny; dx = plan.dx; dy = plan.dy; dt = plan.dt; nsteps = plan.nsteps;
dom = plate_domain(P.plate, nx, ny);
x = dom.x; y = dom.y; X = dom.X; Y = dom.Y;     % w(i, j) lives at (x(j), y(i))
N = nx * ny;
fs = P.output.fs;
notes = {};
if P.fd.cfl > 1
    notes{end + 1} = sprintf('cfl = %g exceeds the stability limit; expect the solution to blow up', P.fd.cfl);
end
rect = strcmp(dom.shape, 'rectangle');
free = any(strcmp(P.bc, {'free', 'pinned-center'}));
use_sparse = free || ~rect;
if free
    unknown = dom.inside;
else
    unknown = dom.interior;
end

% Load on the unknowns only: a supported ring and the outside carry none
s = plate_force_footprint(P.excitation, X, Y, dx * dy, unknown);
tt = (0:nsteps) * dt;
F = plate_force_time(P.excitation, tt);

% Probe: the nearest node whose eight neighbours belong to the plate, so the
% curvatures there are central differences
d2 = (X - P.probe.x).^2 + (Y - P.probe.y).^2;
d2(~dom.full3) = Inf;
[~, kp] = min(d2(:));
ip = mod(kp - 1, ny) + 1;
jp = floor((kp - 1) / ny) + 1;
if abs(x(jp) - P.probe.x) > 0.5 * dx || abs(y(ip) - P.probe.y) > 0.5 * dy
    notes{end + 1} = sprintf('probe moved to the nearest interior node (%.4g, %.4g) m', x(jp), y(ip));
end

% Frames at the steps nearest to nframes equally spaced times in [0, duration]
frame_steps = unique(max(1, min(nsteps, round(linspace(0, P.fd.duration, P.output.nframes) / dt))));

% Listener: every grid node is a patch; the acceleration of a step is added into
% the audio sample that contains it, delayed by the node's travel time.
Lis = plate_listener(P.listener, X, Y, dx * dy, fs);

% Sparse operators on the unknowns, embedded over all N nodes, and the
% projector onto the complement of the rigid-body motions
fixed = zeros(0, 1);
Lm = []; Km = []; mvec = []; Rb = []; Qb = []; nr = 0;
keep = double(unknown);
nunknowns = nnz(unknown);                % the pinned node counts: it is in the system, held at zero
if use_sparse
    [Lu, Ku, Mu, idx, fixed_u] = plate_matrices(dom, P.bc, P.material.nu);
    [ki, kj, kv] = find(Ku);
    Km = sparse(idx(ki), idx(kj), kv, N, N);
    [li, lj, lv] = find(Lu);
    Lm = sparse(idx(li), idx(lj), lv, N, N);
    mvec = ones(N, 1);
    mvec(idx) = full(diag(Mu));
    if ~isempty(fixed_u)
        fixed = idx(fixed_u);
        keep(fixed) = 0;
    end
    if free
        if isempty(fixed)
            [Ru, Qu] = plate_rigid_modes(X(idx), Y(idx), mvec(idx));
            notes{end + 1} = 'free edges: the rigid-body translation and tilts are removed from the displacement after every step';
        else
            [Ru, Qu] = plate_rigid_modes(X(idx), Y(idx), mvec(idx), [X(fixed) Y(fixed)]);
            notes{end + 1} = sprintf('pinned at the node (%.4g, %.4g) m nearest the centre; the two rigid rotations about it are removed after every step', X(fixed), Y(fixed));
        end
        nr = size(Ru, 2);
        Rb = zeros(N, nr); Rb(idx, :) = Ru;
        Qb = zeros(N, nr); Qb(idx, :) = Qu;
    end
end
if ~rect
    notes{end + 1} = sprintf('%s on a staircase of %d grid nodes; the edge is first order in the grid spacing', dom.shape, nnz(dom.inside));
end

S.P = P; S.plan = plan; S.dom = dom;
S.nx = nx; S.ny = ny; S.dx = dx; S.dy = dy; S.x = x; S.y = y; S.X = X; S.Y = Y;
S.dt = dt; S.nsteps = nsteps; S.tt = tt; S.F = F; S.fs = fs;
S.s = s;
S.ip = ip; S.jp = jp;
S.frame_steps = frame_steps; S.nf = numel(frame_steps);
S.Lis = Lis; S.bins = Lis.delay(:) + 1; S.wgt = Lis.weight(:);
S.na = ceil(nsteps * dt * fs) + Lis.nbins + 1;
S.free = free; S.sparse = use_sparse; S.fixed = fixed; S.nunknowns = nunknowns;
S.Lm = Lm; S.Km = Km; S.mvec = mvec; S.keep = keep; S.nr = nr; S.Rb = Rb; S.Qb = Qb;
S.c_old = 1 - P.damping.sigma0 * dt;
S.c_new = 1 / (1 + P.damping.sigma0 * dt);
S.k2dt2 = P.kappa^2 * dt^2;
S.s1dt2 = 2 * P.damping.sigma1 * dt;
S.fdt2 = dt^2 / P.mu;
S.blowup = 1e3 * P.plate.h;              % deflection beyond which the run is declared unstable
S.notes = notes;
end