function [idx, wts, ex, ey] = plate_fem_interp(x, y, params)
% PLATE_FEM_INTERP  Where a point of the plate falls in the mesh, and its shape weights.
%
%   [idx, wts, ex, ey] = plate_fem_interp(x, y, params) finds the element
%   (ex, ey, zero-based) containing the point (x, y) (metres from the
%   plate's centre; a point on the far edge belongs to the last element) and
%   returns the bilinear weights wts (1 x 4) of its four nodes, with idx
%   (1 x 4) the linear index of each node in the node grid stored as an
%   (ny + 1) x (nx + 1) array with rows indexing y.  A field sampled at the
%   nodes is interpolated as wts * field(idx).

nx = params.fem.mesh_x;
ny = params.fem.mesh_y;
W = params.plate.width;
H = params.plate.height;
dx = W / nx;
dy = H / ny;
ux = min(max((x + W / 2) / dx, 0), nx - 1e-12);
uy = min(max((y + H / 2) / dy, 0), ny - 1e-12);
ex = min(floor(ux), nx - 1);
ey = min(floor(uy), ny - 1);
xi = 2 * (ux - ex) - 1;
eta = 2 * (uy - ey) - 1;
wts = plate_fem_shape(xi, eta);
nyn = ny + 1;
% nodes anticlockwise from (ex, ey): rows are y, so the linear index is (ey + 1) + ex * nyn
idx = [(ey + 1) + ex * nyn, (ey + 1) + (ex + 1) * nyn, (ey + 2) + (ex + 1) * nyn, (ey + 2) + ex * nyn];
end