1function basis = maraca_basis(P)
2% MARACA_BASIS Mode set of a complete thin spherical shell.
3%
4% Membrane theory of the full sphere gives two branches per degree l from
5% Omega^4 - (lam + 1 + 3 nu) Omega^2 + (lam - 2)(1 - nu^2) = 0,
6% lam = l(l+1), Omega = omega R / c_p, c_p = sqrt(E / (rho (1 - nu^2))).
7% Checks: l = 0 upper root is the breathing mode Omega^2 = 2(1+nu); the l = 1
8% lower root is exactly 0 (rigid translation, held by the hand and neglected).
9% The lower branch flattens toward Omega^2 = 1 - nu^2 at high l, so a bending
10% term beta^2 (lam - 2)^2, beta^2 = h^2/(12 R^2), is added to it; this is the
11% flat-plate asymptote and is what spreads the high-l modes upward in
12% frequency. The membrane eigenvector is kept for the bending-corrected modes,
13% an approximation that is good while the bending term is a correction.
14%
15% Eigenvector ratio (tangential potential V over radial W, W = 1):
16% V = (2(1+nu) - Omega^2) / (lam (1+nu)),
17% giving the modal mass factor mu = 1 + lam V^2 and modal mass rho h R^2 mu.
18%
19% The closed cavity enters as an air stiffness on the l = 0 (breathing)
20% components only: d(omega^2) = 3 rho0 c^2 / (rho h R), scaled by
21% cavity_coupling. Modes with l >= 1 displace no net volume, and the interior
22% inertial loading they see is of order rho0 R / ((l+1) rho h), well under a
23% percent here, so it is neglected.
24%
25% Damping: zeta = 0.5 (loss_factor + loss_slope omega) plus radiation loss
26% rho0 c sigma_l / (2 rho h omega mu), sigma_l = Re[-i h_l(ka)/h_l'(ka)].
28R = P.shell_radius_m;
29h = P.shell_thickness_m;
30nu = P.poisson_ratio;
31rho = P.density;
32cp = sqrt(P.youngs_modulus / (rho * (1 - nu^2)));
33beta2 = h^2 / (12 * R^2);
34c = P.speed_of_sound;
35rho0 = P.air_density;
37cavShift2 = P.cavity_coupling * 3 * rho0 * c^2 / (rho * h * R); % rad^2/s^2
39maxN = 2 * (P.max_mode_l + 1);
40ll = zeros(maxN, 1);
41ff = zeros(maxN, 1);
42mus = zeros(maxN, 1);
43n = 0;
44breathingDryHz = 0;
45breathingHz = 0;
46for l = 0:P.max_mode_l
47 lam = l * (l + 1);
48 b = lam + 1 + 3 * nu;
49 det2 = (lam - 2) * (1 - nu^2);
50 disc = sqrt(max(0, b^2 - 4 * det2));
51 om2lo = (b - disc) / 2;
52 om2hi = (b + disc) / 2;
53 if l == 0
54 cand = om2hi; % breathing; lower root is negative
55 elseif l == 1
56 cand = om2hi; % lower root 0 = rigid translation
57 else
58 cand = [om2lo + beta2 * (lam - 2)^2, om2hi];
59 end
60 for ci = 1:length(cand)
61 Om2 = cand(ci);
62 om = sqrt(Om2) * cp / R * P.stiffness_trim;
63 if l == 0
64 breathingDryHz = om / (2 * pi);
65 om = sqrt(om^2 + cavShift2); % closed-cavity air stiffness
66 breathingHz = om / (2 * pi);
67 end
68 fm = om / (2 * pi);
69 if fm > P.max_mode_hz
70 continue;
71 end
72 if lam > 0
73 V = (2 * (1 + nu) - Om2) / (lam * (1 + nu));
74 else
75 V = 0;
76 end
77 n = n + 1;
78 ll(n) = l;
79 ff(n) = fm;
80 mus(n) = 1 + lam * V^2;
81 end
82end
83if n == 0
84 error('No shell modes below max_mode_hz -- check geometry/material.');
85end
86[ff, order] = sort(ff(1:n, 1));
87ff = reshape(ff, n, 1);
88lls = reshape(ll(order), n, 1);
89muss = reshape(mus(order), n, 1);
91omega = 2 * pi * ff;
92sigma = zeros(n, 1);
93for m = 1:n
94 ka = omega(m) / c * R;
95 [hh, hhp] = maraca_hankel(lls(m), ka);
96 sigma(m) = real(-1i * hh(lls(m) + 1) / hhp(lls(m) + 1));
97end
98zeta = 0.5 * (P.loss_factor + P.loss_slope * omega) ...
99 + P.radiation_loss * rho0 * c * sigma ./ (2 * rho * h * omega .* muss);
101basis.n = n;
102basis.l = lls;
103basis.f = ff;
104basis.omega = omega;
105basis.zeta = zeta;
106basis.mu = muss;
107basis.massKg = rho * h * R^2 * muss;
108basis.sigma = sigma;
109basis.breathingDryHz = breathingDryHz;
110basis.breathingHz = breathingHz;
111end