function [p, spacingNeed, listenerInside] = woodblock_radiate(P, basis, Qdd, Xdd, fsHi, M)
% WOODBLOCK_RADIATE Retarded free-field pressure at the listener point.
% The top and bottom faces move together, so each strip is a local dipole;
% the strip surface accelerations are computed in one matrix product and
% reused for both faces with opposite sign (woodblock.md section 7).
% spacingNeed is 0 when the element spacing is fine enough, otherwise the
% recommended element count; listenerInside flags |r| < L/2.

N = basis.n;
T = size(Qdd, 2);
L = P.length_m; W = P.width_m;
halfH = 0.5 * P.block_height_m;
patchA = W * L / M;
c = P.speed_of_sound;
rho0 = P.air_density;

sElem = ((1:M).' - 0.5) / M;
phiElem = zeros(M, N);
for n = 1:N
  phiElem(:, n) = interp1(basis.s, basis.shapes(:, n), sElem);
end
xElem = (sElem - 0.5) * L;

lx = P.listener_x_m; ly = P.listener_y_m; lz = P.listener_z_m;
dTop = sqrt((xElem - lx).^2 + ly^2 + (halfH - lz)^2);
dBot = sqrt((xElem - lx).^2 + ly^2 + (-halfH - lz)^2);
slitD = sqrt((P.slit_x_m - lx)^2 + (P.slit_y_m - ly)^2 + (P.slit_z_m - lz)^2);

maxDelay = ceil(max([dTop; dBot; slitD]) / c * fsHi) + 2;
out = zeros(1, T + maxDelay);

wddAll = phiElem * Qdd;                 % M x T strip surface accelerations
for j = 1:M
  aT = patchA * rho0 / (4 * pi * dTop(j));
  aB = -patchA * rho0 / (4 * pi * dBot(j));
  kT = round(dTop(j) / c * fsHi) + 1;
  kB = round(dBot(j) / c * fsHi) + 1;
  out(kT:kT + T - 1) = out(kT:kT + T - 1) + aT * wddAll(j, :);
  out(kB:kB + T - 1) = out(kB:kB + T - 1) + aB * wddAll(j, :);
end

% Slit monopole
aS = P.slit_area_m2 * rho0 / (4 * pi * slitD);
kS = round(slitD / c * fsHi) + 1;
out(kS:kS + T - 1) = out(kS:kS + T - 1) + aS * Xdd;

p = out(1:T);

% Setup checks (check_setup in the reference)
fMax = basis.f(N);
if L / M > 0.25 * c / fMax
  spacingNeed = ceil(L / (0.25 * c / fMax));
else
  spacingNeed = 0;
end
listenerInside = sqrt(lx^2 + ly^2 + lz^2) < 0.5 * L;
end