function [R, Q] = plate_rigid_modes(X, Y, mvec, pivot)
% PLATE_RIGID_MODES  Basis of the rigid-body motions of a free plate and its M-orthogonal projector.
%
%   [R, Q] = plate_rigid_modes(X, Y, mvec) takes the node coordinates X, Y and
%   the diagonal mvec of the mass matrix (the trapezoidal weights of
%   plate_matrices) and returns R = [1, x - xc, y - yc] evaluated at the nodes
%   (n x 3) and Q = M R (R' M R)^-1 (n x 3), so that
%
%     w - R * (Q' * w)
%
%   is the part of w orthogonal to the rigid-body motions in the mass inner
%   product: what remains when the mean displacement and the two mean tilts are
%   taken out.
%
%   [R, Q] = plate_rigid_modes(X, Y, mvec, pivot), with pivot = [x0 y0] the
%   coordinates of a pinned node, keeps only the two rotations about that point,
%   R = [x - x0, y - y0].  Both columns vanish at the pin, so the projection
%   leaves w = 0 there untouched.
%
%   Why remove them.  A completely free plate has three zero-frequency modes,
%   and a strike gives it net momentum and angular momentum: with only the loss
%   sigma0 to stop it, the plate drifts to a displacement of order
%   (impulse / mass) / (2 sigma0), millimetres against the micrometres of the
%   bending motion.  A real free-free measurement rests the plate on foam or
%   hangs it from threads, which absorb exactly this motion and nothing else,
%   so the solvers project it out after every step.  The bending stiffness K
%   annihilates R exactly, so this changes nothing about the flexural dynamics;
%   the sigma1 loss operator does not, and the projection also keeps that term
%   from feeding a rigid tilt back into flexure.

x = X(:);
y = Y(:);
m = mvec(:);
if nargin >= 4
    R = [x - pivot(1), y - pivot(2)];
else
    xc = sum(m .* x) / sum(m);
    yc = sum(m .* y) / sum(m);
    R = [ones(numel(x), 1), x - xc, y - yc];
end
G = R' * (m .* R);                        % 3 x 3 Gram matrix in the mass inner product
Q = (m .* R) / G;
end