function [p, radinfo] = recorder_radiate(P, basis, Zr, Zi, Qin, fsHi)
% RECORDER_RADIATE Free-field pressure at the listener for one note.
%
% Two contributions:
%
% 1. Modal aperture radiation. The resonant part of each modal flow is the
%    narrowband analytic signal a_n(t) = -2 c_n rho_n z_n(t) (see
%    recorder_jet.m), and basis.apW holds each aperture's complex flow per
%    unit loop flow at that mode's centre frequency. Every radiating
%    aperture (the window, each open tone hole, the foot) is a point
%    monopole, so per mode the apertures sum coherently with their relative
%    retardations,
%      G_n = sum_ap (i omega_n rho0 / (4 pi r_ap)) W_ap,n exp(-i k (r_ap - rcen)),
%    and p(t) = sum_n Re(G_n a_n(t)). The common carrier exp(-i k rcen) is
%    stripped: relative phase between the window, the open holes, and the
%    foot is kept (this is where the directivity of a cross-fingered note
%    lives), the overall propagation delay is not. Because each mode is
%    narrowband, applying G_n at the centre frequency is the same
%    approximation the sibling instruments make.
%
% 2. Direct broadband term. The injected jet flow minus the broadband part
%    of the mouth response, (1 - sum c_n) Qin, leaves through the window and
%    radiates as a monopole rho0 qdot / (4 pi r): the breath noise, the
%    attack transient, and everything the modal truncation dropped. Scaled
%    by direct_gain; this is the least principled term of the model, and
%    direct_gain = 0 turns it off.
%
% The instrument: window at the origin, bore along -y (foot down), listener
% at (listener_x_m, listener_y_m, listener_z_m).

rho0 = P.air_density;
c = P.speed_of_sound;
nT = size(Zr, 2);
n = basis.n;

lx = P.listener_x_m;
ly = P.listener_y_m;
lz = P.listener_z_m;

nap = length(basis.apX);
rap = zeros(nap, 1);
for a = 1:nap
  rap(a) = sqrt(lx^2 + (ly + basis.apX(a))^2 + lz^2);
end
listenerInside = min(rap) < 0.02;
rap = max(rap, 0.01);
rcen = rap(1);                        % the window

% modal analytic-signal weights (must mirror recorder_jet.m)
dScale = rho0 * (4 / pi) * sqrt(2 * P.flue_h_m * P.cutup_m) / basis.Sm;
cn = dScale * basis.B;
omd = basis.omega .* sqrt(max(1e-12, 1 - 1 ./ (4 * basis.Q.^2)));
lam = -basis.g / 2 + 1i * omd;
rho = (basis.g .* lam + basis.omega.^2) ./ (2i * omd);

H = zeros(1, n) + 0i;
for m = 1:n
  om = basis.omega(m);
  k = om / c;
  G = sum((1i * om * rho0 ./ (4 * pi * rap)) .* basis.apW(:, m) .* exp(-1i * k * (rap - rcen)));
  H(m) = G * (-2 * cn(m) * rho(m));
end

p = real(H) * Zr - imag(H) * Zi;      % 1 x nT

% direct broadband monopole from the window
dq = zeros(1, nT);
if nT > 2
  dq(2:nT - 1) = (Qin(3:nT) - Qin(1:nT - 2)) * (fsHi / 2);
end
p = p + P.direct_gain * (1 - sum(cn)) * rho0 / (4 * pi * rcen) * dq;

radinfo.gainAbs = abs(H);
radinfo.listenerInside = listenerInside;
radinfo.listenerDistM = rcen;
end