function [p, radinfo] = cymbal_radiate(P, basis, env, tframe, DV, kf, tstr, taustr, fstr, sstr, fsHi)
% CYMBAL_RADIATE Free-field pressure at the listener from the modal envelopes.
%
% Two contributions:
%
% 1. The mode bank. Each retained mode is a normal-dipole sheet on the plate:
%    a surface element of area dS carries a dipole of volume dS * ldip along
%    the local normal (ldip from cymbal_basis.m), and the element dipoles are
%    summed coherently with their relative retardations at the mode's centre
%    frequency,
%      G = sum_e (i om rho0 / 4 pi) V_e (u_e . rhat_e) (i k + 1/r_e)
%            exp(-i k (r_e - rcen)) / r_e,
%    over n_azimuth azimuths and two sample points per meridian element. The
%    common carrier exp(-i k rcen) is stripped, so relative phase across the
%    plate is kept while the overall propagation delay is not double counted.
%    Unlike the far-field form used for the radiated power, this keeps the
%    1/r^2 near-field term and the true distances, which matters: a listener
%    half a metre from a 40 cm plate is nowhere near the far field. Only the
%    normal component of the surface motion radiates; the in-plane component
%    of a thin plate displaces almost no air.
%
%    Because the azimuthal shape is exactly cos(n theta), the sum over
%    azimuths carries the multipole cancellation of the ideal mode: on the
%    axis of symmetry every mode with n >= 1 cancels exactly, and off axis a
%    mode radiates only into the cone where its azimuthal trace is supersonic,
%    n < k a sin(theta). That cutoff is why an unhammered plate is heard as a
%    handful of low-order partials; the hybridization applied below is what
%    gives the rest of the modes a way out. Moving the listener re-runs only
%    this stage.
%
%    The synthesis itself is an oscillator bank rather than the closed-form
%    damped exponentials of the sibling repositories: the cascade makes every
%    amplitude time-varying, so each mode is a fixed-frequency oscillator
%    whose envelope is the cascade solution, interpolated linearly across each
%    control frame. Within a frame the two matrix products below are exact,
%    and the per-frame phase advance is a unit-modulus multiply, so nothing
%    drifts.
%
%    Amplitude comes from the cascade; PHASE comes from a linear impulse
%    response carried alongside it. Each oscillator keeps a unit phasor and the
%    magnitude the mode would have had with no cascade at all; a contact adds
%    its signed modal velocity jump to that complex state and the phasor is
%    renormalized. A mode that was quiet is therefore re-phased by the stroke,
%    which is what makes the attack a coherent click, while a mode already
%    ringing loudly keeps the phase it had, which is also right. What is lost
%    is the intermodulation the real cascade would produce, since the energy
%    it moves into a mode arrives with whatever phase that mode already had.
%
% 2. Direct strike clicks. Modes above max_mode_hz are truncated; their
%    collective off-resonance response is the immediate local motion of the
%    plate under the tip. It is reinstated from the driving-point impedance of
%    an infinite plate, which for a plate is a pure resistance
%    Z = 8 sqrt(D rho t): the contact force divided by Z is the local
%    velocity, which radiates as a dipole over about half a flexural
%    wavelength, evaluated at the pulse's own centre frequency om* = pi / tau.
%    Scaled by direct_gain. This is the least principled term of the model,
%    and direct_gain = 0 turns it off.
%
% Returns p (1 x T, pascals at fsHi) and radinfo diagnostics.

c = P.speed_of_sound;
rho0 = P.air_density;
N = basis.n;
nsurf = length(basis.surfR);
nth = max(8, round(P.n_azimuth));
lx = P.listener_x_m;
ly = P.listener_y_m;
lz = P.listener_z_m;

% ---- surface element geometry, meridian samples x azimuths ----
th = reshape(((0:nth - 1) + 0.5) / nth * 2 * pi, 1, nth);
dth = 2 * pi / nth;
ct = cos(th);
st = sin(th);
ex = basis.surfR * ct;                       % nsurf x nth
ey = basis.surfR * st;
ez = basis.surfZ * ones(1, nth);
nx = (-basis.surfDz) * ct;
ny = (-basis.surfDz) * st;
nz = basis.surfDr * ones(1, nth);
dSa = (basis.surfR .* basis.surfLen) * (dth * ones(1, nth));

dx = lx - ex;
dy = ly - ey;
dz = lz - ez;
re = sqrt(dx.^2 + dy.^2 + dz.^2);
listenerInside = min(min(re)) < 0.02;
re = max(re, 0.01);
cosg = (nx .* dx + ny .* dy + nz .* dz) ./ re;   % nsurf x nth

cx = 0;
cy = 0;
cz = sum(basis.surfR .* basis.surfLen .* basis.surfZ) / sum(basis.surfR .* basis.surfLen);
rcen = sqrt((lx - cx)^2 + (ly - cy)^2 + (lz - cz)^2);

% ---- per-mode listener transfer ----
G = zeros(N, 1) + 0i;
for m = 1:N
  om = basis.omega(m);
  k = om / c;
  azi = cos(basis.ndia(m) * th);                          % 1 x nth
  amp = (basis.ldip(m) * dSa .* (basis.SW(:, m) * ones(1, nth))) .* (ones(nsurf, 1) * azi);
  ker = amp .* cosg .* (1i * k + 1 ./ re) .* exp(-1i * k * (re - rcen)) ./ re;
  G(m) = (1i * om * rho0 / (4 * pi)) * sum(sum(ker));
end
% the plate's real modes are the hammered hybrids, and the transfer is linear
% in the shape, so the mixing of cymbal_basis.m applies directly here
G = basis.mix * G;

% ---- oscillator bank driven by the cascade envelopes ----
T = floor(tframe(end) * fsHi);
p = zeros(1, T);
nf = length(tframe) - 1;
L = max(1, round((tframe(2) - tframe(1)) * fsHi));
omd = basis.omega .* sqrt(max(1e-12, 1 - basis.zeta.^2));
Eb = exp(1i * omd * ((0:L - 1) / fsHi));                  % N x L, hoisted
stepPh = exp(1i * omd * (L / fsHi));                     % per-frame advance
decLin = exp(-basis.zeta .* basis.omega * (L / fsHi));
ramp = (0:L - 1) / L;
sPh = ones(N, 1) + 0i;
aLin = zeros(N, 1);
pos = 1;
for f = 1:nf
  if pos > T
    break;
  end
  hit = find(kf == f);
  for q = 1:length(hit)
    cc = aLin .* sPh + DV(:, hit(q));
    aLin = abs(cc);
    sPh = cc ./ max(1e-300, aLin);
  end
  len = min(L, T - pos + 1);
  A0 = env(:, f);
  A1 = env(:, f + 1);
  c1 = (G .* A0 .* sPh).';
  c2 = (G .* (A1 - A0) .* sPh).';
  seg = real((c1 * Eb(:, 1:len)) + (c2 * Eb(:, 1:len)) .* ramp(1:len));
  p(pos:pos + len - 1) = seg;
  sPh = sPh .* stepPh;
  aLin = aLin .* decLin;
  pos = pos + L;
end

% ---- direct strike clicks ----
E = P.youngs_modulus;
nu = P.poisson_ratio;
rho = P.density;
cL = sqrt(E / (rho * (1 - nu^2)));
I = length(tstr);
ks = min(T, max(1, round(tstr * fsHi) + 1));
ang = P.strike_angle_deg * pi / 180;
directPeak = 0;
for i = 1:I
  [~, ~, thk, drds, dzds] = cymbal_at(basis, sstr(i));
  [rs, zs] = cymbal_point(basis, sstr(i));
  dvx = cos(ang) * dzds + sin(ang) * drds;    % stroke direction at theta = 0
  dvz = -cos(ang) * drds + sin(ang) * dzds;
  rsx = lx - rs;
  rsy = ly;
  rsz = lz - zs;
  rr = max(0.01, sqrt(rsx^2 + rsy^2 + rsz^2));
  cg = (rsx * dvx + rsz * dvz) / rr;
  nP = max(2, round(taustr(i) * fsHi));
  % the click carries the true path delay relative to the carrier that was
  % stripped from G, which is negative when the contact is nearer the
  % listener than the plate centroid
  k0 = max(1, ks(i) + round((rr - rcen) / c * fsHi));
  if k0 + nP > T
    continue;
  end
  D = E * thk^3 / (12 * (1 - nu^2));
  Zp = 8 * sqrt(D * rho * thk);
  omStar = pi / taustr(i);
  kap = thk / sqrt(12);
  kB = sqrt(omStar / (kap * cL));
  lamB = 2 * pi / kB;
  Aeff = pi * (lamB / 4)^2;                   % patch driven by the contact
  ldir = pi / max(kB, omStar / c);
  mu = Aeff * ldir / Zp;                      % dipole volume per unit force
  amp = P.direct_gain * fstr(i) * mu * rho0 * omStar^2 * cg / (4 * pi * rr * c);
  pulse = amp * sin(pi * (0:nP) / nP);
  p(k0:k0 + nP) = p(k0:k0 + nP) + pulse;
  if abs(amp) > directPeak
    directPeak = abs(amp);
  end
end

radinfo.gainAbs = abs(G);
radinfo.directPeakPa = directPeak;
radinfo.listenerInside = listenerInside;
radinfo.listenerDistM = rcen;
end