function basis = recorder_basis(P, note)
% RECORDER_BASIS Modal fit of the passive resonator for one fingering.
%
% The series-loop admittance Y(omega) = 1 / (Z_window + Z_bore) from
% recorder_bore.m relates the acoustic flow U through the mouth to the
% pressure the jet imposes across it; its peaks are the bore resonances the
% jet can lock to. Each peak is fitted with one second-order section,
%   Y_n(omega) = i omega B_n / (omega_n^2 - omega^2 + i omega g_n),
% so that in the time domain U_n'' + g_n U_n' + omega_n^2 U_n = B_n dp/dt.
% omega_n comes from parabolic interpolation of log|Y| on the grid, g_n from
% the half-power bandwidth, and B_n = g_n Re(Y(omega_n)) so the fitted peak
% height matches. The fit discards nothing a recorder player can reach:
% every local maximum below max_mode_hz is kept, up to max_modes, however
% weak (a cross-fingered lattice makes the upper resonances weak but they
% still color the attack).
%
% For radiation the basis stores, per mode, the complex flow out of each
% radiating aperture (the window, each open hole, the foot) per unit loop
% flow, evaluated at the mode's centre frequency, plus the pressure
% standing-wave profile along the bore for the mode-shape plot.

geo = recorder_geometry(P);
note = min(max(1, round(note)), geo.nNotes);
fing = geo.fingerings(note, :);

fgrid = 60:P.grid_hz:P.max_mode_hz;
out = recorder_bore(P, geo, fing, fgrid);
mag = abs(out.Y);

% ---- peak picking with parabolic refinement ----
fpk = [];
qpk = [];
ypk = [];
nf = length(fgrid);
for i = 2:nf - 1
  if ~(mag(i) > mag(i - 1) && mag(i) >= mag(i + 1))
    continue;
  end
  y1 = log(mag(i - 1)); y2 = log(mag(i)); y3 = log(mag(i + 1));
  den = y1 - 2 * y2 + y3;
  sh = 0;
  if abs(den) > 1e-12
    sh = min(0.5, max(-0.5, 0.5 * (y1 - y3) / den));
  end
  fc = fgrid(i) + sh * P.grid_hz;
  yc = exp(y2 - 0.25 * (y1 - y3) * sh);
  % local bandwidth from the two neighbouring grid points: for one section,
  % |Y(f_n + df)|^2 / |Y(f_n)|^2 = g^2 / ((4 pi df)^2 + g^2), so each
  % neighbour gives g = 4 pi |df| sqrt(r / (1 - r)). Staying local keeps the
  % estimate off the skirts of adjacent resonances, which a half-power scan
  % walks onto.
  gs = [];
  for s = [-1, 1]
    df = fgrid(i + s) - fc;
    r = (mag(i + s) / yc)^2;
    if r < 0.999
      gs = [gs, 4 * pi * abs(df) * sqrt(r / (1 - r))];
    end
  end
  if length(gs) == 0
    g = 2 * pi * P.grid_hz;
  else
    g = sum(gs) / length(gs);
  end
  fpk = [fpk, fc];
  qpk = [qpk, 2 * pi * fc / g];
  ypk = [ypk, yc];
end

n = min(length(fpk), round(P.max_modes));
if n == 0
  error('No bore resonances below max_mode_hz -- check the geometry.');
end
fpk = fpk(1:n);
qpk = qpk(1:n);
ypk = ypk(1:n);

% ---- aperture weights and standing-wave profiles at the mode centres ----
mo = recorder_bore(P, geo, fing, fpk);

basis.n = n;
basis.note = note;
basis.f = reshape(fpk, n, 1);
basis.omega = 2 * pi * basis.f;
basis.Q = reshape(qpk, n, 1);
basis.g = basis.omega ./ basis.Q;
basis.B = basis.g .* reshape(ypk, n, 1);
basis.apX = mo.apX;                   % nap x 1 axial positions
basis.apW = mo.W;                     % nap x n complex flow weights
basis.profX = mo.profX;               % 1 x nsta
basis.profP = abs(mo.profP);          % nsta x n standing-wave magnitudes
basis.fTarget = geo.fTarget(note);
basis.Sm = geo.Sm;
basis.gridF = fgrid;
basis.gridY = mag;
end