function result = plate_fd_output(S, Rec)
% PLATE_FD_OUTPUT  The result struct of plate_fd from the record of its time loop.
%
%   result = plate_fd_output(S, Rec), with S from plate_fd_setup and Rec from
%   plate_fd_loop (or its WebGPU port), truncates the record after a blow-up,
%   evaluates the surface stresses at the saved frames and at the probe, sums
%   the listener's delay bins into the radiated pressure, takes the probe
%   spectra and assembles the fields listed in plate_result_schema.  The
%   stresses are computed here, in double precision, from the flexural frames
%   the loop kept, whichever precision the loop ran in; the rigid-body part
%   has no curvature.  Nodes outside the plate are NaN in every field, so
%   they are drawn blank.

P = S.P;
dom = S.dom; mnt = S.mnt;
nx = S.nx; ny = S.ny; dx = S.dx; dy = S.dy; dt = S.dt; fs = S.fs;
nr = S.nr; N = S.N;
frame_steps = S.frame_steps;
notes = [S.notes, Rec.notes];
w_p = Rec.w_p; wxx_p = Rec.wxx_p; wyy_p = Rec.wyy_p; wxy_p = Rec.wxy_p;
C = Rec.C; Fsup = Rec.Fsup; Wf = Rec.Wf; Cf = Rec.Cf; Pbins = Rec.Pbins;

% Complete output samples only: the last one is partly filled, and after a
% blow-up everything past the last finite step is dropped
nsteps_done = S.nsteps;
if Rec.diverged_step > 0
    nsteps_done = Rec.diverged_step - 1;
    keepf = frame_steps < Rec.diverged_step;
    Wf = Wf(:, :, keepf); Cf = Cf(:, keepf); frame_steps = frame_steps(keepf);
end
n = floor(nsteps_done * dt * fs);
n = max(1, min(n, S.na));
t = ((0:n - 1) + 0.5) / fs;
w_p = w_p(1:n); wxx_p = wxx_p(1:n); wyy_p = wyy_p(1:n); wxy_p = wxy_p(1:n);
C = C(:, 1:n); Fsup = Fsup(:, 1:n);
nf = numel(frame_steps);
tf = frame_steps * dt;

% Total departure at the frames, and the stresses from the flexural frames
Wtot = zeros(ny, nx, nf);
VMf = zeros(ny, nx, nf); SPf = zeros(ny, nx, nf);
rect = strcmp(dom.shape, 'rectangle');
for kf = 1:nf
    phi = Wf(:, :, kf);
    if nr > 0
        Wtot(:, :, kf) = phi + reshape(S.R * Cf(:, kf), ny, nx);
    else
        Wtot(:, :, kf) = phi;
    end
    if rect
        [wxx, wyy, wxy] = plate_curvature_free(phi, dx, dy, P.material.nu);
    else
        [wxx, wyy, wxy] = plate_curvature_domain(phi, dom);
    end
    Sf = plate_stress(P, wxx, wyy, wxy);
    VMf(:, :, kf) = Sf.vm; SPf(:, :, kf) = Sf.sp;
end
if isempty(tf)
    env_vm = zeros(ny, nx); env_sp = zeros(ny, nx);
else
    env_vm = max(VMf, [], 3); env_sp = max(abs(SPf), [], 3);
end
out = ~dom.inside;
w0 = reshape(S.w0, ny, nx);
for kf = 1:nf
    tmp = Wf(:, :, kf); tmp(out) = NaN; Wf(:, :, kf) = tmp;
    tmp = Wtot(:, :, kf); tmp(out) = NaN; Wtot(:, :, kf) = tmp;
    tmp = VMf(:, :, kf); tmp(out) = NaN; VMf(:, :, kf) = tmp;
    tmp = SPf(:, :, kf); tmp(out) = NaN; SPf(:, :, kf) = tmp;
end
env_vm(out) = NaN; env_sp(out) = NaN; w0(out) = NaN;

% Audio: sum the delay bins, then trim to the complete samples
na = S.na;
p_audio = zeros(1, na);
for d = 1:S.nbins
    p_audio(d:end) = p_audio(d:end) + Pbins(d, 1:na - d + 1);
end
p_audio = p_audio(1:n);

Sp = plate_stress(P, wxx_p, wyy_p, wxy_p);
[fw, Aw] = plate_spectrum(t, w_p);
[~, Avm] = plate_spectrum(t, Sp.vm);

% Static support forces per support
fstat_sup = zeros(1, S.nsup);
for q = 1:S.nsup
    nd = S.sup_nodes(S.sup_ptr(q):S.sup_ptr(q + 1) - 1);
    if S.sup_kind(q) == 0
        fstat_sup(q) = sum(S.fstat_s(nd));
    elseif S.sup_kind(q) == 1
        fstat_sup(q) = sum(S.fstat_b(nd));
    else
        fstat_sup(q) = sum(S.fstat_t(nd));
    end
end

result.method = 'fd';
result.grid = struct('x', S.x, 'y', S.y);
result.t = t;
result.probe = struct('x', S.x(S.jp), 'y', S.y(S.ip), 'w', w_p, 'sxx', Sp.sxx, 'syy', Sp.syy, 'txy', Sp.txy, 'vm', Sp.vm, 'sp', Sp.sp);
result.frames = struct('t', tf, 'w', Wf, 'w_total', Wtot, 'rigid', Cf, 'vm', VMf, 'sp', SPf);
result.envelope = struct('vm', env_vm, 'sp', env_sp);
result.spectrum = struct('f', fw, 'w', Aw, 'vm', Avm);
result.audio = struct('fs', fs, 'p', p_audio, 't0', S.Lis.t0, 'model', S.Lis.model);
result.rigid = struct('nr', nr, 'names', {mnt.rigid_names}, 'series', C, ...
    'pivot_x', mnt.pivot(1), 'pivot_y', mnt.pivot(2), 'ed', mnt.ed, 'eh', mnt.eh, ...
    'hz', S.mdl.rigid_hz(:)', 'mode_names', {S.mdl.rigid_mode_names});
result.supports = struct('n', S.nsup, 'names', {S.sup_names}, 'units', {S.sup_units}, 'kinds', {S.sup_kinds}, ...
    'series', Fsup, 'static', fstat_sup);
result.static = struct('w0', w0, 'sag', S.mdl.st.sag, 'gx', mnt.gx, 'gy', mnt.gy, 'gz', mnt.gz, 'g', P.g, ...
    'mass', mnt.mass, 'xc', mnt.xc, 'yc', mnt.yc, 'held', double(mnt.held));
if S.mdl.ps.solved
    result.static.Nxx = S.mdl.ps.Nxx; result.static.Nyy = S.mdl.ps.Nyy; result.static.Nxy = S.mdl.ps.Nxy;
end
result.info = struct('nx', nx, 'ny', ny, 'dx', dx, 'dy', dy, 'dt', dt, 'dt_max', S.mdl.dt_max, 'cfl', P.fd.cfl, ...
    'nsteps', S.nsteps, 'ncells', N, 'nunknowns', S.mdl.nunknowns, 'shape', dom.shape, 'mount', mnt.type, ...
    'cell_steps', N * S.nsteps, 'fs', fs, ...
    'stable', Rec.diverged_step == 0, 'diverged_step', Rec.diverged_step, ...
    'loop', Rec.backend, 'loop_s', Rec.loop_s, 'wall_s', toc(S.t_start), 'notes', {notes});
end