aa03146Maraca shake simulator: MATLAB model and interactive web appmagland 1function [p, radinfo] = maraca_radiate(P, basis, timp, jimp, tauimp, fimp, ximp, fsHi)
2% MARACA_RADIATE Modal synthesis and free-field pressure at the listener.
3%
4% Two contributions:
5%
6% 1. Modal ringing. An impact with impulse J at surface point x_i gives every
7% mode (l, m) a velocity jump J S(omega tau) Y_lm(x_i) / m_l, where
8% S(x) = cos(x/2) / (1 - (x/pi)^2) is the spectrum of the half-sine contact
9% pulse (harder or faster grains excite higher modes). At mode frequency
10% the sphere radiates each (l, m) with the exact multipole transfer
11% -i rho0 c h_l(kr) / h_l'(ka); modes are narrow-band, so the transfer is
12% evaluated at the mode centre frequency. Because all m of a degree share
13% one frequency and one transfer, the addition theorem
14% sum_m Y_lm(x_i) Y_lm(x_r) = (2l+1)/(4pi) P_l(cos gamma)
15% collapses the m sum: only one damped resonator per retained mode is
16% synthesized, carrying listener-projected weights. The listener-direction
17% projection therefore lives in this stage, and moving the listener re-runs
18% only this stage. The carrier phase e^{-ikr} is stripped from h_l(kr) and
19% the propagation applied as a common group delay, so envelope and phase
20% are not double counted.
21%
22% 2. Direct impact clicks. Modal truncation at max_mode_hz discards the dense
23% high-l modes whose collective off-resonance response is the immediate
24% forced motion of the wall around the contact point. That broadband
25% component is reinstated with the classic point-driven infinite plate
26% result p = rho0 F(t - r/c) / (2 pi rho_s h r), scaled by direct_gain,
27% weighted (1 + cos gamma)/2 for shadowing of impacts on the far side, and
28% advanced/delayed by -R cos(gamma)/c for the position on the sphere. This
29% is the least principled term of the model: the plate formula ignores the
30% shell's curvature and closure.
31%
32% The synthesis is exact in time: between impacts every resonator evolves as
33% exp(lambda t), evaluated per inter-impact segment; nothing is integrated.
34%
35% Returns p (1 x T, pascals at fsHi) and radinfo diagnostics.
37c = P.speed_of_sound;
38rho0 = P.air_density;
39R = P.shell_radius_m;
40duration = P.stroke_count / P.stroke_rate_hz + P.tail_s;
41T = floor(duration * fsHi);
42N = basis.n;
43I = length(timp);
45lx = P.listener_x_m; ly = P.listener_y_m; lz = P.listener_z_m;
46rmag = sqrt(lx^2 + ly^2 + lz^2);
47listenerInside = rmag < R + 0.02;
48rmag = max(rmag, R + 0.005);
49xhat = [lx; ly; lz] / max(1e-9, sqrt(lx^2 + ly^2 + lz^2));
51% ---- listener geometry per impact: Legendre table P_l(cos gamma) ----
52cosg = xhat.' * ximp; % 1 x I
53lmax = max(basis.l);
54PL = zeros(lmax + 1, I);
55PL(1, :) = ones(1, I);
56if lmax >= 1
57 PL(2, :) = cosg;
58 for l = 1:lmax - 1
59 PL(l + 2, :) = ((2 * l + 1) * cosg .* PL(l + 1, :) - l * PL(l, :)) / (l + 1);
60 end
61end
63% ---- per-mode radiation transfer and complex pole ----
64G = zeros(N, 1) + 0i;
65lam = zeros(N, 1) + 0i;
66phv = zeros(N, 1) + 0i; % velocity phasor per unit jump
67for m = 1:N
68 l = basis.l(m);
69 om = basis.omega(m);
70 zt = basis.zeta(m);
71 k = om / c;
72 [ha, hpa] = maraca_hankel(l, k * R);
73 [hr, hpr] = maraca_hankel(l, k * rmag);
74 hred = hr(l + 1) * exp(1i * k * rmag); % strip carrier: no double delay
75 G(m) = -1i * rho0 * c * hred / hpa(l + 1);
76 omd = om * sqrt(max(1e-12, 1 - zt^2));
77 lam(m) = -zt * om + 1i * omd;
78 phv(m) = 1 + 1i * zt * om / omd; % v(t) = Re[phv e^{lam t}]
79end
81% ---- per-impact modal weights: A(n, i) velocity-jump phasors ----
82x = basis.omega * tauimp; % N x I
83den = 1 - (x / pi).^2;
84S = cos(x / 2) ./ den;
85sing = abs(den) < 1e-6;
86S(sing) = pi / 4; % limit at x -> pi
87Cl = ((2 * basis.l + 1) / (4 * pi)) .* (1 ./ basis.massKg); % N x 1
88A = (phv * jimp) .* S .* (Cl * ones(1, I)) .* PL(basis.l + 1, :);
90% ---- segment synthesis: exact damped-exponential evolution ----
91ks = min(T, max(1, round(timp * fsHi) + 1)); % impact sample indices, sorted
92p = zeros(1, T);
93s = zeros(N, 1) + 0i;
94Gt = G.';
95dt = 1 / fsHi;
96prevK = ks(1);
97i = 1;
98while i <= I
99 % accumulate all impacts landing on this sample
100 s = s + A(:, i);
101 while i < I && ks(i + 1) == ks(i)
102 i = i + 1;
103 s = s + A(:, i);
104 end
105 if i < I
106 nextK = ks(i + 1);
107 else
108 nextK = T + 1;
109 end
110 len = nextK - prevK;
111 if len > 0
112 seg = real(Gt * ((s * ones(1, len)) .* exp(lam * ((0:len - 1) * dt))));
113 p(prevK:nextK - 1) = p(prevK:nextK - 1) + seg;
114 s = s .* exp(lam * (len * dt));
115 end
116 prevK = nextK;
117 i = i + 1;
118end
120% ---- direct impact clicks (point-driven plate approximation) ----
121mpp = P.density * P.shell_thickness_m; % shell mass per unit area
122dirScale = P.direct_gain * rho0 / (2 * pi * mpp * rmag);
123directPeak = 0;
124for i = 1:I
125 nP = max(2, round(tauimp(i) * fsHi));
126 k0 = ks(i) + round(-R * cosg(i) / c * fsHi); % near side arrives earlier
127 if k0 < 1 || k0 + nP > T
128 continue;
129 end
130 amp = dirScale * fimp(i) * 0.5 * (1 + cosg(i));
131 pulse = amp * sin(pi * (0:nP) / nP);
132 p(k0:k0 + nP) = p(k0:k0 + nP) + pulse;
133 if amp > directPeak
134 directPeak = amp;
135 end
136end
138radinfo.gainAbs = abs(G);
139radinfo.directPeakPa = directPeak;
140radinfo.listenerInside = listenerInside;
141radinfo.listenerDistM = rmag;
142end