function result = plate_implicit(params)
% PLATE_IMPLICIT  Implicit finite differences (theta scheme) for the damped Kirchhoff-Love plate.
%
%   result = plate_implicit(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 the same grid and with the same spatial stencils as plate_fd, but the
%   stiffness term is averaged over three time levels,
%
%     (w+ - 2w + w-)/dt^2 + sigma0 (w+ - w-)/dt - sigma1 L (w+ - w-)/dt
%         + kappa^2 B (theta w+ + (1 - 2 theta) w + theta w-) = F s / mu,
%
%   with L and B the sparse Laplacian and biharmonic matrices on the interior
%   nodes (plate_matrices).  Multiplying by dt^2 and collecting terms gives one
%   linear system per step with a constant matrix,
%
%     A w+ = C w - A_old w- + dt^2 F s / mu
%     A     = (1 + sigma0 dt) I - sigma1 dt L + theta kappa^2 dt^2 B
%     C     = 2 I - (1 - 2 theta) kappa^2 dt^2 B
%     A_old = (1 - sigma0 dt) I + sigma1 dt L + theta kappa^2 dt^2 B.
%
%   A is symmetric positive definite, so it is factored once by sparse Cholesky
%   and each step costs two sparse triangular solves plus two sparse products.
%   (R' is formed once; forming it every step would cost more than the solves.)
%   For theta >= 1/4 the scheme is unconditionally stable (von Neumann analysis
%   of the amplification polynomial); theta = 1/4 is the trapezoidal rule.  The
%   step is therefore chosen for accuracy, dt = 1 / (fs steps_per_sample).
%
%   Approximations.  Second order in space and time.  Implicit averaging
%   lengthens the period of a mode by about (omega dt)^2 / 12, so the highest
%   modes are detuned rather than unstable.  Probe, frames, stress envelope and
%   audio follow the conventions of plate_fd (see plate_result_schema).

t_start = tic;
P = plate_setup(params);
plan = plate_implicit_plan(params);
nx = plan.nx; ny = plan.ny; dx = plan.dx; dy = plan.dy; dt = plan.dt; nsteps = plan.nsteps;
theta = P.implicit.theta;
x = linspace(0, P.plate.Lx, nx);
y = linspace(0, P.plate.Ly, ny);
[X, Y] = meshgrid(x, y);
fs = P.output.fs;
notes = {};
if theta < 0.25
    notes{end + 1} = sprintf('theta = %g is below 1/4, so the scheme is only conditionally stable', theta);
end

% ---- Matrices on the interior nodes, factored once
[L, B, mx, my] = plate_matrices(nx, ny, dx, dy, P.bc);
m = mx * my;
I = speye(m);
k2dt2 = P.kappa^2 * dt^2;
A = (1 + P.damping.sigma0 * dt) * I - P.damping.sigma1 * dt * L + theta * k2dt2 * B;
C = 2 * I - (1 - 2 * theta) * k2dt2 * B;
A_old = (1 - P.damping.sigma0 * dt) * I + P.damping.sigma1 * dt * L + theta * k2dt2 * B;
t_factor = tic;
R = chol(A);                             % A = R' * R
Rt = R';                                 % transposed once, not once per step
factor_s = toc(t_factor);

% ---- Load, probe, frames, listener (as in plate_fd)
s = plate_force_footprint(P.excitation, X, Y, dx * dy);
s_int = s(2:end - 1, 2:end - 1);
load = (dt^2 / P.mu) * s_int(:);
tt = (0:nsteps) * dt;
F = plate_force_time(P.excitation, tt);

[~, 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

frame_steps = unique(max(1, min(nsteps, round(linspace(0, P.implicit.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);

Lis = plate_listener(P.listener, X, Y, dx * dy, fs);
bins = Lis.delay(:) + 1;
wgt = Lis.weight(:);
na = ceil(nsteps * dt * fs) + Lis.nbins + 1;
Pbins = zeros(Lis.nbins, na);

% ---- Time stepping on the interior vector
w = zeros(m, 1); wm = zeros(m, 1);
wfull = 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);
blowup = 1e3 * P.plate.h;
diverged_step = 0;
kf = 1;

for nstep = 1:nsteps
    rhs = C * w - A_old * wm + F(nstep) * load;
    wn = R \ (Rt \ rhs);                 % forward then back substitution

    wfull(2:end - 1, 2:end - 1) = reshape(wn, my, mx);
    w_p(nstep + 1) = wfull(ip, jp);
    wxx_p(nstep + 1) = (wfull(ip, jp + 1) - 2 * wfull(ip, jp) + wfull(ip, jp - 1)) / dx^2;
    wyy_p(nstep + 1) = (wfull(ip + 1, jp) - 2 * wfull(ip, jp) + wfull(ip - 1, jp)) / dy^2;
    wxy_p(nstep + 1) = (wfull(ip + 1, jp + 1) - wfull(ip + 1, jp - 1) - wfull(ip - 1, jp + 1) + wfull(ip - 1, jp - 1)) / (4 * dx * dy);

    acc = zeros(ny, nx);
    acc(2:end - 1, 2:end - 1) = reshape((wn - 2 * w + wm) / dt^2, my, mx);
    ia = floor(nstep * dt * fs) + 1;
    Pbins(:, ia) = Pbins(:, ia) + accumarray(bins, wgt .* acc(:), [Lis.nbins 1]) * (dt * fs);

    if kf <= nf && nstep == frame_steps(kf)
        [wxx, wyy, wxy] = plate_curvature(plate_pad(wfull, P.bc), dx, dy);
        S = plate_stress(P, wxx, wyy, wxy);
        Wf(:, :, kf) = wfull; VMf(:, :, kf) = S.vm; SPf(:, :, kf) = S.sp; tf(kf) = nstep * dt;
        kf = kf + 1;
    end

    peak = max(abs(wn));
    if ~(peak < blowup)
        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;
end

if diverged_step > 0
    nkeep = diverged_step;
    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

p_audio = zeros(1, na);
for d = 1:Lis.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 = 'implicit';
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', Lis.t0);
result.info = struct('nx', nx, 'ny', ny, 'dx', dx, 'dy', dy, 'dt', dt, 'theta', theta, ...
    'steps_per_sample', P.implicit.steps_per_sample, 'nunknowns', m, 'nnz_factor', nnz(R), 'factor_s', factor_s, ...
    '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