1function result = plate_modal(params)
2% PLATE_MODAL Modal superposition for a simply supported rectangular plate.
3%
4% result = plate_modal(params) solves the damped Kirchhoff-Love equation
5%
6% w_tt + 2 sigma0 w_t - 2 sigma1 lap(w_t) + kappa^2 lap(lap(w)) = F(t) s(x,y) / mu
7%
8% on 0 < x < a, 0 < y < b with simply supported edges (w = lap(w) = 0). For
9% these edges the eigenfunctions are known exactly:
10%
11% phi_mn(x,y) = sin(m pi x / a) sin(n pi y / b)
12% omega_mn = kappa lambda_mn, lambda_mn = (m pi / a)^2 + (n pi / b)^2
13%
14% Expanding w = sum_mn q_mn(t) phi_mn(x,y) decouples the plate into
15% independent damped oscillators, one per mode,
16%
17% q_mn'' + 2 sigma_mn q_mn' + omega_mn^2 q_mn = G_mn F(t) / M_mn
18%
19% with modal mass M_mn = mu a b / 4, modal force factor G_mn = integral of
20% s(x,y) phi_mn over the footprint, and damping rate sigma_mn = sigma0 +
21% sigma1 lambda_mn. Each oscillator is integrated exactly for
22% piecewise-linear forcing (plate_sdof). Stresses follow from the analytic
23% second derivatives of the modes, and the radiated pressure from the modal
24% accelerations through the Rayleigh integral (plate_listener).
25%
26% The rectangle is the only shape with closed-form modes for these edges;
27% an ellipse is refused, and the finite-difference methods handle it.
28%
29% Approximations. The modal sum is truncated at params.modal.fmax and
30% params.modal.max_modes; the force is sampled at params.output.fs; the
31% footprint integral and the Rayleigh integral use quadrature grids. Stress
32% converges more slowly than displacement as modes are added, and the
33% damping model is a choice, not a measurement.
34%
35% See plate_result_schema for the fields of result.
37t_start = tic;
38P = plate_setup(params);
39if ~strcmp(P.bc, 'simply-supported')
40 error('plate_modal: closed-form modes exist only for simply supported edges (got "%s"); use the finite-difference method for other edges', P.bc);
41end
42if ~strcmp(P.plate.shape, 'rectangle')
43 error('plate_modal: closed-form modes exist only for the rectangle (got "%s"); use a finite-difference method for the ellipse', P.plate.shape);
44end
45a = P.plate.Lx;
46b = P.plate.Ly;
47fs = P.output.fs;
48dt = 1 / fs;
49nt = round(P.output.duration * fs) + 1; % samples spanning [0, duration]
50t = (0:nt - 1) * dt;
51notes = {};
53% ---- Mode selection: every mode below fmax, lowest first, at most max_modes
54omega_max = 2 * pi * P.modal.fmax;
55m_hi = max(1, ceil(a / pi * sqrt(omega_max / P.kappa)));
56n_hi = max(1, ceil(b / pi * sqrt(omega_max / P.kappa)));
57[MM, NN] = meshgrid(1:m_hi, 1:n_hi);
58lam_all = (MM * pi / a).^2 + (NN * pi / b).^2;
59keep = P.kappa * lam_all <= omega_max;
60m = MM(keep); n = NN(keep); lam = lam_all(keep);
61[lam, order] = sort(lam);
62m = m(order); n = n(order);
63truncated = numel(lam) > P.modal.max_modes;
64if truncated
65 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);
66 m = m(1:P.modal.max_modes); n = n(1:P.modal.max_modes); lam = lam(1:P.modal.max_modes);
67end
68m = m(:); n = n(:); lam = lam(:);
69nm = numel(lam);
70omega = P.kappa * lam;
71sigma = P.damping.sigma0 + P.damping.sigma1 * lam;
72over = sigma >= omega;
73if any(over)
74 notes{end + 1} = sprintf('%d modes would be overdamped by this damping model; their damping was capped just below critical', sum(over));
75 sigma(over) = 0.99 * omega(over);
76end
77M = P.mu * a * b / 4; % modal mass (same for every mode)
78kx = m * pi / a; % modal wavenumbers
79ky = n * pi / b;
81% ---- Applied force: total force F(t) and its projection G_mn onto each mode.
82% G_mn = integral s phi_mn dA over the footprint, by midpoint quadrature on a fine local grid.
83F = plate_force_time(P.excitation, t);
84exc = P.excitation;
85nq = 64;
86qx = linspace(max(0, exc.x - exc.radius), min(a, exc.x + exc.radius), nq);
87qy = linspace(max(0, exc.y - exc.radius), min(b, exc.y + exc.radius), nq);
88[QX, QY] = meshgrid(qx, qy);
89dAq = (qx(2) - qx(1)) * (qy(2) - qy(1));
90s = plate_force_footprint(exc, QX, QY, dAq);
91PhiQ = sin(QX(:) * kx') .* sin(QY(:) * ky'); % (nq*nq) x nm
92G = (s(:)' * PhiQ)' * dAq; % nm x 1
94% ---- Mode shapes and their second derivatives at the probe
95xp = P.probe.x; yp = P.probe.y;
96phi_p = sin(kx * xp) .* sin(ky * yp);
97phixx_p = -kx.^2 .* phi_p;
98phiyy_p = -ky.^2 .* phi_p;
99phixy_p = kx .* ky .* cos(kx * xp) .* cos(ky * yp);
101% ---- Display grid for animation frames
102ngx = P.output.grid(1); ngy = P.output.grid(2);
103gx = linspace(0, a, ngx); gy = linspace(0, b, ngy);
104[GX, GY] = meshgrid(gx, gy);
105Phi = sin(GX(:) * kx') .* sin(GY(:) * ky'); % (ngy*ngx) x nm
106Phixx = -Phi .* repmat(kx'.^2, ngx * ngy, 1);
107Phiyy = -Phi .* repmat(ky'.^2, ngx * ngy, 1);
108Phixy = (cos(GX(:) * kx') .* cos(GY(:) * ky')) .* repmat((kx .* ky)', ngx * ngy, 1);
109% Frames at the samples nearest to nframes equally spaced times in [0, duration]
110frame_idx = unique(min(nt, round(linspace(0, P.output.duration, P.output.nframes) * fs) + 1));
111nf = numel(frame_idx);
112Qf = zeros(nm, nf);
114% ---- Listener: Rayleigh integral over a coarse patch grid. Patches are grouped by
115% propagation delay, giving for each mode a short delay kernel W(:, k) [Pa per (m/s^2)].
116npx = P.listener.patches(1); npy = P.listener.patches(2);
117cx = (0.5:npx) * a / npx; cy = (0.5:npy) * b / npy;
118[CX, CY] = meshgrid(cx, cy);
119L = plate_listener(P.listener, CX, CY, (a / npx) * (b / npy), fs);
120PhiL = sin(CX(:) * kx') .* sin(CY(:) * ky'); % npatch x nm
121W = zeros(L.nbins, nm);
122for j = 1:numel(CX)
123 d = L.delay(j) + 1;
124 W(d, :) = W(d, :) + L.weight(j) * PhiL(j, :);
125end
127% ---- Integrate the modes one at a time, accumulating the probe series and the radiated pressure
128w_p = zeros(1, nt); wxx_p = zeros(1, nt); wyy_p = zeros(1, nt); wxy_p = zeros(1, nt);
129p_audio = zeros(1, nt);
130for k = 1:nm
131 pk = (G(k) / M) * F; % modal force per unit modal mass
132 q = plate_sdof(omega(k), sigma(k), pk, dt);
133 w_p = w_p + phi_p(k) * q;
134 wxx_p = wxx_p + phixx_p(k) * q;
135 wyy_p = wyy_p + phiyy_p(k) * q;
136 wxy_p = wxy_p + phixy_p(k) * q;
137 Qf(k, :) = q(frame_idx);
138 % Radiated pressure. The oscillator is linear and time-invariant, so the
139 % delay-weighted response sum_d W(d,k) q''(t - d dt) equals the acceleration
140 % response to the delay-weighted force sum_d W(d,k) pk(t - d dt).
141 pk_delayed = conv(pk, W(:, k)');
142 [~, ~, acc] = plate_sdof(omega(k), sigma(k), pk_delayed(1:nt), dt);
143 p_audio = p_audio + acc;
144end
146% ---- Assemble outputs
147Sp = plate_stress(P, wxx_p, wyy_p, wxy_p);
148Wf = reshape(Phi * Qf, ngy, ngx, nf);
149Sf = plate_stress(P, reshape(Phixx * Qf, ngy, ngx, nf), reshape(Phiyy * Qf, ngy, ngx, nf), reshape(Phixy * Qf, ngy, ngx, nf));
150[fw, Aw] = plate_spectrum(t, w_p);
151[~, Avm] = plate_spectrum(t, Sp.vm);
153result.method = 'modal';
154result.grid = struct('x', gx, 'y', gy);
155result.t = t;
156result.probe = struct('x', xp, 'y', yp, 'w', w_p, 'sxx', Sp.sxx, 'syy', Sp.syy, 'txy', Sp.txy, 'vm', Sp.vm, 'sp', Sp.sp);
157result.frames = struct('t', t(frame_idx), 'w', Wf, 'vm', Sf.vm, 'sp', Sf.sp);
158result.envelope = struct('vm', max(Sf.vm, [], 3), 'sp', max(abs(Sf.sp), [], 3));
159result.spectrum = struct('f', fw, 'w', Aw, 'vm', Avm);
160result.audio = struct('fs', fs, 'p', p_audio, 't0', L.t0);
161result.modes = struct('m', m', 'n', n', 'f_hz', omega' / (2 * pi), 'sigma', sigma', 'G', G', 'truncated', truncated);
162result.info = struct('nmodes', nm, 'fmax', P.modal.fmax, 'fs', fs, 'dt', dt, 'nsteps', nt, ...
163 'ncells', ngx * ngy, 'stable', true, 'diverged_step', 0, 'wall_s', toc(t_start), 'notes', {notes});
164end