function [Zr, Zi, Qin, jinfo] = recorder_jet(P, basis, spanS, noteOrd, fsHi)
% RECORDER_JET Time-domain integration of the jet / bore feedback loop.
%
% This is the one stage of the model that cannot be closed-form: a recorder
% is a self-oscillator, so the loop is integrated sample by sample at fsHi.
% The loop, following the jet-drive model of Verge and Fabre:
%
%   1. Breath pressure p_f(t) (attack/release envelope, vibrato) sets the
%      jet exit velocity Uj = sqrt(2 p_f / rho0) (Bernoulli through the flue).
%   2. The acoustic flow U through the mouth gives a transverse velocity
%      v = U / S_m at the flue exit, which seeds a displacement wave on the
%      jet. The wave convects at jet_conv Uj (so the transit delay varies
%      with breath, which is what makes attacks speak the way they do),
%      grows by exp(jet_growth W / h), and is low-passed (two poles) at
%      jet_stcut Uj / W, the jet's own instability cutoff: frequencies the
%      jet cannot amplify must not drive the loop, and a one-pole rolloff
%      is shallow enough that an all-holes-closed fingering squeaks to its
%      fifth resonance instead of speaking the fundamental. Turbulence
%      enters here as seeded displacement noise scaling with Uj^1.5.
%   3. The displaced jet splits at the labium: the flow injected inside is
%      Qin = w_j b Uj (1 + tanh((eta - y0) / b)), b the Bickley half-width
%      grown over the cutup, y0 the labium offset.
%   4. The injection drives the bore as the jet-drive pressure source
%      dp = (rho0 delta_d / S_m) dQin/dt. Combined with the modal admittance
%      Y_n(s) = B_n s / (s^2 + g_n s + omega_n^2), the flow response of one
%      mode to the injected jet flow is
%        U_n(s) / Qin(s) = c_n s^2 / (s^2 + g_n s + omega_n^2),
%      c_n = B_n rho0 delta_d / S_m,
%      and this is discretized exactly as c_n (Qin - lowpass_n(Qin)) with
%      one impulse-invariant complex pole per mode: no discrete derivative
%      appears anywhere. (Approximating the two s factors by first
%      differences at fsHi is Nyquist-unstable by three orders of magnitude,
%      which is a property of the discretization, not of the physics.) The
%      broadband c_n Qin part is the mouth spitting the injected flow
%      straight back out through the window mass; the complex pole part is
%      the resonant flow the apertures radiate mode by mode.
%   5. The nonlinear vortex loss at the window, a series resistance
%      R = 0.5 rho0 K_v |v| / S_m, enters as amplitude-dependent extra
%      damping B_n R on each mode; with the tanh it caps the amplitude.
%      The total U closes the loop.
%
% Oscillation is not prescribed anywhere: it grows out of the noise when the
% loop gain crosses one, locks near a bore resonance where the phase closes
% (jet transit near a quarter period), and overblows to the second resonance
% when breath pressure carries the transit time out of the first's range.
%
% Returns the complex modal flow states Zr + i Zi (n x nT; the analytic
% signal 2 z_n of each modal flow, which recorder_radiate.m needs for its
% phase-correct aperture sum), the injected jet flow Qin (1 x nT, the direct
% broadband source), and diagnostics including the playing frequency
% estimated from zero crossings of the sustained half of the note.

rho0 = P.air_density;
dt = 1 / fsHi;
nT = max(16, floor(spanS * fsHi));
n = basis.n;

W = P.cutup_m;
h = P.flue_h_m;
Sm = basis.Sm;
bj = 0.41 * (h^2 * W)^(1 / 3);        % Bickley jet half-width at the labium
growth = exp(min(2.5, P.jet_growth * W / h));
dScale = rho0 * (4 / pi) * sqrt(2 * h * W) / Sm;
lossC = 0.5 * rho0 * P.vortex_loss / Sm;   % series resistance per unit |v|
y0 = P.labium_offset_m;
jg = P.jet_gain;

% ---- breath pressure envelope and derived vectors (all precomputed) ----
t = (0:nT - 1) * dt;
tauA = max(1e-4, P.attack_ms * 1e-3) / 2.3;
tauR = max(1e-4, P.release_ms * 1e-3) / 2.3;
env = 1 - exp(-t / tauA);
rel = find(t > P.note_dur_s);
env(rel) = env(rel) .* exp(-(t(rel) - P.note_dur_s) / tauR);
vib = 1 + P.vibrato_depth * sin(2 * pi * P.vibrato_hz * t);
% Breath tracking: players blow higher notes harder. The jet transit time
% must stay near a quarter period of the note, so the pressure that centers
% the lock range grows as f^2; breath_track < 2 undershoots that slightly,
% which is what measured blowing pressures do.
trackFac = (basis.fTarget / 523.25)^P.breath_track;
pf = P.breath_pa * trackFac * env .* vib;
Uj = sqrt(2 * max(pf, 0) / rho0);
UjE = max(Uj, 1.0);                   % guarded: receptivity and cutoff need Uj > 0
tauS = min(0.005, W ./ (P.jet_conv * UjE)) * fsHi;
tauS = max(2, tauS);
aJet = 1 - exp(-2 * pi * (P.jet_stcut * UjE / W) * dt);   % jet instability cutoff
etaAmp = P.noise_level * h * (UjE / 8).^1.5;

% ---- seeded turbulence noise, filtered in the loop ----
state = mod(floor(abs(P.seed)) + 7919 * noteOrd, 4294967296);
[u, state] = recorder_rng(state, nT);
white = 2 * u.' - 1;
aN1 = 1 - exp(-2 * pi * 4000 * dt);   % turbulence band: 400 Hz .. 4 kHz
aN2 = 1 - exp(-2 * pi * 400 * dt);

% ---- per-mode discrete resonator coefficients (impulse invariant) ----
% U_n = c_n (Qin - 2 Re(rho_n z_n)), z_n' = lam_n z_n + Qin, with rho_n the
% residue of (g s + omega^2) / (s^2 + g s + omega^2) at lam_n.
omd = basis.omega .* sqrt(max(1e-12, 1 - 1 ./ (4 * basis.Q.^2)));
lam = -basis.g / 2 + 1i * omd;
Az = exp(lam * dt);
ar = real(Az);
ai = imag(Az);
cn = dScale * basis.B;
rho = (basis.g .* lam + basis.omega.^2) ./ (2i * omd);
% modal flow per state: U_res,n = wr_n zr_n + wi_n zi_n (the -2 c Re(rho z))
wr = -2 * cn .* real(rho);
wi = 2 * cn .* imag(rho);
cSum = sum(cn);

% ---- the loop ----
Zr = zeros(n, nT);
Zi = zeros(n, nT);
Qin = zeros(1, nT);
vh = zeros(1, nT);
zr = zeros(n, 1);
zi = zeros(n, 1);
vprev = 0;
vlp = 0;
vlp2 = 0;
nlp1 = 0;
nlp2 = 0;
for k = 1:nT
  % turbulence: white noise bandpassed 400 Hz .. 4 kHz (two one-poles)
  nlp1 = nlp1 + aN1 * (white(k) - nlp1);
  nlp2 = nlp2 + aN2 * (nlp1 - nlp2);
  % transverse acoustic velocity at the flue exit, delayed by the jet
  % transit and low-passed at the jet's instability cutoff
  dk = k - tauS(k);
  vdel = 0;
  if dk >= 1
    i0 = floor(dk);
    wfr = dk - i0;
    vdel = vh(i0) * (1 - wfr) + vh(min(i0 + 1, nT)) * wfr;
  end
  vlp = vlp + aJet(k) * (vdel - vlp);
  vlp2 = vlp2 + aJet(k) * (vlp - vlp2);
  eta = jg * (h / UjE(k)) * growth * vlp2 + etaAmp(k) * (nlp1 - nlp2);
  qin = P.flue_w_m * bj * Uj(k) * (1 + tanh((eta - y0) / bj));
  damp = 1 - (dt / 2) * (lossC * abs(vprev)) * basis.B;
  zr2 = damp .* (ar .* zr - ai .* zi) + dt * qin;
  zi = damp .* (ar .* zi + ai .* zr);
  zr = zr2;
  Zr(:, k) = zr;
  Zi(:, k) = zi;
  U = cSum * qin + wr.' * zr + wi.' * zi;
  vprev = U / Sm;
  vh(k) = vprev;
  Qin(k) = qin;
end

% ---- diagnostics: playing frequency from the sustained half ----
k0 = max(1, floor(0.35 * P.note_dur_s * fsHi));
k1 = min(nT, floor(0.95 * P.note_dur_s * fsHi));
seg = vh(k0:k1);
mu = sum(seg) / length(seg);
seg = seg - mu;
nzc = 0;
kFirst = 0;
kLast = 0;
for k = 2:length(seg)
  if seg(k - 1) <= 0 && seg(k) > 0
    nzc = nzc + 1;
    if kFirst == 0
      kFirst = k;
    end
    kLast = k;
  end
end
if nzc > 2
  jinfo.fPlay = (nzc - 1) * fsHi / (kLast - kFirst);
else
  jinfo.fPlay = 0;
end
jinfo.fTarget = basis.fTarget;
if jinfo.fPlay > 0
  jinfo.cents = 1200 * log(jinfo.fPlay / basis.fTarget) / log(2);
else
  jinfo.cents = 0;
end
jinfo.vacRms = sqrt(sum(seg.^2) / length(seg));
jinfo.ujMax = max(Uj);
jinfo.qinMean = sum(Qin) / nT;
end