function basis = bass_basis(P, str, stop)
% BASS_BASIS Coupled string-body mode set, one eigensolve per fingering.
%
% For each fingering of the gesture (a string and a stop, deduplicated so a
% walking line that revisits a note reuses its solve), the vibrating string
% and the modal body of bass_body.m are assembled into ONE symmetric
% generalized eigenproblem and diagonalized together, so the coupling
% through the bridge is two-way: the string drives the body, and the moving
% bridge detunes and damps the string in return. This is where wolf-note
% physics lives; a note whose fundamental lands on a strong body resonance
% comes out as a split, fast-beating pair of coupled modes rather than one.
%
% The displacement ansatz is
%   u(x) = sum_k q_k phi_k(x) + y_b x/Lv,     y_b = sum_j a_j xi_j,
% with phi_k the mass-normalized pinned-pinned modes of the stiff string
% (om_k^2 = (k pi/Lv)^2 T/mu + (k pi/Lv)^4 EI/mu; the bending stiffness EI
% comes from the inharmonicity parameter, EI = B T L0^2/pi^2, and is what
% stretches the high partials sharp) and x/Lv the static shape the string
% takes when the bridge end moves. With this split the elastic coupling
% integrals vanish (the ramp has no curvature and phi_k vanishes at both
% ends), so all the coupling sits in the mass matrix, plus the ramp's own
% inertia and the tension load T/Lv on the body block; bass_assemble.m
% builds the matrices. Both are symmetric and M is positive definite, so
% the problem reduces through the Cholesky factor of M and is diagonalized
% by svd, which keeps near-degenerate coupled pairs orthonormal where a
% nonsymmetric eigensolver does not. The shapes come back mass-normalized.
%
% Damping is assigned in the uncoupled coordinates and projected,
% zeta_m = (Phi' C0 Phi)_mm / (2 om_m): string internal loss (loss_factor +
% loss_slope om), the finger termination on stopped notes (finger_loss), the
% resistive residual bridge admittance above the retained body modes
% (residual_y, high-passed so it does not double-count the modal region),
% and the body Qs. Radiation damping from each coupled mode's own monopole
% and dipole strengths is added on top, scaled by radiation_loss.
%
% Stored per coupled mode: frequency, damping, monopole volume Vm and dipole
% volume Dd per unit modal coordinate (what bass_radiate.m radiates), the
% bridge force coefficient Fcoef (what drives the direct term), and the
% projection matrix PM = (M Phi)' restricted to the string block, which
% turns a plucked initial shape into coupled initial amplitudes. Fields are
% concatenated across fingerings with offset indexing (numbl has no struct
% arrays).

rho0 = P.air_density;
cair = P.speed_of_sound;

% ---- deduplicate fingerings ----
K = length(str);
fkey = str * 100 + stop;
uKey = zeros(1, K);
fIdx = zeros(1, K);
nf = 0;
for i = 1:K
  hit = 0;
  for u = 1:nf
    if uKey(u) == fkey(i)
      hit = u;
    end
  end
  if hit == 0
    nf = nf + 1;
    uKey(nf) = fkey(i);
    hit = nf;
  end
  fIdx(i) = hit;
end

% ---- per-fingering eigensolve, concatenated storage ----
f = [];
zeta = [];
Vm = [];
Dd = [];
Fcoef = [];
ybAmp = [];
bodyFrac = [];
PMflat = [];
off = zeros(nf, 1);
nModes = zeros(nf, 1);
pmOff = zeros(nf, 1);
nsArr = zeros(nf, 1);
LvArr = zeros(nf, 1);
muArr = zeros(nf, 1);
TArr = zeros(nf, 1);
f0Arr = zeros(nf, 1);
strArr = zeros(nf, 1);
stopArr = zeros(nf, 1);
body = bass_body(P);

for u = 1:nf
  s = floor(uKey(u) / 100);
  st = uKey(u) - 100 * s;
  [M, Kmat, omS, zS, meta] = bass_assemble(P, s, st);
  ns = meta.ns;
  nb = meta.nb;
  n = meta.n;
  Lv = meta.Lv;
  mu = meta.mu;
  T = meta.T;
  kv = (1:ns)';

  Rc = chol(M);
  Am = (Rc') \ (Kmat / Rc);
  Am = (Am + Am') / 2;
  [Uv, Sv, ~] = svd(Am);
  omc2 = flipud(diag(Sv));            % ascending
  Vv = fliplr(Uv);
  omc = sqrt(max(1e-6, omc2));
  Phi = Rc \ Vv;                      % mass-normalized coupled shapes, n x n

  % projected damping plus per-mode radiation damping
  cvec = 2 * [zS .* omS; body.zeta .* body.omega];
  zc = ((cvec') * (Phi.^2))' ./ (2 * omc);
  PhiB = Phi(ns + 1:n, :);
  bodyDisp = (body.a * ones(1, n)) .* PhiB;          % nb x n physical
  vm = ((body.S .* body.mono)' * bodyDisp)';         % n x 1 monopole volume
  dd = P.dipole_length_m * ((body.S .* (1 - body.mono))' * bodyDisp)';
  zRad = P.radiation_loss * (rho0 * omc.^4 .* vm.^2 / (8 * pi * cair) ...
    + rho0 * omc.^6 .* dd.^2 / (24 * pi * cair^3)) ./ omc.^3;
  zc = zc + zRad;

  yb = ((body.a') * PhiB)';                          % n x 1 bridge displacement
  phiSlope = sqrt(2 / (mu * Lv)) * (kv * pi / Lv) .* ((-1).^kv);
  fc = T * (((phiSlope') * Phi(1:ns, :))' + yb / Lv);
  bfrac = ((ones(1, nb) * (PhiB.^2)))';

  PM = (M * Phi);
  PM = (PM(1:ns, :))';                               % n x ns projection

  off(u) = length(f);
  nModes(u) = n;
  pmOff(u) = length(PMflat);
  nsArr(u) = ns;
  LvArr(u) = Lv;
  muArr(u) = mu;
  TArr(u) = T;
  f0Arr(u) = meta.f0;
  strArr(u) = s;
  stopArr(u) = st;

  f = [f; omc / (2 * pi)];
  zeta = [zeta; zc];
  Vm = [Vm; vm];
  Dd = [Dd; dd];
  Fcoef = [Fcoef; fc];
  ybAmp = [ybAmp; yb];
  bodyFrac = [bodyFrac; bfrac];
  PMflat = [PMflat; reshape(PM, n * ns, 1)];
end

basis.nf = nf;
basis.fIdx = fIdx;
basis.off = off;
basis.nModes = nModes;
basis.pmOff = pmOff;
basis.ns = nsArr;
basis.Lv = LvArr;
basis.mu = muArr;
basis.T = TArr;
basis.f0 = f0Arr;
basis.str = strArr;
basis.stop = stopArr;
basis.f = f;
basis.zeta = zeta;
basis.Vm = Vm;
basis.Dd = Dd;
basis.Fcoef = Fcoef;
basis.ybAmp = ybAmp;
basis.bodyFrac = bodyFrac;
basis.PMflat = PMflat;
basis.n = length(f);
basis.bodyF = body.f;
basis.bodyZeta = body.zeta;
basis.bodyA = body.a;
basis.bodyS = body.S;
basis.bodyMono = body.mono;
end