function [R, Q, G] = plate_rigid_basis(dom, mnt, Mvec, unk)
% PLATE_RIGID_BASIS  Rigid-body coordinates of the plate about the mount's pivot.
%
%   [R, Q, G] = plate_rigid_basis(dom, mnt, Mvec, unk) returns the N x nr
%   basis R of the rigid-body motions the mount leaves free, evaluated at the
%   grid nodes and zero at the nodes that are not unknowns (unk false), the
%   projector Q = M R (R' M R)^-1 with M = diag(Mvec) the node masses, so that
%
%     c = Q' w              are the rigid-body coordinates of a field w and
%     w - R c               its flexural remainder, mass-orthogonal to R,
%
%   and the nr x nr rigid mass matrix G = R' M R.  For nr = 3 the columns are
%
%     1                       translation along the normal (heave, or the
%                             sway of a hanging plate's hole, in metres)
%     (P - pivot) . ed        rotation about the horizontal in-plane axis eh
%                             through the pivot, in radians (the swing of a
%                             hung plate, the tilt of one on pads)
%     (P - pivot) . eh        rotation about the in-plane axis ed through the
%                             pivot (the spin of a hung plate about its thread)
%
%   with mnt.pivot, mnt.ed (in-plane down) and mnt.eh (in-plane horizontal)
%   from plate_mount.  A mount that holds nodes still (the vise) has nr = 0.
%
%   The time stepping carries the two parts separately (plate_fd_loop): the
%   rigid coordinates as three scalar equations with the exact reduced
%   stiffness R' A R, the remainder on the grid.  Bending annihilates R, so
%   the split costs nothing in accuracy and it is what keeps the metre-scale
%   sway and the micrometre-scale flexure apart in single precision.

N = dom.nx * dom.ny;
nr = mnt.nr;
if nr == 0
    R = zeros(N, 0); Q = zeros(N, 0); G = zeros(0, 0);
    return
end
X = dom.X(:); Y = dom.Y(:);
rd = (X - mnt.pivot(1)) * mnt.ed(1) + (Y - mnt.pivot(2)) * mnt.ed(2);
rh = (X - mnt.pivot(1)) * mnt.eh(1) + (Y - mnt.pivot(2)) * mnt.eh(2);
R = [ones(N, 1), rd, rh];
R(~unk(:), :) = 0;
m = Mvec(:) .* double(unk(:));
G = R' * (m .* R);
Q = (m .* R) / G;
end