function S = plate_fd_setup(params)
% PLATE_FD_SETUP  Everything the explicit time loop needs, computed once.
%
%   S = plate_fd_setup(params) takes the model (plate_model: operators,
%   statics, rigid-body basis, time step) and adds what a run needs: the load
%   footprint and force history, the probe node, the frame schedule, the
%   listener geometry, the initial state and the constants of the update.
%   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.
%
%   State.  The loop carries the departure from the static equilibrium w0 as
%   rigid-body coordinates c (nr of them, plate_rigid_basis) and a flexural
%   field phi on the grid; the displacement is w = w0 + R c + phi.  A release
%   (excitation.type = 'release') starts from the static deflection under a
%   force at the pull point, scaled to the asked displacement there and split
%   into c0 and phi0; a strike starts from rest.
%
%   Fields the loop reads (indices are 1-based, arrays over all N = nx ny nodes
%   in column-major order unless said otherwise).
%     nx, ny, dx, dy, dt, nsteps, fs, na   grid, step, count, sample rate and
%                                          number of output samples
%     keep                                 N x 1: 1 at the unknowns, 0 elsewhere
%     Mvec                                 N x 1 node masses (1 where keep is 0)
%     KB, KG, Lm                           sparse N x N: bending, prestress plus
%                                          linear springs, Laplacian per node area
%     Gs, Gscale, WB, WG, WL, ng, inva     the same operators in factored form
%                                          (see below), for the port
%     nr, R, Q, AR, Krr, Ginv              rigid-body basis (N x nr), projector,
%                                          coupling forces, reduced stiffness and
%                                          inverse mass (nr x nr)
%     kb, cb, gb, kt, ct, gt               bottom and top contacts: stiffness
%                                          [N/m], dashpot [N s/m], compression at
%                                          rest [m], per node
%     cs, ks                               linear dashpots and springs (the springs
%                                          for the force report; they are in KG)
%     fstat_s, fstat_b, fstat_t            static force of the springs, the bottom
%                                          and the top contacts per node
%     sup_ptr, sup_nodes, nsup, sup_kind   supports as a CSR list of node indices,
%                                          each reporting one kind of force
%                                          (0 spring, 1 bottom, 2 top)
%     sa                                   force per node per unit F(t) [N/N]
%     F                                    force history at the nsteps + 1 levels
%     ip, jp                               probe row and column
%     frame_steps, nf                      steps at which the field is kept
%     lis_node, lis_bin, lis_wgt, nbins    listener entries (plate_listener)
%     c0, phi0                             initial state
%     wsamp                                na x 1 weight of a step in its sample
%     c_old, c_new, s1f, LR, dt2           constants of the update: the sigma0 loss,
%                                          the sigma1 loss as a nodal force on the
%                                          total velocity (with Lm R for the rigid part)
%     blowup                               deflection that declares the run unstable
%     notes, t_start

S.t_start = tic;
mdl = plate_model(params);
P = mdl.P; dom = mdl.dom; mnt = mdl.mnt;
nx = mdl.nx; ny = mdl.ny; dx = mdl.dx; dy = mdl.dy; N = mdl.N;
dt = mdl.dt; nsteps = mdl.nsteps;
X = dom.X; Y = dom.Y; x = dom.x; y = dom.y;
unk = mdl.unk;
fs = P.output.fs;
notes = mdl.notes;

% Load on the unknowns: force per node per unit of F(t)
s = plate_force_footprint(P.excitation, X, Y, dx * dy, unk);
sa = s(:) .* mdl.a(:);
tt = (0:nsteps) * dt;
F = plate_force_time(P.excitation, tt);

% Probe: the nearest node whose eight neighbours belong to the plate
d2 = (X - P.probe.x).^2 + (Y - P.probe.y).^2;
d2(~dom.full3 | mnt.held) = 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, and the output samples: each series value is the mean over the
% steps whose time falls in the sample, so the weight of a step is one over
% the number of steps its sample contains
Lis = plate_listener(P.listener, dom, dx * dy, fs);
na = floor(nsteps * dt * fs) + 1;
ia_of = floor((0:nsteps - 1) * dt * fs) + 1;
cnt = accumarray(ia_of', 1, [na 1]);
wsamp = 1 ./ max(cnt, 1);

% Initial state
nr = mdl.nr;
c0 = zeros(nr, 1);
phi0 = zeros(N, 1);
if strcmp(P.excitation.type, 'release')
    e = P.excitation;
    sp = plate_force_footprint(e, X, Y, dx * dy, unk);
    f = sp(:) .* mdl.a(:);
    A = mdl.st.Alin;
    if mnt.free_rotation && nr > 0
        % the spin about the thread has no stiffness to speak of: take the
        % spin moment out of the load and pin the spin at one node
        rh = mdl.R(:, 3);
        Mr = mdl.Mvec .* rh;
        f = f - (rh' * f) / (rh' * Mr) * Mr;
        [~, kf] = max(abs(rh));
        A = A + sparse(kf, kf, mdl.Mvec(kf) * 1e6, N, N);
        notes{end + 1} = 'the pull is applied without twisting the plate about its thread';
    end
    iu = find(unk(:));
    wr = zeros(N, 1);
    [Rc, p] = chol(A(iu, iu));
    if p > 0
        error('plate_fd_setup: no static deflection to release from; the plate is not held stably');
    end
    wr(iu) = Rc \ (Rc' \ f(iu));
    wP = (sp(:) .* mdl.a(:))' * wr;                       % displacement at the pull point
    if ~(abs(wP) > 0)
        error('plate_fd_setup: the pull point does not move; is it on a held part of the plate?');
    end
    wr = wr * (e.displacement / wP);
    if nr > 0
        c0 = mdl.Q' * wr;
        if mnt.free_rotation
            c0(3) = 0;
        end
        phi0 = wr - mdl.R * c0;
    else
        phi0 = wr;
    end
    phi0 = phi0 .* double(unk(:));
    notes{end + 1} = sprintf('released from a static pull of %.3g mm at (%.4g, %.4g) m, which needed %.3g N', e.displacement * 1e3, e.x, e.y, e.displacement / wP);
end

% The operators in factored form, for the WebGPU port.  Every operator of the
% loop is a sum of terms G_q' W G_q with G_q a difference stencil, and the
% port applies them that way: strains e = G phi first (each an exact integer
% stencil times a scale), then the weights, then G' back to the nodes.  The
% assembled matrices lose their smooth modes in single precision (the 13-point
% biharmonic of a smooth field is the small residual of large cancelling
% terms, and the rounding of the field alone is a few percent of it at low
% wavenumbers); the factored form keeps them to about 1e-5.  Gs stacks
% Dxx, Dyy, Dxc, Dyc, Dxy, Bx, By, Bxc, Byc; Gscale is the scale of each row
% (the stencil divided by it is integer); WB, WG, WL are the symmetric weight
% matrices with KB = Gs' WB Gs, KG - Ks = Gs' WG Gs and Lm = diag(1/a) Gs' WL Gs
% on the unknowns.
ops = mdl.ops;
ni = ops.n; idx = ops.idx;
Gu = [ops.Dxx; ops.Dyy; ops.Dxc; ops.Dyc; ops.Dxy; ops.Bx; ops.By; ops.Bxc; ops.Byc];
nrow = [size(ops.Dxx, 1), size(ops.Dyy, 1), size(ops.Dxc, 1), size(ops.Dyc, 1), size(ops.Dxy, 1), ...
        size(ops.Bx, 1), size(ops.By, 1), size(ops.Bxc, 1), size(ops.Byc, 1)];
scl = [1 / dx^2, 1 / dy^2, 1 / dx^2, 1 / dy^2, 1 / (dx * dy), 1 / dx, 1 / dy, 1 / (2 * dx), 1 / (2 * dy)];
Gscale = zeros(sum(nrow), 1);
off = cumsum([0 nrow]);
for q = 1:9
    Gscale(off(q) + 1:off(q + 1)) = scl(q);
end
ng = sum(nrow);
[gi, gj, gv] = find(Gu);
Gs = sparse(gi, idx(gj), gv, ng, N);
D = ops.D; nu = ops.nu; area = ops.area;
o = off;
WB = sparse(ng, ng);
WB = WB + sparse(o(1) + (1:nrow(1)), o(1) + (1:nrow(1)), D * ops.mkx, ng, ng);
WB = WB + sparse(o(2) + (1:nrow(2)), o(2) + (1:nrow(2)), D * ops.mky, ng, ng);
WB = WB + sparse(o(3) + (1:nrow(3)), o(4) + (1:nrow(4)), D * nu * ops.mkc, ng, ng);
WB = WB + sparse(o(4) + (1:nrow(4)), o(3) + (1:nrow(3)), D * nu * ops.mkc, ng, ng);
WB = WB + sparse(o(5) + (1:nrow(5)), o(5) + (1:nrow(5)), D * 2 * (1 - nu) * area * ones(nrow(5), 1), ng, ng);
WG = sparse(ng, ng);
if mdl.ps.solved
    WG = WG + sparse(o(6) + (1:nrow(6)), o(6) + (1:nrow(6)), mdl.ps.nxx, ng, ng);
    WG = WG + sparse(o(7) + (1:nrow(7)), o(7) + (1:nrow(7)), mdl.ps.nyy, ng, ng);
    WG = WG + sparse(o(8) + (1:nrow(8)), o(9) + (1:nrow(9)), mdl.ps.nxy, ng, ng);
    WG = WG + sparse(o(9) + (1:nrow(9)), o(8) + (1:nrow(8)), mdl.ps.nxy, ng, ng);
end
WL = sparse(o(6) + (1:nrow(6)), o(6) + (1:nrow(6)), -ops.wx, ng, ng) + sparse(o(7) + (1:nrow(7)), o(7) + (1:nrow(7)), -ops.wy, ng, ng);

% Supports as a CSR list
sup = mnt.supports;
nsup = numel(sup);
sup_ptr = ones(nsup + 1, 1);
sup_nodes = zeros(0, 1);
for k = 1:nsup
    nd = sup(k).nodes(:);
    sup_nodes = [sup_nodes; nd];
    sup_ptr(k + 1) = sup_ptr(k) + numel(nd);
end
sup_names = cell(1, nsup); sup_units = cell(1, nsup); sup_kinds = cell(1, nsup);
sup_kind = zeros(nsup, 1);                      % 0 spring, 1 bottom contact, 2 top contact
for k = 1:nsup
    sup_names{k} = sup(k).name; sup_units{k} = sup(k).unit; sup_kinds{k} = sup(k).kind;
    switch sup(k).kind
        case 'spring'
            sup_kind(k) = 0;
        case 'bottom'
            sup_kind(k) = 1;
        case 'top'
            sup_kind(k) = 2;
        otherwise
            error('plate_fd_setup: unknown support kind "%s"', sup(k).kind);
    end
end

S.mdl = mdl; S.P = P; S.dom = dom; S.mnt = mnt;
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.N = N;
S.dt = dt; S.nsteps = nsteps; S.tt = tt; S.F = F; S.fs = fs; S.na = na; S.wsamp = wsamp;
S.keep = double(unk(:));
S.Mvec = mdl.Mvec; S.a = mdl.a(:);
S.KB = mdl.Kb; S.KG = mdl.Kg + mdl.Ks; S.Lm = mdl.Lm;
S.Gs = Gs; S.Gscale = Gscale; S.WB = WB; S.WG = WG; S.WL = WL; S.ng = ng;
S.inva = double(unk(:)) ./ max(mdl.a(:), eps);
S.nr = nr; S.R = mdl.R; S.Q = mdl.Q; S.AR = mdl.AR; S.Krr = mdl.Krr; S.Ginv = mdl.Ginv;
S.kb = mnt.kb(:); S.cb = mnt.cb(:); S.gb = mdl.st.gb;
S.kt = mnt.kt(:); S.ct = mnt.ct(:); S.gt = mdl.st.gt;
S.cs = mnt.cs(:); S.ks = mnt.ks(:);
S.fstat_s = mdl.st.fstat_s; S.fstat_b = mdl.st.fstat_b; S.fstat_t = mdl.st.fstat_t;
S.w0 = mdl.st.w0;
S.sup_ptr = sup_ptr; S.sup_nodes = sup_nodes; S.nsup = nsup; S.sup_kind = sup_kind;
S.sup_names = sup_names; S.sup_units = sup_units; S.sup_kinds = sup_kinds;
S.sa = sa; S.s = s;
S.ip = ip; S.jp = jp;
S.frame_steps = frame_steps; S.nf = numel(frame_steps);
S.Lis = Lis; S.lis_node = Lis.node; S.lis_bin = Lis.delay + 1; S.lis_wgt = Lis.weight; S.nbins = Lis.nbins;
S.c0 = c0; S.phi0 = phi0;
S.c_old = 1 - P.damping.sigma0 * dt;
S.c_new = 1 / (1 + P.damping.sigma0 * dt);
S.s1f = 2 * P.damping.sigma1 * mdl.Mvec / dt;    % the sigma1 loss as a nodal force per unit Lm(w - wm)
S.LR = mdl.Lm * mdl.R;                          % Laplacian per node of the rigid basis, N x nr
S.dt2 = dt^2;
S.blowup = 1e3 * P.plate.h;
S.notes = notes;
end