function out = recorder_bore(P, geo, fing, fvec)
% RECORDER_BORE Transmission-line acoustics of the bore at given frequencies.
%
% The bore is a chain of short lossy cylindrical segments (the linear taper
% is discretized at 6 mm pitch) with the eight tone holes as lumped shunts:
% an open hole is the air mass of its chimney plus two end corrections and a
% baffled radiation resistance, a closed hole is the compliance of its
% chimney volume, and a pinched thumb hole is an open hole at a reduced
% effective radius. Viscothermal wall losses enter through the standard
% wide-tube propagation constant,
%   Gamma = i k (1 + (1 - i) lossFac),  lossFac = (dv / 2a)(1 + (gam-1)/sqrt(Pr)),
% with dv = sqrt(2 nu / omega) the viscous boundary layer. The foot is an
% unflanged open end; the window is an air mass (mouth end correction) with
% a baffled monopole resistance.
%
% Everything is vectorized over the frequency vector fvec (row, Hz).
%
% Returns, per frequency:
%   out.Y      series-loop admittance 1 / (Z_window + Z_bore): the passive
%              resonator the jet drives, whose peaks are the playable notes
%   out.Zbore  bore input impedance at the window
%   out.W      (nap x nf) complex flow out of each radiating aperture per
%              unit loop flow into the bore (aperture 1 is the window with
%              weight -1, then each open hole, then the foot)
%   out.apX    (1 x nap) axial position of each radiating aperture, m
%   out.profX  (1 x nsta) axial stations along the bore
%   out.profP  (nsta x nf) complex pressure at each station per unit inflow
%              (the standing-wave profile of each resonance when evaluated
%              at its frequency)

c = P.speed_of_sound;
rho0 = P.air_density;
nf = length(fvec);
om = 2 * pi * reshape(fvec, 1, nf);
k = om / c;
nu = 1.5e-5;                          % kinematic viscosity of air
gam = 1.4;
sqPr = sqrt(0.71);
dv = sqrt(2 * nu ./ max(om, 1));      % viscous boundary layer thickness

L = geo.boreLen;
rH = geo.rHead;
rF = geo.rFoot;

% ---- segment list: window (x = 0) to foot (x = L), split at the holes ----
bounds = [0, geo.holePos, L];
segX0 = [];
segX1 = [];
for j = 1:length(bounds) - 1
  span = bounds(j + 1) - bounds(j);
  if span <= 1e-9
    continue;
  end
  m = max(1, ceil(span / 0.006));
  e = bounds(j) + span * (0:m) / m;
  segX0 = [segX0, e(1:m)];
  segX1 = [segX1, e(2:m + 1)];
end
nseg = length(segX0);

% ---- hole shunt impedances ----
% holeZ(j, :) for holes 1..8; holeOpen(j) marks a radiating (open) hole.
holeZ = zeros(8, nf) + 0i;
holeOpen = zeros(1, 8);
for j = 1:8
  b = geo.holeRad(j);
  t = geo.holeChim(j);
  s = fing(j);
  if s >= 0.25
    pinched = s < 0.75;
    if pinched
      b = b * P.thumb_pinch;          % pinched: reduced effective radius
    end
    Sh = pi * b^2;
    teff = t + 1.5 * b;               % inner + outer end corrections
    holeZ(j, :) = 1i * om * rho0 * teff / Sh + rho0 * c * (k * b).^2 / (2 * Sh);
    if pinched
      % The thumb crack is a narrow slit, not a round hole: its viscous
      % resistance is what damps the fundamental (whose pressure maximum
      % the vent sits near) while barely touching the second register,
      % which is the whole point of pinching.
      holeZ(j, :) = holeZ(j, :) + P.pinch_loss * rho0 * c / Sh;
    end
    holeOpen(j) = 1;
  else
    V = pi * b^2 * t;
    holeZ(j, :) = -1i * rho0 * c^2 ./ (max(om, 1) * V);
  end
end

% ---- upward march: input impedance seen from the window ----
aF = rF;
SF = pi * aF^2;
ZcF = rho0 * c / SF;
Z = ZcF * ((k * aF).^2 / 4 + 1i * 0.61 * (k * aF));   % unflanged foot

% precompute per-segment exp(Gamma l) and Zc, reused by both passes
segE = zeros(nseg, nf) + 0i;
segZc = zeros(nseg, 1);
for e = 1:nseg
  xm = (segX0(e) + segX1(e)) / 2;
  a = rH + (rF - rH) * xm / L;
  S = pi * a^2;
  segZc(e) = rho0 * c / S;
  lossFac = (dv / (2 * a)) * (1 + (gam - 1) / sqPr);
  G = 1i * k .* (1 + lossFac) + k .* lossFac;
  segE(e, :) = exp(G * (segX1(e) - segX0(e)));
end

% which hole (if any) sits at the top boundary of each segment
holeAt = zeros(1, nseg);
for j = 1:8
  idx = find(abs(segX0 - geo.holePos(j)) < 1e-9);
  if length(idx) == 1
    holeAt(idx) = j;
  end
end

for e = nseg:-1:1
  Zc = segZc(e);
  E = segE(e, :);
  E2 = E .* E;
  th = (E2 - 1) ./ (E2 + 1);
  Z = Zc * (Z + Zc * th) ./ (Zc + Z .* th);
  j = holeAt(e);
  if j > 0
    Zh = holeZ(j, :);
    Z = Z .* Zh ./ (Z + Zh);
  end
end
Zbore = Z;

% ---- window impedance and the series-loop admittance ----
Zw = 1i * om * rho0 * geo.deltaM / geo.Sm + rho0 * c * k.^2 / (2 * pi);
Y = 1 ./ (Zw + Zbore);

% ---- downward pass: pressure profile and aperture flows per unit inflow ----
nsta = nseg + 1;
profX = zeros(1, nsta);
profP = zeros(nsta, nf) + 0i;
p = Zbore;                            % pressure at the bore entry, U = 1
U = ones(1, nf) + 0i;
profX(1) = 0;
profP(1, :) = p;
apX = 0;
W = -(ones(1, nf) + 0i);              % window: flow out = -inflow
for e = 1:nseg
  % at the top boundary of this segment, shunt flow leaves through a hole
  j = holeAt(e);
  if j > 0
    q = p ./ holeZ(j, :);
    U = U - q;
    if holeOpen(j)
      apX = [apX; segX0(e)];
      W = [W; q];
    end
  end
  Zc = segZc(e);
  E = segE(e, :);
  ch = (E + 1 ./ E) / 2;
  sh = (E - 1 ./ E) / 2;
  p2 = p .* ch - Zc * U .* sh;
  U = U .* ch - (p / Zc) .* sh;
  p = p2;
  profX(e + 1) = segX1(e);
  profP(e + 1, :) = p;
end
apX = [apX; L];
W = [W; U];                           % what remains flows out the foot

out.f = reshape(fvec, 1, nf);
out.Y = Y;
out.Zbore = Zbore;
out.W = W;
out.apX = reshape(apX, length(apX), 1);
out.profX = profX;
out.profP = profP;
end