function [p, radinfo] = maraca_radiate(P, basis, timp, jimp, tauimp, fimp, ximp, fsHi)
% MARACA_RADIATE Modal synthesis and free-field pressure at the listener.
%
% Two contributions:
%
% 1. Modal ringing. An impact with impulse J at surface point x_i gives every
%    mode (l, m) a velocity jump J S(omega tau) Y_lm(x_i) / m_l, where
%    S(x) = cos(x/2) / (1 - (x/pi)^2) is the spectrum of the half-sine contact
%    pulse (harder or faster grains excite higher modes). At mode frequency
%    the sphere radiates each (l, m) with the exact multipole transfer
%    -i rho0 c h_l(kr) / h_l'(ka); modes are narrow-band, so the transfer is
%    evaluated at the mode centre frequency. Because all m of a degree share
%    one frequency and one transfer, the addition theorem
%      sum_m Y_lm(x_i) Y_lm(x_r) = (2l+1)/(4pi) P_l(cos gamma)
%    collapses the m sum: only one damped resonator per retained mode is
%    synthesized, carrying listener-projected weights. The listener-direction
%    projection therefore lives in this stage, and moving the listener re-runs
%    only this stage. The carrier phase e^{-ikr} is stripped from h_l(kr) and
%    the propagation applied as a common group delay, so envelope and phase
%    are not double counted.
%
% 2. Direct impact clicks. Modal truncation at max_mode_hz discards the dense
%    high-l modes whose collective off-resonance response is the immediate
%    forced motion of the wall around the contact point. That broadband
%    component is reinstated with the classic point-driven infinite plate
%    result p = rho0 F(t - r/c) / (2 pi rho_s h r), scaled by direct_gain,
%    weighted (1 + cos gamma)/2 for shadowing of impacts on the far side, and
%    advanced/delayed by -R cos(gamma)/c for the position on the sphere. This
%    is the least principled term of the model: the plate formula ignores the
%    shell's curvature and closure.
%
% The synthesis is exact in time: between impacts every resonator evolves as
% exp(lambda t), evaluated per inter-impact segment; nothing is integrated.
%
% Returns p (1 x T, pascals at fsHi) and radinfo diagnostics.

c = P.speed_of_sound;
rho0 = P.air_density;
R = P.shell_radius_m;
duration = P.stroke_count / P.stroke_rate_hz + P.tail_s;
T = floor(duration * fsHi);
N = basis.n;
I = length(timp);

lx = P.listener_x_m; ly = P.listener_y_m; lz = P.listener_z_m;
rmag = sqrt(lx^2 + ly^2 + lz^2);
listenerInside = rmag < R + 0.02;
rmag = max(rmag, R + 0.005);
xhat = [lx; ly; lz] / max(1e-9, sqrt(lx^2 + ly^2 + lz^2));

% ---- listener geometry per impact: Legendre table P_l(cos gamma) ----
cosg = xhat.' * ximp;                          % 1 x I
lmax = max(basis.l);
PL = zeros(lmax + 1, I);
PL(1, :) = ones(1, I);
if lmax >= 1
  PL(2, :) = cosg;
  for l = 1:lmax - 1
    PL(l + 2, :) = ((2 * l + 1) * cosg .* PL(l + 1, :) - l * PL(l, :)) / (l + 1);
  end
end

% ---- per-mode radiation transfer and complex pole ----
G = zeros(N, 1) + 0i;
lam = zeros(N, 1) + 0i;
phv = zeros(N, 1) + 0i;                        % velocity phasor per unit jump
for m = 1:N
  l = basis.l(m);
  om = basis.omega(m);
  zt = basis.zeta(m);
  k = om / c;
  [ha, hpa] = maraca_hankel(l, k * R);
  [hr, hpr] = maraca_hankel(l, k * rmag);
  hred = hr(l + 1) * exp(1i * k * rmag);       % strip carrier: no double delay
  G(m) = -1i * rho0 * c * hred / hpa(l + 1);
  omd = om * sqrt(max(1e-12, 1 - zt^2));
  lam(m) = -zt * om + 1i * omd;
  phv(m) = 1 + 1i * zt * om / omd;             % v(t) = Re[phv e^{lam t}]
end

% ---- per-impact modal weights: A(n, i) velocity-jump phasors ----
x = basis.omega * tauimp;                      % N x I
den = 1 - (x / pi).^2;
S = cos(x / 2) ./ den;
sing = abs(den) < 1e-6;
S(sing) = pi / 4;                              % limit at x -> pi
Cl = ((2 * basis.l + 1) / (4 * pi)) .* (1 ./ basis.massKg);   % N x 1
A = (phv * jimp) .* S .* (Cl * ones(1, I)) .* PL(basis.l + 1, :);

% ---- segment synthesis: exact damped-exponential evolution ----
ks = min(T, max(1, round(timp * fsHi) + 1));   % impact sample indices, sorted
p = zeros(1, T);
s = zeros(N, 1) + 0i;
Gt = G.';
dt = 1 / fsHi;
prevK = ks(1);
i = 1;
while i <= I
  % accumulate all impacts landing on this sample
  s = s + A(:, i);
  while i < I && ks(i + 1) == ks(i)
    i = i + 1;
    s = s + A(:, i);
  end
  if i < I
    nextK = ks(i + 1);
  else
    nextK = T + 1;
  end
  len = nextK - prevK;
  if len > 0
    seg = real(Gt * ((s * ones(1, len)) .* exp(lam * ((0:len - 1) * dt))));
    p(prevK:nextK - 1) = p(prevK:nextK - 1) + seg;
    s = s .* exp(lam * (len * dt));
  end
  prevK = nextK;
  i = i + 1;
end

% ---- direct impact clicks (point-driven plate approximation) ----
mpp = P.density * P.shell_thickness_m;         % shell mass per unit area
dirScale = P.direct_gain * rho0 / (2 * pi * mpp * rmag);
directPeak = 0;
for i = 1:I
  nP = max(2, round(tauimp(i) * fsHi));
  k0 = ks(i) + round(-R * cosg(i) / c * fsHi); % near side arrives earlier
  if k0 < 1 || k0 + nP > T
    continue;
  end
  amp = dirScale * fimp(i) * 0.5 * (1 + cosg(i));
  pulse = amp * sin(pi * (0:nP) / nP);
  p(k0:k0 + nP) = p(k0:k0 + nP) + pulse;
  if amp > directPeak
    directPeak = amp;
  end
end

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