function M = plate_fem_modes(params)
% PLATE_FEM_MODES  Modes of the hung plate by the Q4 Mindlin-Reissner finite-element method.
%
%   M = plate_fem_modes(params) assembles the mesh (plate_fem_assemble) and
%   solves the generalized eigenproblem K v = omega^2 M v for its lowest
%   modes (plate_fem_eigs), enough of them to pass 0.47 fs.  Three modes lie
%   below 20 Hz: the two rigid rotations about axes through the hole, which
%   nothing restrains, and the pendulum mode of the thread's spring.  They
%   are kept, so that the collision can put momentum into the freely hanging
%   plate, and marked by their frequency; plate_synthesize leaves them out of
%   the sound.  The audible modes are those between 20 Hz and 0.47 fs, at
%   most params.render.modes of them.
%
%   Returns M with method 'fem', freqs (1 x nm, the low ones first), the
%   mass-normalized transverse shapes on the node grid, shapes ((ny + 1) x
%   (nx + 1) x nm, rows index y), the curvatures kxx, kyy, kxy of each shape
%   (from the rotations: w_xx = -d theta_x / dx and so on, element-centre
%   values averaged to the nodes), support_part (the shape at the hole),
%   volume_part (the trapezoidal integral of the shape over the plate),
%   grid_x, grid_y, the mesh sizes and the eigensolver's report.

t_start = tic;
P = plate_setup(params);
A = plate_fem_assemble(params, P);
nx = A.nx;
ny = A.ny;
nyn = ny + 1;
nxn = nx + 1;
notes = {};

% how many modes to ask for: the modal-density estimate plus a margin, then
% enlarge until the wanted set reaches past f_max or the cap
plan = plate_plan_common(params, P);
n_rigid = 3;
lam_max = (2 * pi * P.f_max)^2;
nwant = min(params.render.modes, plan.modes_expected + 12) + n_rigid;
nwant = max(nwant, n_rigid + 4);
tries = 0;
while true
    tries = tries + 1;
    nb = min(A.ndof, nwant + max(nwant, 8));
    [lam, V, einfo] = plate_fem_eigs(A.K, A.M, nwant, nb, 60, 1e-10);
    n_below = nnz(lam < lam_max);
    reached = lam(end) >= lam_max;
    capped = n_below - n_rigid >= params.render.modes;
    if reached || capped || nb >= A.ndof || tries >= 4
        break
    end
    nwant = ceil(1.5 * nwant) + 4;
end
freqs = sqrt(max(lam, 0))' / (2 * pi);
low = freqs < 20;
audible = (lam' > 0) & (freqs >= 20) & (freqs < P.f_max);
n_audible = nnz(audible);
count = min(params.render.modes, n_audible);
sel_a = find(audible);
chosen = [find(low), sel_a(1:count)];
freqs = freqs(chosen);
V = V(:, chosen);
nm = numel(freqs);
% explicit M-normalization, as the reference simulator does
for k = 1:nm
    V(:, k) = V(:, k) / sqrt(V(:, k)' * (A.M * V(:, k)));
end
if ~einfo.converged
    notes{end + 1} = sprintf('the eigensolver stopped after %d iterations without meeting its tolerance; largest residual %.1e', einfo.iterations, einfo.residual);
end
if ~reached && ~capped
    notes{end + 1} = sprintf('the eigensolver found %d modes but did not reach 0.47 fs; some audible modes may be missing', n_audible);
end
if count < n_audible
    notes{end + 1} = sprintf('%d modes lie below 0.47 fs; the first %d are kept (params.render.modes)', n_audible, count);
end
if nnz(low) ~= n_rigid
    notes{end + 1} = sprintf('%d modes below 20 Hz where 3 (two free rotations and the pendulum) were expected', nnz(low));
end

% shapes on the node grid: node ey (nx + 1) + ex -> array (ey + 1, ex + 1)
Wn = V(A.wdofs, :);
shapes = permute(reshape(Wn, nxn, nyn, nm), [2, 1, 3]);
RX = permute(reshape(V(2:3:end, :), nxn, nyn, nm), [2, 1, 3]);
RY = permute(reshape(V(3:3:end, :), nxn, nyn, nm), [2, 1, 3]);
dx = A.dx;
dy = A.dy;
% element-centre derivatives of the rotations, then the average of the elements around each node
drx_dx = ((RX(1:end - 1, 2:end, :) + RX(2:end, 2:end, :)) - (RX(1:end - 1, 1:end - 1, :) + RX(2:end, 1:end - 1, :))) / (2 * dx);
dry_dy = ((RY(2:end, 1:end - 1, :) + RY(2:end, 2:end, :)) - (RY(1:end - 1, 1:end - 1, :) + RY(1:end - 1, 2:end, :))) / (2 * dy);
drx_dy = ((RX(2:end, 1:end - 1, :) + RX(2:end, 2:end, :)) - (RX(1:end - 1, 1:end - 1, :) + RX(1:end - 1, 2:end, :))) / (2 * dy);
dry_dx = ((RY(1:end - 1, 2:end, :) + RY(2:end, 2:end, :)) - (RY(1:end - 1, 1:end - 1, :) + RY(2:end, 1:end - 1, :))) / (2 * dx);
kxx = -plate_fem_nodal(drx_dx, nyn, nxn);
kyy = -plate_fem_nodal(dry_dy, nyn, nxn);
kxy = -plate_fem_nodal(drx_dy + dry_dx, nyn, nxn) / 2;

% participations
support_part = (V' * A.support)';
area = ones(nyn, nxn) * dx * dy;
area([1, end], :) = area([1, end], :) * 0.5;
area(:, [1, end]) = area(:, [1, end]) * 0.5;
volume_part = (reshape(shapes, nyn * nxn, nm)' * area(:))';

M.method = 'fem';
M.freqs = freqs;
M.support_part = support_part;
M.volume_part = volume_part;
M.shapes = shapes;
M.kxx = kxx;
M.kyy = kyy;
M.kxy = kxy;
M.grid_x = A.grid_x;
M.grid_y = A.grid_y;
M.ndof = A.ndof;
M.nn = A.nn;
M.nx = nx;
M.ny = ny;
M.dx = dx;
M.dy = dy;
M.n_rigid = nnz(low);
M.n_found = numel(lam);
M.plate_mass = P.plate_mass;
M.eig_iterations = einfo.iterations;
M.eig_residual = einfo.residual;
M.eig_block = einfo.block;
M.size_note = sprintf('%d x %d Q4 elements, %d unknowns; %d Ritz vectors, %d subspace iterations', nx, ny, A.ndof, einfo.block, einfo.iterations);
M.wall_s = toc(t_start);
M.notes = notes;
end