function [wxx, wyy, wxy] = plate_curvature(wp, dx, dy)
% PLATE_CURVATURE  Second derivatives of a padded displacement field.
%
%   [wxx, wyy, wxy] = plate_curvature(wp, dx, dy) returns second-order central
%   differences at every node of the original ny x nx grid, given the padded
%   (ny+2) x (nx+2) field wp from plate_pad.  Rows index y, columns index x.

c = wp(2:end - 1, 2:end - 1);
wxx = (wp(2:end - 1, 3:end) - 2 * c + wp(2:end - 1, 1:end - 2)) / dx^2;
wyy = (wp(3:end, 2:end - 1) - 2 * c + wp(1:end - 2, 2:end - 1)) / dy^2;
wxy = (wp(3:end, 3:end) - wp(3:end, 1:end - 2) - wp(1:end - 2, 3:end) + wp(1:end - 2, 1:end - 2)) / (4 * dx * dy);
end