1function out = recorder_bore(P, geo, fing, fvec)
2% RECORDER_BORE Transmission-line acoustics of the bore at given frequencies.
3%
4% The bore is a chain of short lossy cylindrical segments (the linear taper
5% is discretized at 6 mm pitch) with the eight tone holes as lumped shunts:
6% an open hole is the air mass of its chimney plus two end corrections and a
7% baffled radiation resistance, a closed hole is the compliance of its
8% chimney volume, and a pinched thumb hole is an open hole at a reduced
9% effective radius. Viscothermal wall losses enter through the standard
10% wide-tube propagation constant,
11% Gamma = i k (1 + (1 - i) lossFac), lossFac = (dv / 2a)(1 + (gam-1)/sqrt(Pr)),
12% with dv = sqrt(2 nu / omega) the viscous boundary layer. The foot is an
13% unflanged open end; the window is an air mass (mouth end correction) with
14% a baffled monopole resistance.
15%
16% Everything is vectorized over the frequency vector fvec (row, Hz).
17%
18% Returns, per frequency:
19% out.Y series-loop admittance 1 / (Z_window + Z_bore): the passive
20% resonator the jet drives, whose peaks are the playable notes
21% out.Zbore bore input impedance at the window
22% out.W (nap x nf) complex flow out of each radiating aperture per
23% unit loop flow into the bore (aperture 1 is the window with
24% weight -1, then each open hole, then the foot)
25% out.apX (1 x nap) axial position of each radiating aperture, m
26% out.profX (1 x nsta) axial stations along the bore
27% out.profP (nsta x nf) complex pressure at each station per unit inflow
28% (the standing-wave profile of each resonance when evaluated
29% at its frequency)
31c = P.speed_of_sound;
32rho0 = P.air_density;
33nf = length(fvec);
34om = 2 * pi * reshape(fvec, 1, nf);
35k = om / c;
36nu = 1.5e-5; % kinematic viscosity of air
37gam = 1.4;
38sqPr = sqrt(0.71);
39dv = sqrt(2 * nu ./ max(om, 1)); % viscous boundary layer thickness
41L = geo.boreLen;
42rH = geo.rHead;
43rF = geo.rFoot;
45% ---- segment list: window (x = 0) to foot (x = L), split at the holes ----
46bounds = [0, geo.holePos, L];
47segX0 = [];
48segX1 = [];
49for j = 1:length(bounds) - 1
50 span = bounds(j + 1) - bounds(j);
51 if span <= 1e-9
52 continue;
53 end
54 m = max(1, ceil(span / 0.006));
55 e = bounds(j) + span * (0:m) / m;
56 segX0 = [segX0, e(1:m)];
57 segX1 = [segX1, e(2:m + 1)];
58end
59nseg = length(segX0);
61% ---- hole shunt impedances ----
62% holeZ(j, :) for holes 1..8; holeOpen(j) marks a radiating (open) hole.
63holeZ = zeros(8, nf) + 0i;
64holeOpen = zeros(1, 8);
65for j = 1:8
66 b = geo.holeRad(j);
67 t = geo.holeChim(j);
68 s = fing(j);
69 if s >= 0.25
70 pinched = s < 0.75;
71 if pinched
72 b = b * P.thumb_pinch; % pinched: reduced effective radius
73 end
74 Sh = pi * b^2;
75 teff = t + 1.5 * b; % inner + outer end corrections
76 holeZ(j, :) = 1i * om * rho0 * teff / Sh + rho0 * c * (k * b).^2 / (2 * Sh);
77 if pinched
78 % The thumb crack is a narrow slit, not a round hole: its viscous
79 % resistance is what damps the fundamental (whose pressure maximum
80 % the vent sits near) while barely touching the second register,
81 % which is the whole point of pinching.
82 holeZ(j, :) = holeZ(j, :) + P.pinch_loss * rho0 * c / Sh;
83 end
84 holeOpen(j) = 1;
85 else
86 V = pi * b^2 * t;
87 holeZ(j, :) = -1i * rho0 * c^2 ./ (max(om, 1) * V);
88 end
89end
91% ---- upward march: input impedance seen from the window ----
92aF = rF;
93SF = pi * aF^2;
94ZcF = rho0 * c / SF;
95Z = ZcF * ((k * aF).^2 / 4 + 1i * 0.61 * (k * aF)); % unflanged foot
97% precompute per-segment exp(Gamma l) and Zc, reused by both passes
98segE = zeros(nseg, nf) + 0i;
99segZc = zeros(nseg, 1);
100for e = 1:nseg
101 xm = (segX0(e) + segX1(e)) / 2;
102 a = rH + (rF - rH) * xm / L;
103 S = pi * a^2;
104 segZc(e) = rho0 * c / S;
105 lossFac = (dv / (2 * a)) * (1 + (gam - 1) / sqPr);
106 G = 1i * k .* (1 + lossFac) + k .* lossFac;
107 segE(e, :) = exp(G * (segX1(e) - segX0(e)));
108end
110% which hole (if any) sits at the top boundary of each segment
111holeAt = zeros(1, nseg);
112for j = 1:8
113 idx = find(abs(segX0 - geo.holePos(j)) < 1e-9);
114 if length(idx) == 1
115 holeAt(idx) = j;
116 end
117end
119for e = nseg:-1:1
120 Zc = segZc(e);
121 E = segE(e, :);
122 E2 = E .* E;
123 th = (E2 - 1) ./ (E2 + 1);
124 Z = Zc * (Z + Zc * th) ./ (Zc + Z .* th);
125 j = holeAt(e);
126 if j > 0
127 Zh = holeZ(j, :);
128 Z = Z .* Zh ./ (Z + Zh);
129 end
130end
131Zbore = Z;
133% ---- window impedance and the series-loop admittance ----
134Zw = 1i * om * rho0 * geo.deltaM / geo.Sm + rho0 * c * k.^2 / (2 * pi);
135Y = 1 ./ (Zw + Zbore);
137% ---- downward pass: pressure profile and aperture flows per unit inflow ----
138nsta = nseg + 1;
139profX = zeros(1, nsta);
140profP = zeros(nsta, nf) + 0i;
141p = Zbore; % pressure at the bore entry, U = 1
142U = ones(1, nf) + 0i;
143profX(1) = 0;
144profP(1, :) = p;
145apX = 0;
146W = -(ones(1, nf) + 0i); % window: flow out = -inflow
147for e = 1:nseg
148 % at the top boundary of this segment, shunt flow leaves through a hole
149 j = holeAt(e);
150 if j > 0
151 q = p ./ holeZ(j, :);
152 U = U - q;
153 if holeOpen(j)
154 apX = [apX; segX0(e)];
155 W = [W; q];
156 end
157 end
158 Zc = segZc(e);
159 E = segE(e, :);
160 ch = (E + 1 ./ E) / 2;
161 sh = (E - 1 ./ E) / 2;
162 p2 = p .* ch - Zc * U .* sh;
163 U = U .* ch - (p / Zc) .* sh;
164 p = p2;
165 profX(e + 1) = segX1(e);
166 profP(e + 1, :) = p;
167end
168apX = [apX; L];
169W = [W; U]; % what remains flows out the foot
171out.f = reshape(fvec, 1, nf);
172out.Y = Y;
173out.Zbore = Zbore;
174out.W = W;
175out.apX = reshape(apX, length(apX), 1);
176out.profX = profX;
177out.profP = profP;
178end