function S = plate_stress(P, wxx, wyy, wxy)
% PLATE_STRESS  Surface stresses of a Kirchhoff-Love plate from its curvatures.
%
%   S = plate_stress(P, wxx, wyy, wxy) takes the second derivatives of the
%   displacement (arrays of any common size) and returns, at the surface z = h/2,
%     S.sxx, S.syy, S.txy   in-plane stress components                  [Pa]
%     S.vm                  von Mises equivalent stress                 [Pa]
%     S.sp                  principal stress of largest magnitude, signed
%                           (tension positive)                          [Pa]
%
%   In thin-plate theory the in-plane strains vary linearly through the
%   thickness, eps_xx = -z w_xx, eps_yy = -z w_yy, gamma_xy = -2 z w_xy, and
%   plane stress gives
%     sigma_xx = -E z / (1 - nu^2) (w_xx + nu w_yy)
%     sigma_yy = -E z / (1 - nu^2) (w_yy + nu w_xx)
%     tau_xy   = -E z / (1 + nu)   w_xy
%   Stresses are largest at the surfaces z = +-h/2 and vanish at the mid-plane;
%   the two surfaces carry equal and opposite stress.  P.Kb and P.Ks are the
%   factors E h / (2 (1 - nu^2)) and E h / (2 (1 + nu)) (see plate_setup).
%
%   Because stress is a second derivative of the primary unknown, any error in
%   the computed displacement is amplified here.

nu = P.material.nu;
S.sxx = -P.Kb * (wxx + nu * wyy);
S.syy = -P.Kb * (wyy + nu * wxx);
S.txy = -P.Ks * wxy;

S.vm = sqrt(S.sxx.^2 - S.sxx .* S.syy + S.syy.^2 + 3 * S.txy.^2);

c = (S.sxx + S.syy) / 2;
r = sqrt(((S.sxx - S.syy) / 2).^2 + S.txy.^2);
s1 = c + r;
s2 = c - r;
S.sp = s1;
pick = abs(s2) > abs(s1);
S.sp(pick) = s2(pick);
end