function result = plate_modal(params)
% PLATE_MODAL  Modal superposition for a simply supported rectangular plate.
%
%   result = plate_modal(params) solves the damped Kirchhoff-Love equation
%
%     w_tt + 2 sigma0 w_t - 2 sigma1 lap(w_t) + kappa^2 lap(lap(w)) = F(t) s(x,y) / mu
%
%   on 0 < x < a, 0 < y < b with simply supported edges (w = lap(w) = 0).  For
%   these edges the eigenfunctions are known exactly:
%
%     phi_mn(x,y) = sin(m pi x / a) sin(n pi y / b)
%     omega_mn    = kappa lambda_mn,   lambda_mn = (m pi / a)^2 + (n pi / b)^2
%
%   Expanding w = sum_mn q_mn(t) phi_mn(x,y) decouples the plate into
%   independent damped oscillators, one per mode,
%
%     q_mn'' + 2 sigma_mn q_mn' + omega_mn^2 q_mn = G_mn F(t) / M_mn
%
%   with modal mass M_mn = mu a b / 4, modal force factor G_mn = integral of
%   s(x,y) phi_mn over the footprint, and damping rate sigma_mn = sigma0 +
%   sigma1 lambda_mn.  Each oscillator is integrated exactly for
%   piecewise-linear forcing (plate_sdof).  Stresses follow from the analytic
%   second derivatives of the modes, and the radiated pressure from the modal
%   accelerations through the Rayleigh integral (plate_listener).
%
%   The rectangle is the only shape with closed-form modes for these edges;
%   an ellipse is refused, and the finite-difference methods handle it.
%
%   Approximations.  The modal sum is truncated at params.modal.fmax and
%   params.modal.max_modes; the force is sampled at params.output.fs; the
%   footprint integral and the Rayleigh integral use quadrature grids.  Stress
%   converges more slowly than displacement as modes are added, and the
%   damping model is a choice, not a measurement.
%
%   See plate_result_schema for the fields of result.

t_start = tic;
P = plate_setup(params);
if ~strcmp(P.bc, 'simply-supported')
    error('plate_modal: closed-form modes exist only for simply supported edges (got "%s"); use the finite-difference method for other edges', P.bc);
end
if ~strcmp(P.plate.shape, 'rectangle')
    error('plate_modal: closed-form modes exist only for the rectangle (got "%s"); use a finite-difference method for the ellipse', P.plate.shape);
end
a = P.plate.Lx;
b = P.plate.Ly;
fs = P.output.fs;
dt = 1 / fs;
nt = round(P.output.duration * fs) + 1;      % samples spanning [0, duration]
t = (0:nt - 1) * dt;
notes = {};

% ---- Mode selection: every mode below fmax, lowest first, at most max_modes
omega_max = 2 * pi * P.modal.fmax;
m_hi = max(1, ceil(a / pi * sqrt(omega_max / P.kappa)));
n_hi = max(1, ceil(b / pi * sqrt(omega_max / P.kappa)));
[MM, NN] = meshgrid(1:m_hi, 1:n_hi);
lam_all = (MM * pi / a).^2 + (NN * pi / b).^2;
keep = P.kappa * lam_all <= omega_max;
m = MM(keep); n = NN(keep); lam = lam_all(keep);
[lam, order] = sort(lam);
m = m(order); n = n(order);
truncated = numel(lam) > P.modal.max_modes;
if truncated
    notes{end + 1} = sprintf('%d modes lie below %g Hz but only the lowest %d were kept', numel(lam), P.modal.fmax, P.modal.max_modes);
    m = m(1:P.modal.max_modes); n = n(1:P.modal.max_modes); lam = lam(1:P.modal.max_modes);
end
m = m(:); n = n(:); lam = lam(:);
nm = numel(lam);
omega = P.kappa * lam;
sigma = P.damping.sigma0 + P.damping.sigma1 * lam;
over = sigma >= omega;
if any(over)
    notes{end + 1} = sprintf('%d modes would be overdamped by this damping model; their damping was capped just below critical', sum(over));
    sigma(over) = 0.99 * omega(over);
end
M = P.mu * a * b / 4;                     % modal mass (same for every mode)
kx = m * pi / a;                          % modal wavenumbers
ky = n * pi / b;

% ---- Applied force: total force F(t) and its projection G_mn onto each mode.
%      G_mn = integral s phi_mn dA over the footprint, by midpoint quadrature on a fine local grid.
F = plate_force_time(P.excitation, t);
exc = P.excitation;
nq = 64;
qx = linspace(max(0, exc.x - exc.radius), min(a, exc.x + exc.radius), nq);
qy = linspace(max(0, exc.y - exc.radius), min(b, exc.y + exc.radius), nq);
[QX, QY] = meshgrid(qx, qy);
dAq = (qx(2) - qx(1)) * (qy(2) - qy(1));
s = plate_force_footprint(exc, QX, QY, dAq);
PhiQ = sin(QX(:) * kx') .* sin(QY(:) * ky');             % (nq*nq) x nm
G = (s(:)' * PhiQ)' * dAq;                                % nm x 1

% ---- Mode shapes and their second derivatives at the probe
xp = P.probe.x; yp = P.probe.y;
phi_p = sin(kx * xp) .* sin(ky * yp);
phixx_p = -kx.^2 .* phi_p;
phiyy_p = -ky.^2 .* phi_p;
phixy_p = kx .* ky .* cos(kx * xp) .* cos(ky * yp);

% ---- Display grid for animation frames
ngx = P.output.grid(1); ngy = P.output.grid(2);
gx = linspace(0, a, ngx); gy = linspace(0, b, ngy);
[GX, GY] = meshgrid(gx, gy);
Phi = sin(GX(:) * kx') .* sin(GY(:) * ky');               % (ngy*ngx) x nm
Phixx = -Phi .* repmat(kx'.^2, ngx * ngy, 1);
Phiyy = -Phi .* repmat(ky'.^2, ngx * ngy, 1);
Phixy = (cos(GX(:) * kx') .* cos(GY(:) * ky')) .* repmat((kx .* ky)', ngx * ngy, 1);
% Frames at the samples nearest to nframes equally spaced times in [0, duration]
frame_idx = unique(min(nt, round(linspace(0, P.output.duration, P.output.nframes) * fs) + 1));
nf = numel(frame_idx);
Qf = zeros(nm, nf);

% ---- Listener: Rayleigh integral over a coarse patch grid.  Patches are grouped by
%      propagation delay, giving for each mode a short delay kernel W(:, k) [Pa per (m/s^2)].
npx = P.listener.patches(1); npy = P.listener.patches(2);
cx = (0.5:npx) * a / npx; cy = (0.5:npy) * b / npy;
[CX, CY] = meshgrid(cx, cy);
L = plate_listener(P.listener, CX, CY, (a / npx) * (b / npy), fs);
PhiL = sin(CX(:) * kx') .* sin(CY(:) * ky');              % npatch x nm
W = zeros(L.nbins, nm);
for j = 1:numel(CX)
    d = L.delay(j) + 1;
    W(d, :) = W(d, :) + L.weight(j) * PhiL(j, :);
end

% ---- Integrate the modes one at a time, accumulating the probe series and the radiated pressure
w_p = zeros(1, nt); wxx_p = zeros(1, nt); wyy_p = zeros(1, nt); wxy_p = zeros(1, nt);
p_audio = zeros(1, nt);
for k = 1:nm
    pk = (G(k) / M) * F;                                  % modal force per unit modal mass
    q = plate_sdof(omega(k), sigma(k), pk, dt);
    w_p = w_p + phi_p(k) * q;
    wxx_p = wxx_p + phixx_p(k) * q;
    wyy_p = wyy_p + phiyy_p(k) * q;
    wxy_p = wxy_p + phixy_p(k) * q;
    Qf(k, :) = q(frame_idx);
    % Radiated pressure.  The oscillator is linear and time-invariant, so the
    % delay-weighted response sum_d W(d,k) q''(t - d dt) equals the acceleration
    % response to the delay-weighted force sum_d W(d,k) pk(t - d dt).
    pk_delayed = conv(pk, W(:, k)');
    [~, ~, acc] = plate_sdof(omega(k), sigma(k), pk_delayed(1:nt), dt);
    p_audio = p_audio + acc;
end

% ---- Assemble outputs
Sp = plate_stress(P, wxx_p, wyy_p, wxy_p);
Wf = reshape(Phi * Qf, ngy, ngx, nf);
Sf = plate_stress(P, reshape(Phixx * Qf, ngy, ngx, nf), reshape(Phiyy * Qf, ngy, ngx, nf), reshape(Phixy * Qf, ngy, ngx, nf));
[fw, Aw] = plate_spectrum(t, w_p);
[~, Avm] = plate_spectrum(t, Sp.vm);

result.method = 'modal';
result.grid = struct('x', gx, 'y', gy);
result.t = t;
result.probe = struct('x', xp, 'y', yp, 'w', w_p, 'sxx', Sp.sxx, 'syy', Sp.syy, 'txy', Sp.txy, 'vm', Sp.vm, 'sp', Sp.sp);
result.frames = struct('t', t(frame_idx), 'w', Wf, 'vm', Sf.vm, 'sp', Sf.sp);
result.envelope = struct('vm', max(Sf.vm, [], 3), 'sp', max(abs(Sf.sp), [], 3));
result.spectrum = struct('f', fw, 'w', Aw, 'vm', Avm);
result.audio = struct('fs', fs, 'p', p_audio, 't0', L.t0);
result.modes = struct('m', m', 'n', n', 'f_hz', omega' / (2 * pi), 'sigma', sigma', 'G', G', 'truncated', truncated);
result.info = struct('nmodes', nm, 'fmax', P.modal.fmax, 'fs', fs, 'dt', dt, 'nsteps', nt, ...
    'ncells', ngx * ngy, 'stable', true, 'diverged_step', 0, 'wall_s', toc(t_start), 'notes', {notes});
end