function Rec = plate_fd_loop(S)
% PLATE_FD_LOOP  The explicit time stepping of plate_fd.
%
%   Rec = plate_fd_loop(S), with S from plate_fd_setup, advances the departure
%   from the static equilibrium, w = R c + phi, for S.nsteps steps: the rigid-
%   body coordinates c as nr scalar equations and the flexural field phi on
%   the grid,
%
%     cacc = Ginv (R' f - Krr c - AR' phi)
%     cn   = c_new (2 c - c_old cm + dt^2 cacc)
%     acc  = (f - KB phi - KG phi - AR c) / M - R cacc
%     phin = c_new (2 phi - c_old phim + dt^2 acc)
%
%   where f is the nodal force at this step: the strike F(n) sa, the contact
%   and dashpot forces about rest, and the sigma1 loss 2 sigma1 M Lm (w - wm) / dt
%   on the backward-difference velocity of the total field w = R c + phi.
%   c_old = 1 - sigma0 dt and c_new = 1 / (1 + sigma0 dt) carry the
%   frequency-independent loss.  Every force enters both parts through f, so
%   the sum R cn + phin is exactly the leapfrog step of the whole system and
%   the losses stay dissipative; the term - R cacc keeps phi free of
%   rigid-body content, so the two parts stay what they are called.  (A
%   sigma1 loss applied to phi alone is not symmetric on the whole system and
%   was found to pump the swing of a hanging plate slowly.)
%   The scheme is second order in space and time (Bilbao, Numerical Sound
%   Synthesis, 2009, ch. 12).
%
%   Contacts.  A bottom contact with compression gb at rest pushes up by
%   kb (gb - w) while gb - w > 0, a top contact with compression gt pushes
%   down by kt (gt + w) while gt + w > 0; the static parts are subtracted so
%   f is zero at rest, and the dashpots act on the backward-difference
%   velocity while the contact is engaged.
%
%   What the loop records, and returns in Rec.  The series are averaged onto
%   the output samples (S.na of them at S.fs): each step adds its value,
%   weighted by S.wsamp, to the sample that contains its time, so a sample is
%   the mean over its steps.
%     w_p, wxx_p, wyy_p, wxy_p   displacement (total departure) and curvatures
%                                (flexural) at the probe node, 1 x na
%     C                          rigid-body coordinates, nr x na
%     Fsup                       force of each support, nsup x na [N]: the
%                                sum over its nodes of the force of its kind
%                                (spring, bottom contact or top contact),
%                                static part included
%     Pbins                      radiated pressure per listener delay bin and
%                                sample, nbins x na
%     Wf, Cf                     phi at the frame steps, ny x nx x nf, and c there, nr x nf
%     diverged_step              0, or the step at which the field passed
%                                S.blowup (the run stops there)
%     notes, backend, loop_s
%
%   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 phi as (phi, phi - phim) and c likewise, the same arithmetic
%   rearranged so that the rounding of the field 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; na = S.na;
N = S.N; nr = S.nr;
KB = S.KB; KG = S.KG; Lm = S.Lm; Mvec = S.Mvec; keep = S.keep;
R = S.R; AR = S.AR; Krr = S.Krr; Ginv = S.Ginv;
kb = S.kb; cb = S.cb; gb = S.gb; kt = S.kt; ct = S.ct; gt = S.gt; cs = S.cs; ks = S.ks;
fstat_s = S.fstat_s; fstat_b = S.fstat_b; fstat_t = S.fstat_t;
sa = S.sa; F = S.F; ip = S.ip; jp = S.jp; kp = ip + (jp - 1) * ny;
frame_steps = S.frame_steps; nf = S.nf;
ln = S.lis_node; lb = S.lis_bin; lw = S.lis_wgt; nbins = S.nbins;
sup_ptr = S.sup_ptr; sup_nodes = S.sup_nodes; nsup = S.nsup; sup_kind = S.sup_kind;
c_old = S.c_old; c_new = S.c_new; s1f = S.s1f; LR = S.LR; dt2 = S.dt2; blowup = S.blowup;
wsamp = S.wsamp;
gb0 = max(gb, 0); gt0 = max(gt, 0);
hasb = any(kb > 0); hast = any(kt > 0); hascs = any(cs > 0);
rmax = max(abs(R), [], 1);

% State
phi = S.phi0; phim = phi; c = S.c0; cm = c;
Lphi_m = Lm * phi;
w_p = zeros(1, na); wxx_p = zeros(1, na); wyy_p = zeros(1, na); wxy_p = zeros(1, na);
C = zeros(nr, na); Fsup = zeros(nsup, na);
Pbins = zeros(nbins, na);
Wf = zeros(ny, nx, nf); Cf = zeros(nr, nf);
diverged_step = 0;
notes = {};
kf = 1;

for nstep = 1:nsteps
    ia = floor((nstep - 1) * dt * fs) + 1;
    ws = wsamp(ia);
    Lphi = Lm * phi;
    if nr > 0
        w = R * c + phi;
        v = ((phi - phim) + R * (c - cm)) / dt;
    else
        w = phi;
        v = (phi - phim) / dt;
    end

    % The sigma1 loss on the total velocity, as a nodal force: 2 sigma1 M Lm (w - wm) / dt
    dLw = Lphi - Lphi_m;
    if nr > 0
        dLw = dLw + LR * (c - cm);
    end
    f = F(nstep) * sa + s1f .* dLw;

    % Contact and dashpot forces about rest
    fsS = fstat_s - ks .* w;                 % springs: their static force plus the departure
    fsB = fstat_b; fsT = fstat_t;
    if hasb
        compb = gb - w;
        actb = kb > 0 & compb > 0;
        fb = kb .* max(compb, 0) - kb .* gb0 - (cb .* v) .* actb;
        f = f + fb; fsB = fsB + fb;
    end
    if hast
        compt = gt + w;
        actt = kt > 0 & compt > 0;
        ft = -kt .* max(compt, 0) + kt .* gt0 - (ct .* v) .* actt;
        f = f + ft; fsT = fsT + ft;
    end
    if hascs
        fl = -cs .* v;
        f = f + fl; fsS = fsS + fl;
    end

    % Series of the current state, box-averaged onto the output samples
    w_p(ia) = w_p(ia) + w(kp) * ws;
    wxx_p(ia) = wxx_p(ia) + (phi(kp + ny) - 2 * phi(kp) + phi(kp - ny)) / dx^2 * ws;
    wyy_p(ia) = wyy_p(ia) + (phi(kp + 1) - 2 * phi(kp) + phi(kp - 1)) / dy^2 * ws;
    wxy_p(ia) = wxy_p(ia) + (phi(kp + ny + 1) - phi(kp + ny - 1) - phi(kp - ny + 1) + phi(kp - ny - 1)) / (4 * dx * dy) * ws;
    if nr > 0
        C(:, ia) = C(:, ia) + c * ws;
    end
    for q = 1:nsup
        nd = sup_nodes(sup_ptr(q):sup_ptr(q + 1) - 1);
        if sup_kind(q) == 0
            fq = fsS(nd);
        elseif sup_kind(q) == 1
            fq = fsB(nd);
        else
            fq = fsT(nd);
        end
        Fsup(q, ia) = Fsup(q, ia) + sum(fq) * ws;
    end

    % Rigid-body coordinates
    if nr > 0
        cacc = Ginv * (R' * f - Krr * c - AR' * phi);
        cn = c_new * (2 * c - c_old * cm + dt2 * cacc);
        acc = (f - KB * phi - KG * phi - AR * c) ./ Mvec - R * cacc;
    else
        cn = c;
        acc = (f - KB * phi - KG * phi) ./ Mvec;
    end
    phin = c_new * (2 * phi - c_old * phim + dt2 * acc);
    phin = phin .* keep;

    % Radiated pressure: acceleration of the total field at this step, weighted and delayed per entry
    accw = ((phin - phi) - (phi - phim)) / dt2;
    if nr > 0
        accw = accw + R * (((cn - c) - (c - cm)) / dt2);
    end
    Pbins(:, ia) = Pbins(:, ia) + accumarray(lb, lw .* accw(ln), [nbins 1]) * ws;

    if kf <= nf && nstep == frame_steps(kf)
        Wf(:, :, kf) = reshape(phin, ny, nx);
        if nr > 0
            Cf(:, kf) = cn;
        end
        kf = kf + 1;
    end

    peak = max(abs(phin));
    if nr > 0
        peak = peak + sum(abs(cn(:)') .* rmax);
    end
    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
    phim = phi; phi = phin; cm = c; c = cn; Lphi_m = Lphi;
end

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