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 of plate_matrices
%   on the plate's footprint (plate_domain), which also supplies the diagonal
%   mass matrix M: the identity on the interior nodes for simply supported or
%   clamped edges, the trapezoidal node weights on every node for free edges.
%   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 M s / mu
%     A     = (1 + sigma0 dt) M - sigma1 dt L + theta kappa^2 dt^2 B
%     C     = 2 M - (1 - 2 theta) kappa^2 dt^2 B
%     A_old = (1 - sigma0 dt) M + sigma1 dt L + theta kappa^2 dt^2 B.
%
%   For free edges the rigid-body motions are projected out of w+ after every
%   solve (plate_rigid_modes); for the centre-pinned plate the pinned node is
%   zeroed first.
%
%   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); on an
%   ellipse the edge is a staircase and first order, and the nodes outside the
%   plate are NaN in the fields.

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;
dom = plate_domain(P.plate, nx, ny);
x = dom.x; y = dom.y; X = dom.X; Y = dom.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
rect = strcmp(dom.shape, 'rectangle');

% ---- Matrices on the unknowns, factored once
free = any(strcmp(P.bc, {'free', 'pinned-center'}));
[L, B, M, idx, fixed] = plate_matrices(dom, P.bc, P.material.nu);
m = numel(idx);
mvec = full(diag(M));
unknown = false(ny, nx); unknown(idx) = true;
k2dt2 = P.kappa^2 * dt^2;
A = (1 + P.damping.sigma0 * dt) * M - P.damping.sigma1 * dt * L + theta * k2dt2 * B;
C = 2 * M - (1 - 2 * theta) * k2dt2 * B;
A_old = (1 - P.damping.sigma0 * dt) * M + 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, unknown);
load = (dt^2 / P.mu) * (mvec .* s(idx));   % M s on the unknowns
if free && isempty(fixed)
    [Rb, Qb] = plate_rigid_modes(X(idx), Y(idx), mvec);
    notes{end + 1} = 'free edges: the rigid-body translation and tilts are removed from the displacement after every step';
elseif free
    [Rb, Qb] = plate_rigid_modes(X(idx), Y(idx), mvec, [X(idx(fixed)) Y(idx(fixed))]);
    notes{end + 1} = sprintf('pinned at the node (%.4g, %.4g) m nearest the centre; the two rigid rotations about it are removed after every step', X(idx(fixed)), Y(idx(fixed)));
end
if ~rect
    notes{end + 1} = sprintf('%s on a staircase of %d grid nodes; the edge is first order in the grid spacing', dom.shape, nnz(dom.inside));
end
tt = (0:nsteps) * dt;
F = plate_force_time(P.excitation, tt);

% Probe: the nearest node whose eight neighbours belong to the plate
d2 = (X - P.probe.x).^2 + (Y - P.probe.y).^2;
d2(~dom.full3) = Inf;
[~, kp] = min(d2(:));
ip = mod(kp - 1, ny) + 1;
jp = floor((kp - 1) / 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 vector of unknowns
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
    if free
        if ~isempty(fixed)
            wn(fixed) = 0;
        end
        wn = wn - Rb * (Qb' * wn);
    end

    wfull(idx) = wn;
    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(idx) = (wn - 2 * w + wm) / dt^2;
    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)
        if ~rect
            [wxx, wyy, wxy] = plate_curvature_domain(wfull, dom);
        elseif free
            [wxx, wyy, wxy] = plate_curvature_free(wfull, dx, dy, P.material.nu);
        else
            [wxx, wyy, wxy] = plate_curvature(plate_pad(wfull, P.bc), dx, dy);
        end
        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
nf = numel(tf);
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
    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

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);
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', 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, 'shape', dom.shape, '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});
if ~isempty(fixed)
    result.info.support_x = X(idx(fixed));
    result.info.support_y = Y(idx(fixed));
end
end