function [p, radinfo] = triangle_radiate(P, basis, tstr, vstr, jstr, taustr, fstr, sstr, fsHi)
% TRIANGLE_RADIATE Modal synthesis and free-field pressure at the listener.
%
% Two contributions:
%
% 1. Modal ringing. A strike with impulse J at arclength position s gives
%    every mode a modal-velocity jump J S(omega tau) (phi(s) . dhat), where
%    S(x) = cos(x/2) / (1 - (x/pi)^2) is the spectrum of the half-sine
%    contact pulse and dhat the strike direction (strike_angle_deg between
%    the leg's in-plane normal and the plane normal, so one strike excites
%    both mode families). The rod radiates through its transverse motion:
%    each element is a point dipole of moment Vdip u_perp qdot (Vdip from
%    triangle_basis.m), and per mode the element dipoles are summed
%    coherently with their relative retardations,
%      G = sum_e (i om rho0 / 4 pi) Vdip_e (u_e . rhat_e)
%            (i k + 1/r_e) exp(-i k (r_e - rcen)) / r_e,
%    evaluated at the mode's centre frequency (modes are narrow-band). The
%    common carrier exp(-i k rcen) is stripped, so relative phase across the
%    triangle is kept while the overall propagation delay is not double
%    counted. This sum is where the interference across the frame lives:
%    moving the listener re-runs only this stage.
%
% 2. Direct strike clicks. Modes above max_mode_hz are truncated; their
%    collective off-resonance response is the immediate local motion of the
%    rod around the contact. It is reinstated from the classic driving-point
%    impedance of an infinite Euler beam, Z = 2 rho A c_B (1 + i) with
%    c_B = sqrt(om kappa cL), evaluated at the pulse's own centre frequency
%    om* = pi / tau: the contact force divided by |Z| is the local velocity,
%    which radiates as a dipole of length half a flexural wavelength along
%    the strike direction. Scaled by direct_gain. This is the least
%    principled term of the model, and direct_gain = 0 turns it off.
%
% The synthesis is exact in time: between strikes every mode evolves as
% exp(lambda t), evaluated in chunks over each inter-strike segment (the
% chunking only bounds memory; nothing is integrated).
%
% Returns p (1 x T, pascals at fsHi) and radinfo diagnostics.

c = P.speed_of_sound;
rho0 = P.air_density;
K = max(1, round(P.strike_count));
duration = (K - 1) / max(1e-3, P.strike_rate_hz) + P.tail_s;
T = floor(duration * fsHi);
N = basis.n;
I = length(tstr);
nel = length(basis.elemLen);
nn = length(basis.sfrac);

lx = P.listener_x_m;
ly = P.listener_y_m;
lz = P.listener_z_m;

% ---- element-to-listener geometry ----
rex = lx - basis.elemX;              % nel x 1
rey = ly - basis.elemY;
rez = lz * ones(nel, 1);
re = sqrt(rex.^2 + rey.^2 + rez.^2);
listenerInside = min(re) < 0.02;
re = max(re, 0.01);
rhx = rex ./ re;
rhy = rey ./ re;
rhz = rez ./ re;
cx = sum(basis.elemX .* basis.elemLen) / sum(basis.elemLen);
cy = sum(basis.elemY .* basis.elemLen) / sum(basis.elemLen);
rcen = sqrt((lx - cx)^2 + (ly - cy)^2 + lz^2);

% ---- per-mode radiation transfer and complex pole ----
a = P.rod_diameter_m / 2;
Vdip = 2 * pi * a^2 * basis.elemLen; % nel x 1
G = zeros(N, 1) + 0i;
lam = zeros(N, 1) + 0i;
phv = zeros(N, 1) + 0i;              % modal-velocity phasor per unit jump
for m = 1:N
  om = basis.omega(m);
  zt = basis.zeta(m);
  k = om / c;
  proj = basis.EUpx(:, m) .* rhx + basis.EUpy(:, m) .* rhy + basis.EUpz(:, m) .* rhz;
  ph = exp(-1i * k * (re - rcen));
  G(m) = (1i * om * rho0 / (4 * pi)) * sum(Vdip .* proj .* (1i * k + 1 ./ re) .* ph ./ re);
  omd = om * sqrt(max(1e-12, 1 - zt^2));
  lam(m) = -zt * om + 1i * omd;
  phv(m) = 1 + 1i * zt * om / omd;   % qdot(t) = Re[phv e^{lam t}] per unit jump
end

% ---- per-strike modal weights: A(n, i) modal-velocity-jump phasors ----
A = zeros(N, I) + 0i;
dirx = zeros(1, I);
diry = zeros(1, I);
dirz = zeros(1, I);
hitx = zeros(1, I);
hity = zeros(1, I);
ang = P.strike_angle_deg * pi / 180;
for i = 1:I
  s = min(0.999, max(0.001, sstr(i)));
  idx = sum(basis.sfrac <= s);
  idx = min(max(idx, 1), nn - 1);
  w = (s - basis.sfrac(idx)) / (basis.sfrac(idx + 1) - basis.sfrac(idx));
  uxh = basis.NUx(idx, :) + w * (basis.NUx(idx + 1, :) - basis.NUx(idx, :));
  uyh = basis.NUy(idx, :) + w * (basis.NUy(idx + 1, :) - basis.NUy(idx, :));
  uzh = basis.NUz(idx, :) + w * (basis.NUz(idx + 1, :) - basis.NUz(idx, :));
  % strike direction: the leg's inward in-plane normal, tilted out of plane
  ex = basis.elemTanX(idx);
  ey = basis.elemTanY(idx);
  dirx(i) = cos(ang) * (-ey);
  diry(i) = cos(ang) * ex;
  dirz(i) = sin(ang);
  hitx(i) = basis.nodeX(idx) + w * (basis.nodeX(idx + 1) - basis.nodeX(idx));
  hity(i) = basis.nodeY(idx) + w * (basis.nodeY(idx + 1) - basis.nodeY(idx));
  proj = uxh * dirx(i) + uyh * diry(i) + uzh * dirz(i);   % 1 x N
  x = basis.omega * taustr(i);       % N x 1
  den = 1 - (x / pi).^2;
  S = cos(x / 2) ./ den;
  sing = abs(den) < 1e-6;
  S(sing) = pi / 4;                  % limit at x -> pi
  A(:, i) = phv .* (jstr(i) * S) .* reshape(proj, N, 1);
end

% ---- segment synthesis: exact damped-exponential evolution, chunked ----
ks = min(T, max(1, round(tstr * fsHi) + 1));
p = zeros(1, T);
s = zeros(N, 1) + 0i;
Gt = G.';
dt = 1 / fsHi;
CH = 32768;                          % chunk length: bounds the exp matrix
prevK = ks(1);
i = 1;
while i <= I
  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
  pos = prevK;
  while pos < nextK
    len = min(CH, nextK - pos);
    seg = real(Gt * ((s * ones(1, len)) .* exp(lam * ((0:len - 1) * dt))));
    p(pos:pos + len - 1) = p(pos:pos + len - 1) + seg;
    s = s .* exp(lam * (len * dt));
    pos = pos + len;
  end
  prevK = nextK;
  i = i + 1;
end

% ---- direct strike clicks (driving-point velocity radiating as a dipole) ----
Aarea = pi * (2 * a)^2 / 4;
kapcL = (a / 2) * sqrt(P.youngs_modulus / P.density);   % kappa * cL
directPeak = 0;
for i = 1:I
  rsx = lx - hitx(i);
  rsy = ly - hity(i);
  rs = max(0.01, sqrt(rsx^2 + rsy^2 + lz^2));
  cosg = (rsx * dirx(i) + rsy * diry(i) + lz * dirz(i)) / rs;
  nP = max(2, round(taustr(i) * fsHi));
  k0 = ks(i) + round((rs - rcen) / c * fsHi);
  if k0 < 1 || k0 + nP > T
    continue;
  end
  omStar = pi / taustr(i);
  cB = sqrt(omStar * kapcL);
  Zb = 2 * sqrt(2) * P.density * Aarea * cB;    % |2 rho A cB (1 + i)|
  lamB = 2 * pi * cB / omStar;                  % flexural wavelength at om*
  mu = 2 * pi * a^2 * (lamB / 2) / Zb;          % dipole moment per unit force
  amp = P.direct_gain * fstr(i) * mu * rho0 * omStar^2 * cosg / (4 * pi * rs * 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