1function [p, radinfo] = cymbal_radiate(P, basis, env, tframe, DV, kf, tstr, taustr, fstr, sstr, fsHi)
2% CYMBAL_RADIATE Free-field pressure at the listener from the modal envelopes.
3%
4% Two contributions:
5%
6% 1. The mode bank. Each retained mode is a normal-dipole sheet on the plate:
7% a surface element of area dS carries a dipole of volume dS * ldip along
8% the local normal (ldip from cymbal_basis.m), and the element dipoles are
9% summed coherently with their relative retardations at the mode's centre
10% frequency,
11% G = sum_e (i om rho0 / 4 pi) V_e (u_e . rhat_e) (i k + 1/r_e)
12% exp(-i k (r_e - rcen)) / r_e,
13% over n_azimuth azimuths and two sample points per meridian element. The
14% common carrier exp(-i k rcen) is stripped, so relative phase across the
15% plate is kept while the overall propagation delay is not double counted.
16% Unlike the far-field form used for the radiated power, this keeps the
17% 1/r^2 near-field term and the true distances, which matters: a listener
18% half a metre from a 40 cm plate is nowhere near the far field. Only the
19% normal component of the surface motion radiates; the in-plane component
20% of a thin plate displaces almost no air.
21%
22% Because the azimuthal shape is exactly cos(n theta), the sum over
23% azimuths carries the multipole cancellation of the mode: on the axis of
24% symmetry every mode with n >= 1 cancels exactly, and off axis the high-n
25% modes are weak until their wavelength is short. Moving the listener
26% re-runs only this stage.
27%
28% The synthesis itself is an oscillator bank rather than the closed-form
29% damped exponentials of the sibling repositories: the cascade makes every
30% amplitude time-varying, so each mode is a fixed-frequency oscillator
31% whose envelope is the cascade solution, interpolated linearly across each
32% control frame. Within a frame the two matrix products below are exact,
33% and the per-frame phase advance is a unit-modulus multiply, so nothing
34% drifts.
35%
36% Amplitude comes from the cascade; PHASE comes from a linear impulse
37% response carried alongside it. Each oscillator keeps a unit phasor and the
38% magnitude the mode would have had with no cascade at all; a contact adds
39% its signed modal velocity jump to that complex state and the phasor is
40% renormalized. A mode that was quiet is therefore re-phased by the stroke,
41% which is what makes the attack a coherent click, while a mode already
42% ringing loudly keeps the phase it had, which is also right. What is lost
43% is the intermodulation the real cascade would produce, since the energy
44% it moves into a mode arrives with whatever phase that mode already had.
45%
46% 2. Direct strike clicks. Modes above max_mode_hz are truncated; their
47% collective off-resonance response is the immediate local motion of the
48% plate under the tip. It is reinstated from the driving-point impedance of
49% an infinite plate, which for a plate is a pure resistance
50% Z = 8 sqrt(D rho t): the contact force divided by Z is the local
51% velocity, which radiates as a dipole over about half a flexural
52% wavelength, evaluated at the pulse's own centre frequency om* = pi / tau.
53% Scaled by direct_gain. This is the least principled term of the model,
54% and direct_gain = 0 turns it off.
55%
56% Returns p (1 x T, pascals at fsHi) and radinfo diagnostics.
58c = P.speed_of_sound;
59rho0 = P.air_density;
60N = basis.n;
61nsurf = length(basis.surfR);
62nth = max(8, round(P.n_azimuth));
63lx = P.listener_x_m;
64ly = P.listener_y_m;
65lz = P.listener_z_m;
67% ---- surface element geometry, meridian samples x azimuths ----
68th = reshape(((0:nth - 1) + 0.5) / nth * 2 * pi, 1, nth);
69dth = 2 * pi / nth;
70ct = cos(th);
71st = sin(th);
72ex = basis.surfR * ct; % nsurf x nth
73ey = basis.surfR * st;
74ez = basis.surfZ * ones(1, nth);
75nx = (-basis.surfDz) * ct;
76ny = (-basis.surfDz) * st;
77nz = basis.surfDr * ones(1, nth);
78dSa = (basis.surfR .* basis.surfLen) * (dth * ones(1, nth));
80dx = lx - ex;
81dy = ly - ey;
82dz = lz - ez;
83re = sqrt(dx.^2 + dy.^2 + dz.^2);
84listenerInside = min(min(re)) < 0.02;
85re = max(re, 0.01);
86cosg = (nx .* dx + ny .* dy + nz .* dz) ./ re; % nsurf x nth
88cx = 0;
89cy = 0;
90cz = sum(basis.surfR .* basis.surfLen .* basis.surfZ) / sum(basis.surfR .* basis.surfLen);
91rcen = sqrt((lx - cx)^2 + (ly - cy)^2 + (lz - cz)^2);
93% ---- per-mode listener transfer ----
94G = zeros(N, 1) + 0i;
95for m = 1:N
96 om = basis.omega(m);
97 k = om / c;
98 azi = cos(basis.ndia(m) * th); % 1 x nth
99 amp = (basis.ldip(m) * dSa .* (basis.SW(:, m) * ones(1, nth))) .* (ones(nsurf, 1) * azi);
100 ker = amp .* cosg .* (1i * k + 1 ./ re) .* exp(-1i * k * (re - rcen)) ./ re;
101 G(m) = (1i * om * rho0 / (4 * pi)) * sum(sum(ker));
102end
104% ---- oscillator bank driven by the cascade envelopes ----
105T = floor(tframe(end) * fsHi);
106p = zeros(1, T);
107nf = length(tframe) - 1;
108L = max(1, round((tframe(2) - tframe(1)) * fsHi));
109omd = basis.omega .* sqrt(max(1e-12, 1 - basis.zeta.^2));
110Eb = exp(1i * omd * ((0:L - 1) / fsHi)); % N x L, hoisted
111stepPh = exp(1i * omd * (L / fsHi)); % per-frame advance
112decLin = exp(-basis.zeta .* basis.omega * (L / fsHi));
113ramp = (0:L - 1) / L;
114sPh = ones(N, 1) + 0i;
115aLin = zeros(N, 1);
116pos = 1;
117for f = 1:nf
118 if pos > T
119 break;
120 end
121 hit = find(kf == f);
122 for q = 1:length(hit)
123 cc = aLin .* sPh + DV(:, hit(q));
124 aLin = abs(cc);
125 sPh = cc ./ max(1e-300, aLin);
126 end
127 len = min(L, T - pos + 1);
128 A0 = env(:, f);
129 A1 = env(:, f + 1);
130 c1 = (G .* A0 .* sPh).';
131 c2 = (G .* (A1 - A0) .* sPh).';
132 seg = real((c1 * Eb(:, 1:len)) + (c2 * Eb(:, 1:len)) .* ramp(1:len));
133 p(pos:pos + len - 1) = seg;
134 sPh = sPh .* stepPh;
135 aLin = aLin .* decLin;
136 pos = pos + L;
137end
139% ---- direct strike clicks ----
140E = P.youngs_modulus;
141nu = P.poisson_ratio;
142rho = P.density;
143cL = sqrt(E / (rho * (1 - nu^2)));
144I = length(tstr);
145ks = min(T, max(1, round(tstr * fsHi) + 1));
146ang = P.strike_angle_deg * pi / 180;
147directPeak = 0;
148for i = 1:I
149 [~, ~, thk, drds, dzds] = cymbal_at(basis, sstr(i));
150 [rs, zs] = cymbal_point(basis, sstr(i));
151 dvx = cos(ang) * dzds + sin(ang) * drds; % stroke direction at theta = 0
152 dvz = -cos(ang) * drds + sin(ang) * dzds;
153 rsx = lx - rs;
154 rsy = ly;
155 rsz = lz - zs;
156 rr = max(0.01, sqrt(rsx^2 + rsy^2 + rsz^2));
157 cg = (rsx * dvx + rsz * dvz) / rr;
158 nP = max(2, round(taustr(i) * fsHi));
159 % the click carries the true path delay relative to the carrier that was
160 % stripped from G, which is negative when the contact is nearer the
161 % listener than the plate centroid
162 k0 = max(1, ks(i) + round((rr - rcen) / c * fsHi));
163 if k0 + nP > T
164 continue;
165 end
166 D = E * thk^3 / (12 * (1 - nu^2));
167 Zp = 8 * sqrt(D * rho * thk);
168 omStar = pi / taustr(i);
169 kap = thk / sqrt(12);
170 kB = sqrt(omStar / (kap * cL));
171 lamB = 2 * pi / kB;
172 Aeff = pi * (lamB / 4)^2; % patch driven by the contact
173 ldir = pi / max(kB, omStar / c);
174 mu = Aeff * ldir / Zp; % dipole volume per unit force
175 amp = P.direct_gain * fstr(i) * mu * rho0 * omStar^2 * cg / (4 * pi * rr * c);
176 pulse = amp * sin(pi * (0:nP) / nP);
177 p(k0:k0 + nP) = p(k0:k0 + nP) + pulse;
178 if abs(amp) > directPeak
179 directPeak = abs(amp);
180 end
181end
183radinfo.gainAbs = abs(G);
184radinfo.directPeakPa = directPeak;
185radinfo.listenerInside = listenerInside;
186radinfo.listenerDistM = rcen;
187end