function [p, radinfo] = bass_radiate(P, basis, tOn, str, stop, force, fsHi)
% BASS_RADIATE Pluck projection, closed-form synthesis, and free-field
% pressure at the listener.
%
% The pluck is an initial condition, not a force history. Releasing a
% string pulled with force F at fraction pluck_pos of the vibrating length
% from the bridge leaves the classic triangular shape, whose sine
% coefficients are b_k = 2 F Lv sin(k pi xp/Lv) / (T pi^2 k^2); projected
% through PM = (M Phi)' from bass_basis.m these become initial amplitudes
% of the coupled modes, with zero initial velocity. The finite release from
% the fingertip is a first-order rolloff 1/sqrt(1 + (om tau)^2) applied to
% the coupled amplitudes (finger_release_ms), which is what separates a
% fleshy jazz pizzicato from a hard snap. From there every coupled mode is a
% damped exponential; nothing is integrated. When the next note lands, the
% previous note is choked: its poles pick up an extra decay 1/tau_choke
% (choke_ms), standing in for the player stopping the string, and the last
% note rings through the full tail.
%
% Radiation, per coupled mode at its own centre frequency (the modes are
% narrow-band):
%   - the monopole volume Vm radiates from the two f-holes, half each, as
%     point monopoles with their own retardations;
%   - the dipole volume Dd radiates from the top-plate centre as a point
%     dipole along the plate normal (+z), with the near-field 1/r term kept;
%   - the direct term: above the retained body modes the bridge sees the
%     resistive residual admittance residual_y, so the bridge force (Fcoef
%     per unit modal amplitude) moves the plate broadband, and that motion
%     radiates as a small dipole. This is how the high partials reach the
%     listener at all, and it is the least principled term of the model;
%     direct_gain scales it and 0 turns it off.
% The common carrier exp(-i k rc) to the bridge is stripped, so relative
% phase between the f-holes, the plate, and the listener position is kept
% while the overall propagation delay is not double counted.
%
% Returns p (1 x T, pascals at fsHi) and radinfo diagnostics.

c = P.speed_of_sound;
rho0 = P.air_density;
K = length(tOn);
beat = 60 / max(20, P.tempo_bpm);
duration = (K - 1) * beat + P.tail_s;
T = floor(duration * fsHi);
dt = 1 / fsHi;
tauRel = P.finger_release_ms * 1e-3;
tauChoke = max(1e-3, P.choke_ms * 1e-3);

% ---- source and listener geometry (origin at the bridge) ----
lx = P.listener_x_m;
ly = P.listener_y_m;
lz = P.listener_z_m;
rc = sqrt(lx^2 + ly^2 + lz^2);
fhX = 0.095; fhY = -0.02; fhZ = 0.11;     % the two f-holes
dpX = 0.0; dpY = -0.10; dpZ = 0.11;       % top-plate dipole centre
r1 = max(0.05, sqrt((lx - fhX)^2 + (ly - fhY)^2 + (lz - fhZ)^2));
r2 = max(0.05, sqrt((lx + fhX)^2 + (ly - fhY)^2 + (lz - fhZ)^2));
rd = max(0.05, sqrt((lx - dpX)^2 + (ly - dpY)^2 + (lz - dpZ)^2));
cosg = (lz - dpZ) / rd;                    % dipole axis is the plate normal +z
listenerInside = rc < 0.35;
Sdir = 0.004 * P.rad_area_scale;

% ---- per-fingering radiation transfer, poles, and unit pluck amplitudes ----
nf = basis.nf;
ntot = basis.n;
G = zeros(ntot, 1) + 0i;
lam = zeros(ntot, 1) + 0i;
eta0u = zeros(ntot, 1);                    % per unit pluck force
for u = 1:nf
  n = basis.nModes(u);
  ns = basis.ns(u);
  Lv = basis.Lv(u);
  Tt = basis.T(u);
  mu = basis.mu(u);
  rows = basis.off(u) + (1:n);
  om = 2 * pi * basis.f(rows);
  zt = basis.zeta(rows);

  ph1 = exp(-1i * (om / c) * (r1 - rc));
  ph2 = exp(-1i * (om / c) * (r2 - rc));
  phd = exp(-1i * (om / c) * (rd - rc));
  near = 1i * (om / c) + 1 / rd;
  Gm = -(om.^2 * rho0 / (4 * pi)) .* (basis.Vm(rows) / 2) .* (ph1 / r1 + ph2 / r2);
  Gd = -(om.^2 * rho0 / (4 * pi)) .* basis.Dd(rows) .* near .* phd * (cosg / rd);
  Gr = (1i * om * rho0 / (4 * pi)) * (Sdir * P.residual_y * P.dipole_length_m) ...
    .* basis.Fcoef(rows) .* near .* phd * (cosg / rd);
  G(rows) = Gm + Gd + P.direct_gain * Gr;

  omd = om .* sqrt(max(1e-12, 1 - zt.^2));
  lam(rows) = -zt .* om + 1i * omd;

  % unit-force pluck projection through the coupled shapes
  xp = Lv * (1 - min(0.95, max(0.05, P.pluck_pos)));
  kv = (1:ns)';
  bk = (2 * Lv / (Tt * pi^2)) * sin(kv * pi * xp / Lv) ./ kv.^2;
  q0 = bk * sqrt(mu * Lv / 2);
  PM = reshape(basis.PMflat(basis.pmOff(u) + (1:n * ns)), n, ns);
  e0 = PM * q0;
  e0 = e0 ./ sqrt(1 + (om * tauRel).^2);   % finite finger release, first order
  eta0u(rows) = e0;
end

% ---- per-note synthesis: damped exponentials, choked at the next onset ----
ks = min(T, max(1, round(tOn * fsHi) + 1));
p = zeros(1, T);
CH = 32768;
chokeSpan = ceil(10 * tauChoke * fsHi);
peakEta = 0;
for i = 1:K
  u = basis.fIdx(i);
  n = basis.nModes(u);
  rows = basis.off(u) + (1:n);
  om = 2 * pi * basis.f(rows);
  zt = basis.zeta(rows);
  omd = om .* sqrt(max(1e-12, 1 - zt.^2));
  cAmp = (force(i) * eta0u(rows)) .* (1 - 1i * (zt .* om ./ omd));
  s = G(rows) .* cAmp;
  lamN = lam(rows);
  peakEta = max(peakEta, force(i) * max(abs(eta0u(rows))));

  if i < K
    segEnd = max(ks(i), ks(i + 1));
    stopAt = min(T + 1, segEnd + chokeSpan);
  else
    segEnd = T + 1;
    stopAt = T + 1;
  end

  pos = ks(i);
  while pos < segEnd
    len = min(CH, segEnd - pos);
    seg = real((s.') * exp(lamN * ((0:len - 1) * dt)));
    p(pos:pos + len - 1) = p(pos:pos + len - 1) + seg;
    s = s .* exp(lamN * (len * dt));
    pos = pos + len;
  end
  if pos < stopAt
    lamC = lamN - 1 / tauChoke;
    while pos < stopAt
      len = min(CH, stopAt - pos);
      seg = real((s.') * exp(lamC * ((0:len - 1) * dt)));
      p(pos:pos + len - 1) = p(pos:pos + len - 1) + seg;
      s = s .* exp(lamC * (len * dt));
      pos = pos + len;
    end
  end
end

% ---- diagnostics: the fundamental of the first note ----
u1 = basis.fIdx(1);
rows1 = basis.off(u1) + (1:basis.nModes(u1));
a1 = abs(eta0u(rows1));
best = 1;
bv = 0;
for m = 1:length(rows1)
  if a1(m) > bv
    bv = a1(m);
    best = m;
  end
end
f1 = basis.f(rows1(best));
z1 = basis.zeta(rows1(best));
radinfo.fundHz = f1;
radinfo.fundT60 = log(1000) / max(1e-9, z1 * 2 * pi * f1);
radinfo.listenerDistM = rc;
radinfo.listenerInside = listenerInside;
radinfo.gainAbs = abs(G);
radinfo.peakEta = peakEta;
end