function P = plate_setup(params)
% PLATE_SETUP  Derived quantities for a Kirchhoff-Love plate problem.
%
%   P = plate_setup(params) returns a copy of params with the derived fields
%     P.D      bending stiffness        D = E h^3 / (12 (1 - nu^2))   [N m]
%     P.mu     mass per unit area       mu = rho h                    [kg/m^2]
%     P.kappa  stiffness parameter      kappa = sqrt(D / mu)          [m^2/s]
%     P.Kb     bending-stress factor at the surface z = h/2          [N/m]
%     P.Ks     shear-stress factor at the surface z = h/2            [N/m]
%     P.g      gravitational acceleration [m/s^2], params.gravity.g (0 switches
%              gravity off)
%
%   The plate is params.plate.shape, 'rectangle' or 'ellipse', filling the
%   box [0, Lx] x [0, Ly]: for the ellipse Lx and Ly are the diameters along
%   x and y, and Lx = Ly is a disk.  A missing shape means a rectangle.
%
%   The transverse displacement w(x,y,t) of the plate obeys
%
%     mu w_tt + 2 mu sigma0 w_t - 2 mu sigma1 lap(w_t) + D lap(lap(w))
%         - div(N grad w) = f(x,y,t) + mu g_n + support forces
%
%   where f is the applied pressure [N/m^2], sigma0 [1/s] is a frequency-
%   independent loss rate, sigma1 [m^2/s] adds losses that grow with spatial
%   frequency (Bilbao, Numerical Sound Synthesis, 2009, ch. 12), N is the
%   in-plane stress resultant of the plate's own weight (plate_prestress) and
%   g_n the component of gravity along the plate normal.  How the plate is
%   held (params.mount) decides the support forces and the direction of
%   gravity in the plate's own coordinates: see plate_mount.  The surface
%   stresses follow from the curvatures (plate_stress).

P = params;
pl = params.plate;
mat = params.material;

if ~isfield(pl, 'shape')
    pl.shape = 'rectangle';
    P.plate.shape = 'rectangle';
end
if ~any(strcmp(pl.shape, {'rectangle', 'ellipse'}))
    error('plate_setup: unknown plate shape "%s" (expected rectangle or ellipse)', pl.shape);
end
if pl.Lx <= 0 || pl.Ly <= 0 || pl.h <= 0
    error('plate_setup: plate dimensions must be positive');
end
if pl.h > 0.1 * min(pl.Lx, pl.Ly)
    warning('plate_setup: thickness is more than a tenth of the shorter side; Kirchhoff-Love thin-plate theory is not reliable here');
end
if ~isfield(params, 'mount') || ~isfield(params.mount, 'type')
    error('plate_setup: params.mount.type is required (thread, threads, pads, cymbal, vise or pinched)');
end
if ~any(strcmp(params.mount.type, {'thread', 'threads', 'pads', 'cymbal', 'vise', 'pinched'}))
    error('plate_setup: unknown mount type "%s"', params.mount.type);
end
if isfield(params, 'gravity') && isfield(params.gravity, 'g')
    P.g = params.gravity.g;
else
    P.g = 9.81;
end
if P.g < 0
    error('plate_setup: g cannot be negative; turn the plate over instead');
end

P.D = mat.E * pl.h^3 / (12 * (1 - mat.nu^2));
P.mu = mat.rho * pl.h;
P.kappa = sqrt(P.D / P.mu);
P.Kb = mat.E * pl.h / (2 * (1 - mat.nu^2));
P.Ks = mat.E * pl.h / (2 * (1 + mat.nu));
end