function [wxx, wyy, wxy] = plate_curvature_free(w, dx, dy, nu)
% PLATE_CURVATURE_FREE  Second derivatives of the displacement of a plate with free edges.
%
%   [wxx, wyy, wxy] = plate_curvature_free(w, dx, dy, nu) returns the curvatures
%   at every node of the ny x nx field w (rows index y, columns index x).
%   Inside they are the usual central differences.  On the edges there is no
%   ghost layer to lean on, so the free-edge conditions themselves supply the
%   missing values:
%
%     no bending moment on an edge   w_nn = -nu w_tt   (n normal, t tangent)
%     no twisting moment at a corner w_xy = 0, hence also w_xx = w_yy = 0 there
%
%   The tangential curvature along an edge is a central difference along the
%   edge, and the twist on an edge is one-sided in the normal direction.  Both
%   direct stresses therefore vanish at a corner, and on an edge the normal
%   stress vanishes while the tangential stress is -E z w_tt.

[ny, nx] = size(w);
wxx = zeros(ny, nx);
wyy = zeros(ny, nx);
wxy = zeros(ny, nx);
wxx(:, 2:end - 1) = (w(:, 3:end) - 2 * w(:, 2:end - 1) + w(:, 1:end - 2)) / dx^2;
wyy(2:end - 1, :) = (w(3:end, :) - 2 * w(2:end - 1, :) + w(1:end - 2, :)) / dy^2;
% Edge columns and rows: the moment condition, using the tangential curvature just computed
wxx(2:end - 1, [1 end]) = -nu * wyy(2:end - 1, [1 end]);
wyy([1 end], 2:end - 1) = -nu * wxx([1 end], 2:end - 1);
% Twist: central inside, one-sided in the normal direction on the edges, zero at the corners
wxy(2:end - 1, 2:end - 1) = (w(3:end, 3:end) - w(3:end, 1:end - 2) - w(1:end - 2, 3:end) + w(1:end - 2, 1:end - 2)) / (4 * dx * dy);
wxy(2:end - 1, 1) = (w(3:end, 2) - w(3:end, 1) - w(1:end - 2, 2) + w(1:end - 2, 1)) / (2 * dx * dy);
wxy(2:end - 1, end) = (w(3:end, end) - w(3:end, end - 1) - w(1:end - 2, end) + w(1:end - 2, end - 1)) / (2 * dx * dy);
wxy(1, 2:end - 1) = (w(2, 3:end) - w(2, 1:end - 2) - w(1, 3:end) + w(1, 1:end - 2)) / (2 * dx * dy);
wxy(end, 2:end - 1) = (w(end, 3:end) - w(end, 1:end - 2) - w(end - 1, 3:end) + w(end - 1, 1:end - 2)) / (2 * dx * dy);
end