1function [p, radinfo] = triangle_radiate(P, basis, tstr, vstr, jstr, taustr, fstr, sstr, fsHi)
2% TRIANGLE_RADIATE Modal synthesis and free-field pressure at the listener.
3%
4% Two contributions:
5%
6% 1. Modal ringing. A strike with impulse J at arclength position s gives
7% every mode a modal-velocity jump J S(omega tau) (phi(s) . dhat), where
8% S(x) = cos(x/2) / (1 - (x/pi)^2) is the spectrum of the half-sine
9% contact pulse and dhat the strike direction (strike_angle_deg between
10% the leg's in-plane normal and the plane normal, so one strike excites
11% both mode families). The rod radiates through its transverse motion:
12% each element is a point dipole of moment Vdip u_perp qdot (Vdip from
13% triangle_basis.m), and per mode the element dipoles are summed
14% coherently with their relative retardations,
15% G = sum_e (i om rho0 / 4 pi) Vdip_e (u_e . rhat_e)
16% (i k + 1/r_e) exp(-i k (r_e - rcen)) / r_e,
17% evaluated at the mode's centre frequency (modes are narrow-band). The
18% common carrier exp(-i k rcen) is stripped, so relative phase across the
19% triangle is kept while the overall propagation delay is not double
20% counted. This sum is where the interference across the frame lives:
21% moving the listener re-runs only this stage.
22%
23% 2. Direct strike clicks. Modes above max_mode_hz are truncated; their
24% collective off-resonance response is the immediate local motion of the
25% rod around the contact. It is reinstated from the classic driving-point
26% impedance of an infinite Euler beam, Z = 2 rho A c_B (1 + i) with
27% c_B = sqrt(om kappa cL), evaluated at the pulse's own centre frequency
28% om* = pi / tau: the contact force divided by |Z| is the local velocity,
29% which radiates as a dipole of length half a flexural wavelength along
30% the strike direction. Scaled by direct_gain. This is the least
31% principled term of the model, and direct_gain = 0 turns it off.
32%
33% The synthesis is exact in time: between strikes every mode evolves as
34% exp(lambda t), evaluated in chunks over each inter-strike segment (the
35% chunking only bounds memory; nothing is integrated).
36%
37% Returns p (1 x T, pascals at fsHi) and radinfo diagnostics.
39c = P.speed_of_sound;
40rho0 = P.air_density;
41K = max(1, round(P.strike_count));
42duration = (K - 1) / max(1e-3, P.strike_rate_hz) + P.tail_s;
43T = floor(duration * fsHi);
44N = basis.n;
45I = length(tstr);
46nel = length(basis.elemLen);
47nn = length(basis.sfrac);
49lx = P.listener_x_m;
50ly = P.listener_y_m;
51lz = P.listener_z_m;
53% ---- element-to-listener geometry ----
54rex = lx - basis.elemX; % nel x 1
55rey = ly - basis.elemY;
56rez = lz * ones(nel, 1);
57re = sqrt(rex.^2 + rey.^2 + rez.^2);
58listenerInside = min(re) < 0.02;
59re = max(re, 0.01);
60rhx = rex ./ re;
61rhy = rey ./ re;
62rhz = rez ./ re;
63cx = sum(basis.elemX .* basis.elemLen) / sum(basis.elemLen);
64cy = sum(basis.elemY .* basis.elemLen) / sum(basis.elemLen);
65rcen = sqrt((lx - cx)^2 + (ly - cy)^2 + lz^2);
67% ---- per-mode radiation transfer and complex pole ----
68a = P.rod_diameter_m / 2;
69Vdip = 2 * pi * a^2 * basis.elemLen; % nel x 1
70G = zeros(N, 1) + 0i;
71lam = zeros(N, 1) + 0i;
72phv = zeros(N, 1) + 0i; % modal-velocity phasor per unit jump
73for m = 1:N
74 om = basis.omega(m);
75 zt = basis.zeta(m);
76 k = om / c;
77 proj = basis.EUpx(:, m) .* rhx + basis.EUpy(:, m) .* rhy + basis.EUpz(:, m) .* rhz;
78 ph = exp(-1i * k * (re - rcen));
79 G(m) = (1i * om * rho0 / (4 * pi)) * sum(Vdip .* proj .* (1i * k + 1 ./ re) .* ph ./ re);
80 omd = om * sqrt(max(1e-12, 1 - zt^2));
81 lam(m) = -zt * om + 1i * omd;
82 phv(m) = 1 + 1i * zt * om / omd; % qdot(t) = Re[phv e^{lam t}] per unit jump
83end
85% ---- per-strike modal weights: A(n, i) modal-velocity-jump phasors ----
86A = zeros(N, I) + 0i;
87dirx = zeros(1, I);
88diry = zeros(1, I);
89dirz = zeros(1, I);
90hitx = zeros(1, I);
91hity = zeros(1, I);
92ang = P.strike_angle_deg * pi / 180;
93for i = 1:I
94 s = min(0.999, max(0.001, sstr(i)));
95 idx = sum(basis.sfrac <= s);
96 idx = min(max(idx, 1), nn - 1);
97 w = (s - basis.sfrac(idx)) / (basis.sfrac(idx + 1) - basis.sfrac(idx));
98 uxh = basis.NUx(idx, :) + w * (basis.NUx(idx + 1, :) - basis.NUx(idx, :));
99 uyh = basis.NUy(idx, :) + w * (basis.NUy(idx + 1, :) - basis.NUy(idx, :));
100 uzh = basis.NUz(idx, :) + w * (basis.NUz(idx + 1, :) - basis.NUz(idx, :));
101 % strike direction: the leg's inward in-plane normal, tilted out of plane
102 ex = basis.elemTanX(idx);
103 ey = basis.elemTanY(idx);
104 dirx(i) = cos(ang) * (-ey);
105 diry(i) = cos(ang) * ex;
106 dirz(i) = sin(ang);
107 hitx(i) = basis.nodeX(idx) + w * (basis.nodeX(idx + 1) - basis.nodeX(idx));
108 hity(i) = basis.nodeY(idx) + w * (basis.nodeY(idx + 1) - basis.nodeY(idx));
109 proj = uxh * dirx(i) + uyh * diry(i) + uzh * dirz(i); % 1 x N
110 x = basis.omega * taustr(i); % N x 1
111 den = 1 - (x / pi).^2;
112 S = cos(x / 2) ./ den;
113 sing = abs(den) < 1e-6;
114 S(sing) = pi / 4; % limit at x -> pi
115 A(:, i) = phv .* (jstr(i) * S) .* reshape(proj, N, 1);
116end
118% ---- segment synthesis: exact damped-exponential evolution, chunked ----
119ks = min(T, max(1, round(tstr * fsHi) + 1));
120p = zeros(1, T);
121s = zeros(N, 1) + 0i;
122Gt = G.';
123dt = 1 / fsHi;
124CH = 32768; % chunk length: bounds the exp matrix
125prevK = ks(1);
126i = 1;
127while i <= I
128 s = s + A(:, i);
129 while i < I && ks(i + 1) == ks(i)
130 i = i + 1;
131 s = s + A(:, i);
132 end
133 if i < I
134 nextK = ks(i + 1);
135 else
136 nextK = T + 1;
137 end
138 pos = prevK;
139 while pos < nextK
140 len = min(CH, nextK - pos);
141 seg = real(Gt * ((s * ones(1, len)) .* exp(lam * ((0:len - 1) * dt))));
142 p(pos:pos + len - 1) = p(pos:pos + len - 1) + seg;
143 s = s .* exp(lam * (len * dt));
144 pos = pos + len;
145 end
146 prevK = nextK;
147 i = i + 1;
148end
150% ---- direct strike clicks (driving-point velocity radiating as a dipole) ----
151Aarea = pi * (2 * a)^2 / 4;
152kapcL = (a / 2) * sqrt(P.youngs_modulus / P.density); % kappa * cL
153directPeak = 0;
154for i = 1:I
155 rsx = lx - hitx(i);
156 rsy = ly - hity(i);
157 rs = max(0.01, sqrt(rsx^2 + rsy^2 + lz^2));
158 cosg = (rsx * dirx(i) + rsy * diry(i) + lz * dirz(i)) / rs;
159 nP = max(2, round(taustr(i) * fsHi));
160 k0 = ks(i) + round((rs - rcen) / c * fsHi);
161 if k0 < 1 || k0 + nP > T
162 continue;
163 end
164 omStar = pi / taustr(i);
165 cB = sqrt(omStar * kapcL);
166 Zb = 2 * sqrt(2) * P.density * Aarea * cB; % |2 rho A cB (1 + i)|
167 lamB = 2 * pi * cB / omStar; % flexural wavelength at om*
168 mu = 2 * pi * a^2 * (lamB / 2) / Zb; % dipole moment per unit force
169 amp = P.direct_gain * fstr(i) * mu * rho0 * omStar^2 * cosg / (4 * pi * rs * c);
170 pulse = amp * sin(pi * (0:nP) / nP);
171 p(k0:k0 + nP) = p(k0:k0 + nP) + pulse;
172 if abs(amp) > directPeak
173 directPeak = abs(amp);
174 end
175end
177radinfo.gainAbs = abs(G);
178radinfo.directPeakPa = directPeak;
179radinfo.listenerInside = listenerInside;
180radinfo.listenerDistM = rcen;
181end