function [b, bxx, byy, bxy] = plate_rr_basis(x, y, order, width, height)
% PLATE_RR_BASIS  The tensor-Legendre basis and its curvatures at points of the plate.
%
%   [b, bxx, byy, bxy] = plate_rr_basis(x, y, order, width, height) evaluates
%   the (order + 1)^2 products P_i(2 x / W) P_j(2 y / H), i, j = 0..order, at
%   the paired points (x, y) (any common shape, in metres from the plate's
%   centre), as rows of b, together with the physical second derivatives
%   w_xx, w_yy and w_xy of each basis function.  The index k of the function
%   (i, j) runs with j fastest: k = i (order + 1) + j + 1.
%
%   Legendre products are a complete, well-conditioned basis for a free
%   rectangle: no boundary condition is built in, which is right for free
%   edges, and the natural conditions follow from minimizing the energy.

xi = 2 * x(:)' / width;
eta = 2 * y(:)' / height;
[px, dx, dxx] = plate_rr_legendre(order, xi);
[py, dy, dyy] = plate_rr_legendre(order, eta);
n = (order + 1)^2;
pts = numel(xi);
b = zeros(n, pts);
bxx = zeros(n, pts);
byy = zeros(n, pts);
bxy = zeros(n, pts);
sx = 2 / width;
sy = 2 / height;
k = 0;
for i = 1:order + 1
    for j = 1:order + 1
        k = k + 1;
        b(k, :) = px(i, :) .* py(j, :);
        bxx(k, :) = dxx(i, :) .* py(j, :) * sx^2;
        byy(k, :) = px(i, :) .* dyy(j, :) * sy^2;
        bxy(k, :) = dx(i, :) .* dy(j, :) * sx * sy;
    end
end
end