function basis = maraca_basis(P)
% MARACA_BASIS Mode set of a complete thin spherical shell.
%
% Membrane theory of the full sphere gives two branches per degree l from
%   Omega^4 - (lam + 1 + 3 nu) Omega^2 + (lam - 2)(1 - nu^2) = 0,
% lam = l(l+1), Omega = omega R / c_p, c_p = sqrt(E / (rho (1 - nu^2))).
% Checks: l = 0 upper root is the breathing mode Omega^2 = 2(1+nu); the l = 1
% lower root is exactly 0 (rigid translation, held by the hand and neglected).
% The lower branch flattens toward Omega^2 = 1 - nu^2 at high l, so a bending
% term beta^2 (lam - 2)^2, beta^2 = h^2/(12 R^2), is added to it; this is the
% flat-plate asymptote and is what spreads the high-l modes upward in
% frequency. The membrane eigenvector is kept for the bending-corrected modes,
% an approximation that is good while the bending term is a correction.
%
% Eigenvector ratio (tangential potential V over radial W, W = 1):
%   V = (2(1+nu) - Omega^2) / (lam (1+nu)),
% giving the modal mass factor mu = 1 + lam V^2 and modal mass rho h R^2 mu.
%
% The closed cavity enters as an air stiffness on the l = 0 (breathing)
% components only: d(omega^2) = 3 rho0 c^2 / (rho h R), scaled by
% cavity_coupling. Modes with l >= 1 displace no net volume, and the interior
% inertial loading they see is of order rho0 R / ((l+1) rho h), well under a
% percent here, so it is neglected.
%
% Damping: zeta = 0.5 (loss_factor + loss_slope omega) plus radiation loss
% rho0 c sigma_l / (2 rho h omega mu), sigma_l = Re[-i h_l(ka)/h_l'(ka)].

R = P.shell_radius_m;
h = P.shell_thickness_m;
nu = P.poisson_ratio;
rho = P.density;
cp = sqrt(P.youngs_modulus / (rho * (1 - nu^2)));
beta2 = h^2 / (12 * R^2);
c = P.speed_of_sound;
rho0 = P.air_density;

cavShift2 = P.cavity_coupling * 3 * rho0 * c^2 / (rho * h * R);   % rad^2/s^2

maxN = 2 * (P.max_mode_l + 1);
ll = zeros(maxN, 1);
ff = zeros(maxN, 1);
mus = zeros(maxN, 1);
n = 0;
breathingDryHz = 0;
breathingHz = 0;
for l = 0:P.max_mode_l
  lam = l * (l + 1);
  b = lam + 1 + 3 * nu;
  det2 = (lam - 2) * (1 - nu^2);
  disc = sqrt(max(0, b^2 - 4 * det2));
  om2lo = (b - disc) / 2;
  om2hi = (b + disc) / 2;
  if l == 0
    cand = om2hi;                       % breathing; lower root is negative
  elseif l == 1
    cand = om2hi;                       % lower root 0 = rigid translation
  else
    cand = [om2lo + beta2 * (lam - 2)^2, om2hi];
  end
  for ci = 1:length(cand)
    Om2 = cand(ci);
    om = sqrt(Om2) * cp / R * P.stiffness_trim;
    if l == 0
      breathingDryHz = om / (2 * pi);
      om = sqrt(om^2 + cavShift2);      % closed-cavity air stiffness
      breathingHz = om / (2 * pi);
    end
    fm = om / (2 * pi);
    if fm > P.max_mode_hz
      continue;
    end
    if lam > 0
      V = (2 * (1 + nu) - Om2) / (lam * (1 + nu));
    else
      V = 0;
    end
    n = n + 1;
    ll(n) = l;
    ff(n) = fm;
    mus(n) = 1 + lam * V^2;
  end
end
if n == 0
  error('No shell modes below max_mode_hz -- check geometry/material.');
end
[ff, order] = sort(ff(1:n, 1));
ff = reshape(ff, n, 1);
lls = reshape(ll(order), n, 1);
muss = reshape(mus(order), n, 1);

omega = 2 * pi * ff;
sigma = zeros(n, 1);
for m = 1:n
  ka = omega(m) / c * R;
  [hh, hhp] = maraca_hankel(lls(m), ka);
  sigma(m) = real(-1i * hh(lls(m) + 1) / hhp(lls(m) + 1));
end
zeta = 0.5 * (P.loss_factor + P.loss_slope * omega) ...
  + P.radiation_loss * rho0 * c * sigma ./ (2 * rho * h * omega .* muss);

basis.n = n;
basis.l = lls;
basis.f = ff;
basis.omega = omega;
basis.zeta = zeta;
basis.mu = muss;
basis.massKg = rho * h * R^2 * muss;
basis.sigma = sigma;
basis.breathingDryHz = breathingDryHz;
basis.breathingHz = breathingHz;
end