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]
%
%   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
%
%     w_tt + 2 sigma0 w_t - 2 sigma1 lap(w_t) + kappa^2 lap(lap(w)) = f(x,y,t) / mu
%
%   where f is the applied pressure [N/m^2], sigma0 [1/s] is a frequency-independent
%   loss rate and sigma1 [m^2/s] adds losses that grow with spatial frequency
%   (Bilbao, Numerical Sound Synthesis, 2009, ch. 12).  The surface stresses follow
%   from the curvatures (see 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 ~any(strcmp(params.bc, {'simply-supported', 'clamped', 'free', 'pinned-center'}))
    error('plate_setup: unknown boundary condition "%s" (expected simply-supported, clamped, free or pinned-center)', params.bc);
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