function R = plate_fd_loop(S)
% PLATE_FD_LOOP  The explicit time stepping of plate_fd.
%
%   R = plate_fd_loop(S), with S from plate_fd_setup, advances
%
%     wn = c_new * (2 w - c_old wm - k2dt2 B + s1dt2 (Lw - Lw_m) + fdt2 F(n) s)
%
%   from rest for S.nsteps steps, where Lw is the five-point Laplacian of the
%   current field, B the biharmonic (the Laplacian applied twice, the 13-point
%   stencil) and Lw_m the Laplacian of the previous field.  On a rectangle with
%   supported or clamped edges the operators are the dense stencils with one
%   layer of ghost nodes (plate_pad) and the edge nodes are held at zero.
%   Otherwise (S.sparse) they are the sparse operators S.Lm, S.Km of
%   plate_matrices on every node, divided by the node weights S.mvec; the
%   nodes that are not unknowns (off the plate, on a supported ring, the pin)
%   are zeroed with S.keep, and the S.nr rigid-body motions of a free plate
%   are projected out of the new field with S.Rb, S.Qb.
%
%   What the loop records, and returns in R:
%     w_p, wxx_p, wyy_p, wxy_p   displacement and curvatures at the probe node,
%                                1 x (nsteps + 1); entry n + 1 is step n (t = n dt)
%     Wf                         the field at the frame steps, ny x nx x nf
%     Pbins                      radiated pressure per listener delay bin and
%                                audio sample, nbins x na: the acceleration of
%                                every node at every step, weighted by
%                                S.wgt and box-averaged onto audio samples
%     diverged_step              0, or the step at which max |w| reached
%                                S.blowup (the run stops there)
%     notes                      remarks about the run
%     backend                    'numbl': where the loop ran
%     loop_s                     wall time of the loop by tic/toc
%
%   The lab has a WebGPU port of this function (packages/core/src/gpu/fdLoop.ts)
%   that computes the same record from the same S, in single precision.  It
%   carries the state as (w, w - wm) rather than (w, wm), which is the same
%   arithmetic rearranged so that the rounding of w does not feed back into the
%   velocity; keep the two implementations in step.

t_loop = tic;
nx = S.nx; ny = S.ny; dx = S.dx; dy = S.dy; dt = S.dt; nsteps = S.nsteps; fs = S.fs;
bc = S.P.bc; use_sparse = S.sparse; nr = S.nr;
s = S.s; F = S.F; ip = S.ip; jp = S.jp;
frame_steps = S.frame_steps; nf = S.nf;
bins = S.bins; wgt = S.wgt; nbins = S.Lis.nbins; na = S.na;
c_old = S.c_old; c_new = S.c_new; k2dt2 = S.k2dt2; s1dt2 = S.s1dt2; fdt2 = S.fdt2; blowup = S.blowup;
if use_sparse
    Lm = S.Lm; Km = S.Km; mvec = S.mvec; keep = S.keep; Rb = S.Rb; Qb = S.Qb;
end

% State
w = zeros(ny, nx); wm = zeros(ny, nx); Lw_m = zeros(ny, nx);
w_p = zeros(1, nsteps + 1); wxx_p = zeros(1, nsteps + 1); wyy_p = zeros(1, nsteps + 1); wxy_p = zeros(1, nsteps + 1);
Wf = zeros(ny, nx, nf);
Pbins = zeros(nbins, na);
diverged_step = 0;
notes = {};
kf = 1;

for nstep = 1:nsteps
    if use_sparse
        Lw = reshape((Lm * w(:)) ./ mvec, ny, nx);
        B = reshape((Km * w(:)) ./ mvec, ny, nx);
    else
        Lw = plate_laplacian(plate_pad(w, bc), dx, dy);
        B = zeros(ny, nx);
        B(2:end - 1, 2:end - 1) = (Lw(3:end, 2:end - 1) - 2 * Lw(2:end - 1, 2:end - 1) + Lw(1:end - 2, 2:end - 1)) / dy^2 ...
                                + (Lw(2:end - 1, 3:end) - 2 * Lw(2:end - 1, 2:end - 1) + Lw(2:end - 1, 1:end - 2)) / dx^2;
    end
    wn = c_new * (2 * w - c_old * wm - k2dt2 * B + s1dt2 * (Lw - Lw_m) + fdt2 * F(nstep) * s);
    if use_sparse
        wn = wn .* keep;
        if nr > 0
            wn = wn - reshape(Rb * (Qb' * wn(:)), ny, nx);
        end
    else
        wn(1, :) = 0; wn(end, :) = 0; wn(:, 1) = 0; wn(:, end) = 0;
    end

    % Probe displacement and curvatures from the 3 x 3 neighbourhood
    w_p(nstep + 1) = wn(ip, jp);
    wxx_p(nstep + 1) = (wn(ip, jp + 1) - 2 * wn(ip, jp) + wn(ip, jp - 1)) / dx^2;
    wyy_p(nstep + 1) = (wn(ip + 1, jp) - 2 * wn(ip, jp) + wn(ip - 1, jp)) / dy^2;
    wxy_p(nstep + 1) = (wn(ip + 1, jp + 1) - wn(ip + 1, jp - 1) - wn(ip - 1, jp + 1) + wn(ip - 1, jp - 1)) / (4 * dx * dy);

    % Radiated pressure: acceleration at this step, weighted and delayed per node
    acc = (wn - 2 * w + wm) / dt^2;
    ia = floor(nstep * dt * fs) + 1;
    Pbins(:, ia) = Pbins(:, ia) + accumarray(bins, wgt .* acc(:), [nbins 1]) * (dt * fs);

    if kf <= nf && nstep == frame_steps(kf)
        Wf(:, :, kf) = wn;
        kf = kf + 1;
    end

    peak = max(abs(wn(:)));
    if ~(peak < blowup)                  % also catches NaN and Inf
        diverged_step = nstep;
        notes{end + 1} = sprintf('solution diverged at step %d (t = %.3g s); results are truncated there', nstep, nstep * dt);
        break
    end
    wm = w; w = wn; Lw_m = Lw;
end

R.w_p = w_p; R.wxx_p = wxx_p; R.wyy_p = wyy_p; R.wxy_p = wxy_p;
R.Wf = Wf;
R.Pbins = Pbins;
R.diverged_step = diverged_step;
R.notes = notes;
R.backend = 'numbl';
R.loop_s = toc(t_loop);
end