function mdl = plate_model(params)
% PLATE_MODEL  The discrete plate in its mount: operators, statics and time step.
%
%   mdl = plate_model(params) assembles everything about the problem that
%   does not depend on the excitation: the grid (plate_domain), the mount
%   (plate_mount), the bending and Laplacian operators (plate_bending), the
%   self-weight prestress (plate_prestress), the rigid-body basis
%   (plate_rigid_basis), the static equilibrium (plate_static), the reduced
%   rigid-body system and the stable time step.  plate_fd_plan reports its
%   numbers before a run and plate_fd_setup builds the run on it.
%
%   The semi-discrete equation on the unknown nodes is
%
%     M w_tt + 2 sigma0 M w_t - 2 sigma1 mu Le w_t + A0 w + f_contact(w, w_t) = F(t) s + M g_n
%     A0 = Kb + Kg + Ks
%
%   with M the node masses, Kb the bending stiffness, Kg the prestress
%   operator, Ks the linear springs of the mount (threads, fingers), Le the
%   Laplacian energy operator, and f_contact the one-sided pads and felts
%   with their dashpots.  Writing w = w0 + R c + phi, with w0 the static
%   deflection, R the rigid-body basis and c its coordinates, and phi the
%   flexural remainder (mass-orthogonal to R), the equation splits into
%
%     G (c_tt + 2 sigma0 c_t) = R' f - Krr c - AR' phi,          G = R' M R,
%     phi_tt + 2 sigma0 phi_t = M^-1 (f - A0 phi - AR c) - R (c_tt + 2 sigma0 c_t),
%
%   where f collects the applied and contact forces about rest and the sigma1
%   loss 2 sigma1 mu Le w_t of the total velocity (so that it enters both
%   parts and stays dissipative), AR = A0 R is the force of a unit rigid
%   displacement (with the thread's torsion added as a force on the spin) and
%   Krr = R' AR the reduced stiffness of the rigid motions (the pendulum
%   stiffness, the thread springs, the torsion).  Since Kb R = 0, AR is
%   computed from Kg + Ks alone; both are exact to rounding.  The sum
%   R c + phi then obeys exactly the leapfrog of the whole system, and the
%   rigid content of phi, Q' phi, stays zero: every force acts on the whole
%   plate, and the rigid equations are its projection.  (Leaving AR c out of
%   the flexural equation, an early version of this code, made the two parts
%   pump each other and the swing of a hung plate grow without bound.)
%
%   Stable time step.  For the bending operator alone the von Neumann bound of
%   the explicit scheme is dt <= 2 / (lambda (sigma1 + sqrt(kappa^2 + sigma1^2)))
%   with lambda = 4/dx^2 + 4/dy^2 the largest eigenvalue of the five-point
%   Laplacian.  The extra stiffness of the prestress, the springs and the
%   engaged contacts adds a rate Omega^2 bounded by Gershgorin's theorem on
%   M^-1 (Kg + Ks + K_contact), and the bound becomes
%
%     dt <= ( -2 sigma1 lambda + sqrt(4 sigma1^2 lambda^2 + 4 (kappa^2 lambda^2 + Omega^2)) )
%           / (kappa^2 lambda^2 + Omega^2),
%
%   which reduces to the first form when Omega = 0.  The dashpots are treated
%   explicitly and need dt < 2 M / c; the plan keeps dt below half of that.
%   The largest eigenvalue of the energy-method bending operator lies just
%   below the interior bound (about 0.97 to 0.99 of it), so the step is safe
%   and slightly conservative.
%
%   Fields: P, dom, mnt, unk, Mvec, Kb, Le, Lm, Kg, Ks, A0, R, Q, G, Ginv,
%   Krr, AR, ps (prestress solution), st (statics), rigid_hz (frequencies of
%   the rigid-body modes, ascending, NaN for an unstable one) and
%   rigid_mode_names (the coordinate that dominates each), lam_max, Om2, dt_max, dt,
%   nsteps, ncells, nunknowns, notes.

P = plate_setup(params);
nx = P.fd.nx; ny = P.fd.ny;
if nx < 5 || ny < 5
    error('plate_model: the grid needs at least 5 nodes in each direction');
end
dom = plate_domain(P.plate, nx, ny);
dx = dom.dx; dy = dom.dy;
N = nx * ny;
mnt = plate_mount(P, dom);
unk = dom.inside & ~mnt.held;
a = dom.m * dx * dy;
Mvec = P.mu * a(:);
Mvec(~unk(:)) = 1;                                    % dummy mass at nodes outside the system
[Kb, Le, ops] = plate_bending(dom, mnt.held, P.material.nu, P.D);
[Kg, ps] = plate_prestress(P, dom, mnt, ops);
ksv = mnt.ks(:) .* double(unk(:));
Ks = spdiags(ksv, 0, N, N);
A0 = Kb + Kg + Ks;
Lm = spdiags(1 ./ max(a(:), eps) .* double(unk(:)), 0, N, N) * Le;
[R, Q, G] = plate_rigid_basis(dom, mnt, Mvec, unk);
st = plate_static(P, dom, mnt, A0, Mvec, unk);
nr = mnt.nr;
notes = [mnt.notes, st.notes];

% Reduced rigid-body system: bending annihilates R, so only Kg and Ks act on it
if nr > 0
    AR = (Kg + Ks) * R;
    if mnt.ktor > 0
        % the thread's torsion as a force on the plate: -ktor c3 M r_h / (r_h' M r_h),
        % which changes only the spin coordinate, so that Krr = R' AR stays exact
        AR(:, 3) = AR(:, 3) + mnt.ktor * (Mvec .* R(:, 3)) / G(3, 3);
    end
    Krr = R' * AR;
    Krr = (Krr + Krr') / 2;
    Ginv = inv(G);
    % modes of the reduced system in the mass-normalized coordinates, so that
    % the dominant coordinate of each is meaningful across metres and radians
    Rg = chol(G);                                    % G = Rg' Rg
    Bm = (Rg' \ Krr) / Rg;
    [V, Dm] = eig((Bm + Bm') / 2);
    lam = real(diag(Dm));
    [lam, order] = sort(lam);
    V = real(V(:, order));
    rigid_hz = sqrt(max(lam, 0)) / (2 * pi);
    rigid_hz(lam < -1e-9 * max(abs(lam))) = NaN;
    rigid_mode_names = cell(1, nr);
    for k = 1:nr
        [~, kd] = max(abs(V(:, k)));
        rigid_mode_names{k} = mnt.rigid_names{kd};
    end
    if any(isnan(rigid_hz))
        notes{end + 1} = 'a rigid-body motion has negative stiffness: the plate is in an unstable position and the run will diverge';
    end
else
    AR = zeros(N, 0); Krr = zeros(0, 0); Ginv = zeros(0, 0); rigid_hz = zeros(0, 1); rigid_mode_names = {};
end

% Stable time step
lam_max = 4 / dx^2 + 4 / dy^2;
kc = mnt.kb(:) .* double(st.actb) + mnt.kt(:) .* double(st.actt);
Kx = abs(Kg) + spdiags(abs(ksv) + kc, 0, N, N);
row = sum(Kx, 2);
Om2 = max(row(unk(:)) ./ Mvec(unk(:)));
if isempty(Om2)
    Om2 = 0;
end
s1 = P.damping.sigma1;
om2b = P.kappa^2 * lam_max^2 + Om2;
dt_max = (-2 * s1 * lam_max + sqrt(4 * s1^2 * lam_max^2 + 4 * om2b)) / om2b;
cv = (mnt.cb(:) + mnt.ct(:) + mnt.cs(:)) ./ Mvec;
cmax = max(cv(unk(:)));
if cmax > 0 && dt_max > 0.5 / cmax
    dt_max = 0.5 / cmax;
    notes{end + 1} = 'the support damping, not the bending, limits the time step';
end
dt = P.fd.cfl * dt_max;
nsteps = ceil(P.fd.duration / dt);

mdl.P = P; mdl.dom = dom; mdl.mnt = mnt; mdl.unk = unk;
mdl.nx = nx; mdl.ny = ny; mdl.dx = dx; mdl.dy = dy; mdl.N = N;
mdl.a = a; mdl.Mvec = Mvec;
mdl.Kb = Kb; mdl.Le = Le; mdl.Lm = Lm; mdl.Kg = Kg; mdl.Ks = Ks; mdl.A0 = A0;
mdl.R = R; mdl.Q = Q; mdl.G = G; mdl.Ginv = Ginv; mdl.Krr = Krr; mdl.AR = AR; mdl.nr = nr;
mdl.ps = ps; mdl.st = st; mdl.ops = ops;
mdl.rigid_hz = rigid_hz; mdl.rigid_mode_names = rigid_mode_names;
mdl.lam_max = lam_max; mdl.Om2 = Om2;
mdl.dt_max = dt_max; mdl.dt = dt; mdl.nsteps = nsteps;
mdl.ncells = N; mdl.nunknowns = nnz(unk);
mdl.notes = notes;
end