function basis = cymbal_basis(P)
% CYMBAL_BASIS Mode set of the cymbal: one shell eigensolve per nodal-diameter
% order, plus damping and the radiation data the later stages need.
%
% Because the plate is a shell of revolution, the modes separate exactly in
% the azimuth (see cymbal_shell.m), so the two-dimensional eigenproblem breaks
% into n_diameters + 1 one-dimensional ones. Each is reduced to symmetric
% standard form through the Cholesky factor of the consistent mass matrix and
% diagonalized by svd rather than eig: the reduced matrix is symmetric
% positive semidefinite, so its singular vectors are its eigenvectors, and svd
% keeps them orthonormal for the closely spaced pairs this geometry produces.
% The eigenvectors come back mass-normalized, phi' M phi = I. Rigid-body modes
% (two at n = 0, two at n = 1, none above) sit at zero frequency and are
% dropped with everything below min_mode_hz; modes above max_mode_hz are
% truncated, and the broadband part they carry is reinstated by the direct
% term in cymbal_radiate.m.
%
% Only the cos(n theta) member of each degenerate pair is retained. That is
% exact for the model as posed: every contact sits on the meridian theta = 0,
% where the sin member has a node.
%
% HAMMERING. A perfectly axisymmetric plate is a bad cymbal, and the model has
% to say so. Its modes are exact cos(n theta) shapes, and a mode with n nodal
% diameters radiates only into the cone where its azimuthal trace speed is
% supersonic, n < k a sin(theta). Below a few kilohertz that leaves only the
% handful of n = 0, 1, 2 modes audible at any listening position, so an ideal
% cymbal sounds like a small bell. Real cymbals are hammered and lathed, and
% that is not decoration: the irregularity it leaves destroys the exact
% separation, and every mode of the real plate is a hybrid of the ideal ones
% that lie near it in frequency. A hybrid carries low-order content, so it
% radiates.
%
% The hybridization is first-order degenerate perturbation theory, applied as
% a seeded orthogonal mixing of modes lying within hammer_bw of each other in
% frequency (a sequence of Givens rotations on neighbouring pairs, which is
% orthogonal by construction and therefore conserves energy exactly), plus a
% seeded relative detune of hammer_detune, which is what makes near neighbours
% beat. The mixing matrix is stored as basis.mix and applied downstream
% wherever the shapes enter linearly: to the stroke's modal projection in
% cymbal_cascade.m and to the radiation transfer in cymbal_radiate.m. Because
% the window is relative, the sparse low modes stay pure and only the dense
% region above a kilohertz hybridizes, which is where a real cymbal's wash
% lives. hammer_mix = 0 recovers the ideal axisymmetric plate.
%
% Damping has three parts.
%   internal   damp_rate / om + 0.5 (loss_factor + loss_slope om). The first
%              term is a frequency-independent decay RATE, which is what a
%              resistive path out of the plate gives and what keeps the low
%              modes of a cymbal ringing for tens of seconds; the other two are
%              the usual constant and frequency-proportional loss factors, and
%              they are what makes the shimmer die well before the hum does;
%   mount      a point dashpot of mount_loss N s/m on the axial motion at the
%              inner edge, standing in for the felt washer and the stand. With
%              mass-normalized shapes its modal damping is c phi_axial^2 / 2 om,
%              which falls off steeply with n because the shapes of high
%              nodal-diameter modes vanish at the hole. This is what kills the
%              low n = 0 gong modes a mounted cymbal does not sustain;
%   radiation  computed from the actual radiated power, below.
%
% Radiated power is evaluated in the far field, where it is defined, rather
% than by summing element powers as the sibling triangle does (an incoherent
% element sum is mesh dependent for a continuous surface: halving the element
% size halves the estimate). The plate is treated as a sheet of normal dipoles
% of volume dS * ldip, and for a cos(n theta) source on a surface of revolution
% the azimuthal integral of the far-field kernel is a Bessel function, so the
% far field is exactly
%   p(th_o, ph_o) = -(om^2 rho0 / (4 pi c R)) e^{-ikR} A_n(th_o) cos(n ph_o)
%   A_n = ldip sum_e r_e L_e W_e e^{i k z_e cos th_o}
%         [ 2 pi i^n J_n(al) r'_e cos th_o
%           - pi (i^{n+1} J_{n+1}(al) + i^{n-1} J_{n-1}(al)) z'_e sin th_o ]
%   al  = k r_e sin th_o
% and the power integral over the sphere collapses to one quadrature in
% cos th_o. The cos(n ph_o) factor is the model's multipole structure made
% explicit: a mode with n nodal diameters is a 2^n-pole, so high-n modes
% radiate very little at low frequency, and none of them radiate on the axis.
%
% ldip is the one phenomenological length in the radiation model: the distance
% over which the air short-circuits from one face of the plate to the other.
% It is taken as half the shorter of the flexural and acoustic wavelengths,
% capped at the plate radius, and scaled by dipole_scale. Below the plate's
% coincidence frequency (about 18 kHz for 1 mm bronze) it is the flexural
% wavelength that sets it, which is why a cymbal is such a poor radiator in
% the low register and gets steadily better as the cascade fills the top.

G = cymbal_profile(P);
E = P.youngs_modulus;
nu = P.poisson_ratio;
rho = P.density;
c = P.speed_of_sound;
rho0 = P.air_density;
nn = G.nn;
nel = nn - 1;
Lref = G.S / nel;
nmax = max(0, round(P.n_diameters));

% ---- radiation sample points: two per element, Hermite-evaluated ----
nsub = 2;
nsurf = nel * nsub;
surfS = zeros(nsurf, 1);
surfLen = zeros(nsurf, 1);
surfEl = zeros(nsurf, 1);
surfXi = zeros(nsurf, 1);
for e = 1:nel
  L = G.snode(e + 1) - G.snode(e);
  for j = 1:nsub
    idx = (e - 1) * nsub + j;
    xi = (j - 0.5) / nsub;
    surfS(idx) = G.snode(e) + xi * L;
    surfLen(idx) = L / nsub;
    surfEl(idx) = e;
    surfXi(idx) = xi;
  end
end
surfR = reshape(interp1(G.s, G.r, surfS), nsurf, 1);
surfZ = reshape(interp1(G.s, G.z, surfS), nsurf, 1);
surfDr = reshape(interp1(G.s, G.drds, surfS), nsurf, 1);
surfDz = reshape(interp1(G.s, G.dzds, surfS), nsurf, 1);

% ---- eigensolve per circumferential order ----
fs = [];
nd = [];
NUall = [];
NWall = [];
NTall = [];
SWall = [];
omLo = 2 * pi * P.min_mode_hz;
omHi = 2 * pi * P.max_mode_hz;
for n = 0:nmax
  [K, M] = cymbal_shell(G, n, E, nu, rho, Lref);
  Rc = chol(M);
  Am = (Rc') \ (K / Rc);
  Am = (Am + Am') / 2;
  [Uv, Sv, ~] = svd(Am);
  om2 = flipud(diag(Sv));
  Vv = fliplr(Uv);
  om = sqrt(max(0, om2)) * P.stiffness_trim;
  keep = find(om > omLo & om < omHi);
  if length(keep) == 0
    continue;
  end
  Phi = Rc \ Vv(:, keep);
  nk = length(keep);
  % svd fixes each singular vector only up to sign, and different linear
  % algebra libraries make different choices. On its own that would not
  % matter, since flipping one mode flips both its excitation and its
  % radiation and the two cancel; but the hammering below forms SUMS of
  % modes, where they do not, and the render would then depend on which
  % machine it was computed on. Canonicalize instead: the largest component
  % of every shape is positive, with ties broken by index so that the rule
  % itself is deterministic.
  for k = 1:nk
    av = abs(Phi(:, k));
    hit = find(av >= max(av) * (1 - 1e-9));
    if Phi(hit(1), k) < 0
      Phi(:, k) = -Phi(:, k);
    end
  end
  Uz = Phi(1:4:4 * nn, :);              % meridional, nn x nk
  Wz = Phi(3:4:4 * nn, :);              % normal
  Tz = Phi(4:4:4 * nn, :) * Lref;       % dW/ds, unscaled back to physical
  % normal displacement at the radiation sample points (Hermite)
  SW = zeros(nsurf, nk);
  for idx = 1:nsurf
    e = surfEl(idx);
    xi = surfXi(idx);
    L = G.snode(e + 1) - G.snode(e);
    h1 = 1 - 3 * xi^2 + 2 * xi^3;
    h2 = L * (xi - 2 * xi^2 + xi^3);
    h3 = 3 * xi^2 - 2 * xi^3;
    h4 = L * (-xi^2 + xi^3);
    SW(idx, :) = h1 * Wz(e, :) + h2 * Tz(e, :) + h3 * Wz(e + 1, :) + h4 * Tz(e + 1, :);
  end
  fs = [fs; om(keep) / (2 * pi)];
  nd = [nd; n * ones(nk, 1)];
  NUall = [NUall, Uz];
  NWall = [NWall, Wz];
  NTall = [NTall, Tz];
  SWall = [SWall, SW];
end
N = length(fs);
if N == 0
  error('No cymbal modes between min_mode_hz and max_mode_hz -- check geometry/material.');
end
[fs, order] = sort(fs);
fs = reshape(fs, N, 1);
nd = reshape(nd(order), N, 1);
NUall = NUall(:, order);
NWall = NWall(:, order);
NTall = NTall(:, order);
SWall = SWall(:, order);

% ---- hammering: detune, then hybridize near-degenerate modes ----
hstate = mod(floor(abs(P.hammer_seed)), 4294967296);
[hu, hstate] = cymbal_rng(hstate, 2 * N);
hz = sqrt(-2 * log(hu(1:N))) .* cos(2 * pi * hu(N + 1:2 * N));
fs = fs .* (1 + P.hammer_detune * hz);
[fs, ord2] = sort(fs);
fs = reshape(fs, N, 1);
nd = reshape(nd(ord2), N, 1);
NUall = NUall(:, ord2);
NWall = NWall(:, ord2);
NTall = NTall(:, ord2);
SWall = SWall(:, ord2);

HPASS = 4;
mix = eye(N);
if P.hammer_mix > 0 && N > 1
  [hu2, hstate] = cymbal_rng(hstate, HPASS * (N - 1));
  q = 0;
  for pass = 1:HPASS
    for i = 1:N - 1
      q = q + 1;
      if (fs(i + 1) - fs(i)) / fs(i) > P.hammer_bw
        continue;
      end
      th = P.hammer_mix * (hu2(q) - 0.5) * pi;
      ca = cos(th);
      sa = sin(th);
      ri = mix(i, :);
      rj = mix(i + 1, :);
      mix(i, :) = ca * ri - sa * rj;
      mix(i + 1, :) = sa * ri + ca * rj;
    end
  end
end
mix2 = mix.^2;

omega = 2 * pi * fs;

% ---- the acoustic short-circuit length, per mode ----
tMean = 2 * pi * trapz(G.s, G.r .* G.thk) / G.area;
cL = sqrt(E / (rho * (1 - nu^2)));
kap = tMean / sqrt(12);
kFlex = sqrt(omega / (kap * cL));
kAc = omega / c;
ldip = P.dipole_scale * min(P.radius_m, pi ./ max(kFlex, kAc));

% ---- radiated power, and the radiation damping that follows from it ----
nMu = 24;
mu = reshape(((1:nMu) - 0.5) / nMu * 2 - 1, 1, nMu);   % cos(theta_o) midpoints
dmu = 2 / nMu;
sth = sqrt(max(0, 1 - mu.^2));
zetaRad = zeros(N, 1);
sigmaEff = zeros(N, 1);
for m = 1:N
  om = omega(m);
  k = om / c;
  nn_d = nd(m);
  alpha = (k * surfR) * sth;                       % nsurf x nMu
  Jn = besselj(nn_d, alpha);
  Jp1 = besselj(nn_d + 1, alpha);
  if nn_d == 0
    Jm1 = -Jp1;
  else
    Jm1 = besselj(nn_d - 1, alpha);
  end
  % i^(n+1) J_{n+1} + i^(n-1) J_{n-1} = i^n * i * (J_{n+1} - J_{n-1})
  ipow = 1i^nn_d;
  Tk = ipow * (2 * pi * Jn .* (surfDr * mu) ...
               - 1i * pi * (Jp1 - Jm1) .* (surfDz * sth));
  carrier = exp(1i * (k * surfZ) * mu);
  wgt = ldip(m) * surfR .* surfLen .* SWall(:, m);
  An = sum((wgt * ones(1, nMu)) .* Tk .* carrier, 1);   % 1 x nMu
  Kn = sum(abs(An).^2) * dmu;
  if nn_d == 0
    piN = 2 * pi;
  else
    piN = pi;
  end
  zetaRad(m) = om^3 * rho0 * piN * Kn / (32 * pi^2 * c^3);
  % radiation efficiency, for the diagnostics only
  vv = sum(surfR .* surfLen .* SWall(:, m).^2) * piN;
  sigmaEff(m) = 2 * zetaRad(m) * om / max(1e-30, rho0 * c * vv);
end

% ---- mount: a point dashpot on the axial motion at the inner edge ----
axHole = NUall(1, :).' * G.dzdsnode(1) + NWall(1, :).' * G.drdsnode(1);
zetaMount = P.mount_loss * axHole.^2 ./ (2 * omega);

% The shape-dependent damping terms belong to the ideal modes, so they carry
% over to the hybridized ones as the diagonal of the transformed damping
% matrix, which to first order is the mixing weights squared.
zetaRad = mix2 * zetaRad;
zetaMount = mix2 * zetaMount;
sigmaEff = mix2 * sigmaEff;

zeta = P.damp_rate ./ omega ...
     + 0.5 * (P.loss_factor + P.loss_slope * omega) ...
     + P.radiation_loss * zetaRad + zetaMount;

basis.n = N;
basis.f = fs;
basis.omega = omega;
basis.zeta = zeta;
basis.zetaRad = zetaRad;
basis.zetaMount = zetaMount;
basis.sigma = sigmaEff;
basis.ndia = nd;
basis.ldip = ldip;
basis.mix = mix;
basis.NU = NUall;
basis.NW = NWall;
basis.NT = NTall;
basis.rnode = reshape(G.rnode, nn, 1);
basis.znode = reshape(G.znode, nn, 1);
basis.snode = reshape(G.snode, nn, 1);
basis.drdsnode = reshape(G.drdsnode, nn, 1);
basis.dzdsnode = reshape(G.dzdsnode, nn, 1);
basis.thknode = reshape(G.thknode, nn, 1);
basis.csnode = reshape(interp1(G.s, G.cs, G.snode), nn, 1);
basis.ctnode = reshape(interp1(G.s, G.ct, G.snode), nn, 1);
basis.surfR = surfR;
basis.surfZ = surfZ;
basis.surfDr = surfDr;
basis.surfDz = surfDz;
basis.surfLen = surfLen;
basis.SW = SWall;
basis.S = G.S;
basis.area = G.area;
basis.mass = G.mass;
basis.thkMean = tMean;
end