function S = plate_stress(P, wxx, wyy, wxy)
% PLATE_STRESS  Surface stresses of a plate from its curvatures.
%
%   S = plate_stress(P, wxx, wyy, wxy) takes the second derivatives of the
%   transverse 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]
%
%   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)).
%
%   For the Mindlin-Reissner method the curvatures are the gradients of the
%   rotations, -d(theta_x)/dx and so on; the same formula then gives the
%   bending stress, without the transverse-shear stress, which is small at
%   these thicknesses.

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