function result = plate_fd_output(S, R)
% PLATE_FD_OUTPUT  The result struct of plate_fd from the record of its time loop.
%
%   result = plate_fd_output(S, R), with S from plate_fd_setup and R 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 displacement
%   frames the loop kept, whichever precision the loop ran in.  On an ellipse
%   the nodes outside the plate are set to NaN in every field, so they are
%   drawn blank.

P = S.P;
dom = S.dom;
nx = S.nx; ny = S.ny; dx = S.dx; dy = S.dy; dt = S.dt; fs = S.fs;
tt = S.tt;
frame_steps = S.frame_steps;
w_p = R.w_p; wxx_p = R.wxx_p; wyy_p = R.wyy_p; wxy_p = R.wxy_p;
Wf = R.Wf;
notes = [S.notes, R.notes];
rect = strcmp(dom.shape, 'rectangle');

% Truncate after a blow-up so the outputs stay finite
if R.diverged_step > 0
    nkeep = R.diverged_step;             % keep the last finite level
    tt = tt(1:nkeep); w_p = w_p(1:nkeep); wxx_p = wxx_p(1:nkeep); wyy_p = wyy_p(1:nkeep); wxy_p = wxy_p(1:nkeep);
    keepf = frame_steps < R.diverged_step;
    Wf = Wf(:, :, keepf); frame_steps = frame_steps(keepf);
end
nf = numel(frame_steps);
tf = frame_steps * dt;

% Stress at the saved frames, from the curvatures of each frame
VMf = zeros(ny, nx, nf); SPf = zeros(ny, nx, nf);
for kf = 1:nf
    if ~rect
        [wxx, wyy, wxy] = plate_curvature_domain(Wf(:, :, kf), dom);
    elseif S.free
        [wxx, wyy, wxy] = plate_curvature_free(Wf(:, :, kf), dx, dy, P.material.nu);
    else
        [wxx, wyy, wxy] = plate_curvature(plate_pad(Wf(:, :, kf), P.bc), dx, dy);
    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
if ~rect
    % Off the plate there is nothing to show
    out = ~dom.inside;
    for kf = 1:nf
        tmp = Wf(:, :, kf); tmp(out) = NaN; Wf(:, :, 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;
end

% Audio: sum the delay bins, then trim to the simulated duration
na = S.na;
p_audio = zeros(1, na);
for d = 1:S.Lis.nbins
    p_audio(d:end) = p_audio(d:end) + R.Pbins(d, 1:na - d + 1);
end
n_audio = min(na, ceil(numel(tt) * dt * fs));
p_audio = p_audio(1:n_audio);

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

result.method = 'fd';
result.grid = struct('x', S.x, 'y', S.y);
result.t = tt;
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, '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);
result.info = struct('nx', nx, 'ny', ny, 'dx', dx, 'dy', dy, 'dt', dt, 'dt_max', S.plan.dt_max, 'cfl', P.fd.cfl, ...
    'nsteps', S.nsteps, 'ncells', nx * ny, 'nunknowns', S.nunknowns, 'shape', dom.shape, 'cell_steps', S.plan.cell_steps, 'fs', fs, ...
    'stable', R.diverged_step == 0, 'diverged_step', R.diverged_step, ...
    'loop', R.backend, 'loop_s', R.loop_s, 'wall_s', toc(S.t_start), 'notes', {notes});
if ~isempty(S.fixed)
    result.info.support_x = S.X(S.fixed);
    result.info.support_y = S.Y(S.fixed);
end
end