function result = plate_fd(params)
% PLATE_FD  Explicit finite differences for the damped Kirchhoff-Love plate.
%
%   result = plate_fd(params) advances
%
%     w_tt + 2 sigma0 w_t - 2 sigma1 lap(w_t) + kappa^2 lap(lap(w)) = F(t) s(x,y) / mu
%
%   on a uniform nx x ny grid (nodes include the edges) with the scheme
%
%     (1 + sigma0 dt) w+ = 2 w - (1 - sigma0 dt) w- - dt^2 kappa^2 L(L w)
%                          + 2 sigma1 dt (L w - L w-) + dt^2 F s / mu
%
%   where w-, w, w+ are consecutive time levels and L is the five-point
%   Laplacian.  Applying L twice gives the 13-point biharmonic stencil.  The
%   scheme is second order in space and time (Bilbao, Numerical Sound
%   Synthesis, 2009, ch. 12).  Edge conditions enter through one layer of ghost
%   nodes (plate_pad): simply supported or clamped.
%
%   The step is dt = cfl * dt_max with dt_max from plate_fd_plan; the scheme
%   is unstable for cfl > 1 and the run stops as soon as the field diverges,
%   reporting the step at which it did (info.diverged_step).
%
%   Approximations.  Second-order truncation error in dx, dy, dt; the probe is
%   snapped to the nearest grid node; stress fields are evaluated only at the
%   saved frames, so the envelope is a sample of the true peak; the audio is
%   the Rayleigh integral box-averaged onto audio samples.
%
%   See plate_result_schema for the fields of result.

t_start = tic;
P = plate_setup(params);
plan = plate_fd_plan(params);
nx = plan.nx; ny = plan.ny; dx = plan.dx; dy = plan.dy; dt = plan.dt; nsteps = plan.nsteps;
x = linspace(0, P.plate.Lx, nx);
y = linspace(0, P.plate.Ly, ny);
[X, Y] = meshgrid(x, y);                 % w(i, j) lives at (x(j), y(i))
fs = P.output.fs;
notes = {};
if P.fd.cfl > 1
    notes{end + 1} = sprintf('cfl = %g exceeds the stability limit; expect the solution to blow up', P.fd.cfl);
end

s = plate_force_footprint(P.excitation, X, Y, dx * dy);
tt = (0:nsteps) * dt;
F = plate_force_time(P.excitation, tt);

% Probe: nearest interior node
[~, jp] = min(abs(x - P.probe.x)); jp = min(max(jp, 2), nx - 1);
[~, ip] = min(abs(y - P.probe.y)); ip = min(max(ip, 2), 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))));
nf = numel(frame_steps);
Wf = zeros(ny, nx, nf); VMf = zeros(ny, nx, nf); SPf = zeros(ny, nx, nf);
tf = zeros(1, nf);

% Listener: every grid node is a patch; the acceleration of a step is added into
% the audio sample that contains it, delayed by the node's travel time.
L = plate_listener(P.listener, X, Y, dx * dy, fs);
bins = L.delay(:) + 1;
wgt = L.weight(:);
na = ceil(nsteps * dt * fs) + L.nbins + 1;
Pbins = zeros(L.nbins, na);

% 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);
c_old = 1 - P.damping.sigma0 * dt;
c_new = 1 / (1 + P.damping.sigma0 * dt);
k2dt2 = P.kappa^2 * dt^2;
s1dt2 = 2 * P.damping.sigma1 * dt;
fdt2 = dt^2 / P.mu;
blowup = 1e3 * P.plate.h;                % deflection beyond which the run is declared unstable
diverged_step = 0;
kf = 1;

for nstep = 1:nsteps
    Lw = plate_laplacian(plate_pad(w, P.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;
    wn = c_new * (2 * w - c_old * wm - k2dt2 * B + s1dt2 * (Lw - Lw_m) + fdt2 * F(nstep) * s);
    wn(1, :) = 0; wn(end, :) = 0; wn(:, 1) = 0; wn(:, end) = 0;

    % 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(:), [L.nbins 1]) * (dt * fs);

    if kf <= nf && nstep == frame_steps(kf)
        [wxx, wyy, wxy] = plate_curvature(plate_pad(wn, P.bc), dx, dy);
        S = plate_stress(P, wxx, wyy, wxy);
        Wf(:, :, kf) = wn; VMf(:, :, kf) = S.vm; SPf(:, :, kf) = S.sp; tf(kf) = nstep * dt;
        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

% Truncate after a blow-up so the outputs stay finite
if diverged_step > 0
    nkeep = 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 < diverged_step;
    Wf = Wf(:, :, keepf); VMf = VMf(:, :, keepf); SPf = SPf(:, :, keepf); tf = tf(keepf);
end

% Audio: sum the delay bins, then trim to the simulated duration
p_audio = zeros(1, na);
for d = 1:L.nbins
    p_audio(d:end) = p_audio(d:end) + 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', x, 'y', y);
result.t = tt;
result.probe = struct('x', x(jp), 'y', y(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);
if isempty(tf)
    result.envelope = struct('vm', zeros(ny, nx), 'sp', zeros(ny, nx));
else
    result.envelope = struct('vm', max(VMf, [], 3), 'sp', max(abs(SPf), [], 3));
end
result.spectrum = struct('f', fw, 'w', Aw, 'vm', Avm);
result.audio = struct('fs', fs, 'p', p_audio, 't0', L.t0);
result.info = struct('nx', nx, 'ny', ny, 'dx', dx, 'dy', dy, 'dt', dt, 'dt_max', plan.dt_max, 'cfl', P.fd.cfl, ...
    'nsteps', nsteps, 'ncells', nx * ny, 'cell_steps', plan.cell_steps, 'fs', fs, ...
    'stable', diverged_step == 0, 'diverged_step', diverged_step, 'wall_s', toc(t_start), 'notes', {notes});
end