1function [p, spacingNeed, listenerInside] = woodblock_radiate(P, basis, Qdd, Xdd, fsHi, M)
2% WOODBLOCK_RADIATE Retarded free-field pressure at the listener point.
3% The top and bottom faces move together, so each strip is a local dipole;
4% the strip surface accelerations are computed in one matrix product and
5% reused for both faces with opposite sign (woodblock.md section 7).
6% spacingNeed is 0 when the element spacing is fine enough, otherwise the
7% recommended element count; listenerInside flags |r| < L/2.
9N = basis.n;
10T = size(Qdd, 2);
11L = P.length_m; W = P.width_m;
12halfH = 0.5 * P.block_height_m;
13patchA = W * L / M;
14c = P.speed_of_sound;
15rho0 = P.air_density;
17sElem = ((1:M).' - 0.5) / M;
18phiElem = zeros(M, N);
19for n = 1:N
20 phiElem(:, n) = interp1(basis.s, basis.shapes(:, n), sElem);
21end
22xElem = (sElem - 0.5) * L;
24lx = P.listener_x_m; ly = P.listener_y_m; lz = P.listener_z_m;
25dTop = sqrt((xElem - lx).^2 + ly^2 + (halfH - lz)^2);
26dBot = sqrt((xElem - lx).^2 + ly^2 + (-halfH - lz)^2);
27slitD = sqrt((P.slit_x_m - lx)^2 + (P.slit_y_m - ly)^2 + (P.slit_z_m - lz)^2);
29maxDelay = ceil(max([dTop; dBot; slitD]) / c * fsHi) + 2;
30out = zeros(1, T + maxDelay);
32wddAll = phiElem * Qdd; % M x T strip surface accelerations
33for j = 1:M
34 aT = patchA * rho0 / (4 * pi * dTop(j));
35 aB = -patchA * rho0 / (4 * pi * dBot(j));
36 kT = round(dTop(j) / c * fsHi) + 1;
37 kB = round(dBot(j) / c * fsHi) + 1;
38 out(kT:kT + T - 1) = out(kT:kT + T - 1) + aT * wddAll(j, :);
39 out(kB:kB + T - 1) = out(kB:kB + T - 1) + aB * wddAll(j, :);
40end
42% Slit monopole
43aS = P.slit_area_m2 * rho0 / (4 * pi * slitD);
44kS = round(slitD / c * fsHi) + 1;
45out(kS:kS + T - 1) = out(kS:kS + T - 1) + aS * Xdd;
47p = out(1:T);
49% Setup checks (check_setup in the reference)
50fMax = basis.f(N);
51if L / M > 0.25 * c / fMax
52 spacingNeed = ceil(L / (0.25 * c / fMax));
53else
54 spacingNeed = 0;
55end
56listenerInside = sqrt(lx^2 + ly^2 + lz^2) < 0.5 * L;
57end