function M = plate_rr_modes(params)
% PLATE_RR_MODES  Modes of the hung plate by the Rayleigh-Ritz method.
%
%   M = plate_rr_modes(params) minimizes the Kirchhoff-Love bending energy
%
%     U = D/2 integral ( w_xx^2 + w_yy^2 + 2 nu w_xx w_yy + 2 (1 - nu) w_xy^2 ) dA
%
%   over the tensor-Legendre basis of order params.rr.basis_order
%   (plate_rr_basis), with the kinetic energy rho h / 2 integral w_t^2 dA, both
%   integrated by Gauss-Legendre quadrature of params.rr.quadrature_order
%   points per side.  The hole subtracts its area's mass at its position
%   (its effect on the stiffness is negligible at 0.8 mm) and the thread
%   adds the pendulum spring m g / L there.  The generalized eigenproblem
%   K a = omega^2 Mm a is solved densely after a Cholesky reduction; the
%   eigenvectors are mass-normalized.
%
%   Kept: modes above 20 Hz and below 0.47 fs, at most params.render.modes
%   of them.  The pendulum mode of the thread and the two rigid rotations
%   are dropped here (this method has no rigid-body content; the
%   finite-element method keeps them during the collision).
%
%   Returns M with method 'rr', freqs (1 x nm), vecs (n x nm) the Ritz
%   coefficients, order, support_part (the shape at the hole) and
%   volume_part (the integral of the shape over the plate) per mode, and the
%   shapes and curvatures on a display grid of params.rr.grid points per
%   side: shapes, kxx, kyy, kxy (ny x nx x nm; rows index y), grid_x, grid_y.

t_start = tic;
P = plate_setup(params);
pl = params.plate;
W = pl.width;
H = pl.height;
order = params.rr.basis_order;
nq = params.rr.quadrature_order;
mu = params.material.rho * pl.thickness;
nu = params.material.nu;
notes = {};

% quadrature over the plate
[qx, wx] = plate_gauss_legendre(nq);
[XX, YY] = ndgrid(qx, qx);
WW = wx(:) * wx(:)';
ww = WW(:)' * W * H / 4;
[b, bxx, byy, bxy] = plate_rr_basis(XX(:) * W / 2, YY(:) * H / 2, order, W, H);

Mm = mu * ((b .* ww) * b');
K = P.D * ((bxx .* ww) * bxx' + (byy .* ww) * byy' ...
    + nu * ((bxx .* ww) * byy' + (byy .* ww) * bxx') ...
    + 2 * (1 - nu) * ((bxy .* ww) * bxy'));

% the hole's mass and the thread's spring at the hole
ph = plate_rr_basis(P.hole_x, P.hole_y, order, W, H);
Mm = Mm - mu * P.hole_area * (ph * ph');
K = K + P.thread_k * (ph * ph');

% generalized symmetric eigenproblem by Cholesky reduction
Mm = (Mm + Mm') / 2;
K = (K + K') / 2;
R = chol(Mm);
Cm = R' \ K / R;
Cm = (Cm + Cm') / 2;
[V, Dm] = eig(Cm);
[vals, order_] = sort(diag(Dm)');
V = V(:, order_);
vecs = R \ V;

keep = vals > 1e-5;
vals = vals(keep);
vecs = vecs(:, keep);
freqs = sqrt(max(vals, 0)) / (2 * pi);
audible = (freqs > 20) & (freqs < P.f_max);
n_audible = nnz(audible);
count = min(params.render.modes, n_audible);
sel = find(audible);
sel = sel(1:count);
freqs = freqs(sel);
vecs = vecs(:, sel);
nm = numel(freqs);
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

% participation of the hole and of the plate's volume (surface quadrature of order max(36, nq))
support_part = (vecs' * ph)';
ns = max(36, nq);
[qs, ws] = plate_gauss_legendre(ns);
[XS, YS] = ndgrid(qs, qs);
WS = ws(:) * ws(:)';
area_w = WS(:) * W * H / 4;
bs = plate_rr_basis(XS(:) * W / 2, YS(:) * H / 2, order, W, H);
volume_part = ((vecs' * bs) * area_w)';

% shapes and curvatures on the display grid
ng = params.rr.grid;
grid_x = linspace(-W / 2, W / 2, ng);
grid_y = linspace(-H / 2, H / 2, ng);
[YG, XG] = ndgrid(grid_y, grid_x);
[bg, bgxx, bgyy, bgxy] = plate_rr_basis(XG(:), YG(:), order, W, H);
shapes = reshape(bg' * vecs, ng, ng, nm);
kxx = reshape(bgxx' * vecs, ng, ng, nm);
kyy = reshape(bgyy' * vecs, ng, ng, nm);
kxy = reshape(bgxy' * vecs, ng, ng, nm);

n = (order + 1)^2;
M.method = 'rr';
M.freqs = freqs;
M.vecs = vecs;
M.order = order;
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 = grid_x;
M.grid_y = grid_y;
M.ndof = n;
M.n_rigid = 0;
M.n_found = numel(vals);
M.plate_mass = P.plate_mass;
M.size_note = sprintf('%d Legendre coefficients (order %d), %d x %d quadrature, dense eigensolve', n, order, nq, nq);
M.wall_s = toc(t_start);
M.notes = notes;
end