function plan = plate_implicit_plan(params)
% PLATE_IMPLICIT_PLAN  Grid, time step and size of the implicit finite-difference scheme.
%
%   plan = plate_implicit_plan(params) returns nx, ny, dx, dy, dt, nsteps,
%   ncells, nunknowns and cell_steps for plate_implicit.  There is no stability
%   bound for theta >= 1/4; the step is set for accuracy as
%     dt = 1 / (fs * steps_per_sample)
%   and the trade-off is period elongation: for theta = 1/4 a mode of angular
%   frequency omega runs slow by a factor of about 1 + (omega dt)^2 / 12.  At one
%   step per audio sample a 10 kHz mode is detuned by about 17 percent; at two
%   steps per sample by about 4 percent.

P = plate_setup(params);
plan.nx = P.implicit.nx;
plan.ny = P.implicit.ny;
if plan.nx < 5 || plan.ny < 5
    error('plate_implicit_plan: the grid needs at least 5 nodes in each direction');
end
plan.dx = P.plate.Lx / (plan.nx - 1);
plan.dy = P.plate.Ly / (plan.ny - 1);
plan.dt = 1 / (P.output.fs * P.implicit.steps_per_sample);
plan.nsteps = ceil(P.implicit.duration / plan.dt);
plan.ncells = plan.nx * plan.ny;
dom = plate_domain(P.plate, plan.nx, plan.ny);
if any(strcmp(P.bc, {'free', 'pinned-center'}))
    plan.nunknowns = dom.nunknowns_free;        % every node of the plate is in the system; see plate_matrices
else
    plan.nunknowns = dom.nunknowns_dirichlet;   % the boundary ring carries w = 0
end
plan.cell_steps = plan.ncells * plan.nsteps;
plan.theta = P.implicit.theta;
end