% Validate the recorder model. There is no external reference implementation,
% so the checks are (a) closed-form identities the passive acoustics must
% satisfy, (b) internal identities of the modal fit, (c) emergent behavior
% the self-oscillating jet must reproduce (register selection, overblowing),
% (d) exact reproducibility of the seeded turbulence, and (e) pinned
% end-to-end diagnostics from the default configuration, recorded when the
% model was built and treated as a regression reference since.
% Run with: numbl run recorder_validate.m

fails = 0;
P = recorder_defaults();

% ---- 1. uniform open-open cylinder against the closed form ----
% With the tone holes shrunk to nothing and the taper removed, the bore
% admittance peaks must sit at n c / 2 (L + 0.61 a): the open-open pipe with
% the unflanged end correction at the foot (the window end is where Zbore is
% measured, so it carries no correction of its own). Viscothermal losses
% shift these by well under a percent.
Q = P;
Q.bore_r_head_m = 0.006;
Q.bore_r_foot_m = 0.006;
geo = recorder_geometry(Q);
geo.holeRad = 1e-9 * ones(1, 8);
f = 100:1:2500;
out = recorder_bore(Q, geo, zeros(1, 8), f);
mag = abs(1 ./ out.Zbore);
pk = [];
for i = 2:length(f) - 1
  if mag(i) > mag(i - 1) && mag(i) >= mag(i + 1)
    pk = [pk, f(i)];
  end
end
c = Q.speed_of_sound;
Leff = geo.boreLen + 0.61 * 0.006;
fClosed = (1:3) * c / (2 * Leff);
% viscothermal dispersion: the phase velocity is c / (1 + lossFac), with
% lossFac = (dv / 2a)(1 + (gam - 1) / sqrt(Pr)) evaluated at each resonance
dv = sqrt(2 * 1.5e-5 ./ (2 * pi * fClosed));
fClosed = fClosed ./ (1 + (dv / 0.012) * (1 + 0.4 / sqrt(0.71)));
errC = max(abs(pk(1:3) - fClosed) ./ fClosed);
ok = errC < 3e-3;
fails = fails + (~ok);
fprintf('open-open cylinder: %.1f %.1f %.1f Hz vs closed form %.1f %.1f %.1f (err %.1e) %s\n', ...
  pk(1), pk(2), pk(3), fClosed(1), fClosed(2), fClosed(3), errC, rv_verdict(ok));

% ---- 2. the modal fit reconstructs the admittance it was fitted to ----
% At each fitted centre frequency the modal sum must reproduce the computed
% series-loop admittance. The check covers the first four resonances, which
% carry the oscillation; above them the per-peak fit overcounts by the
% skirts of neighbouring modes (up to ~20% on the weakest peaks), which is
% inherent to independent-peak fitting and acoustically negligible there.
basis = recorder_basis(P, 5);
geo = recorder_geometry(P);
bo = recorder_bore(P, geo, geo.fingerings(5, :), reshape(basis.f, 1, basis.n));
eFit = 0;
for m = 1:min(4, basis.n)
  om = basis.omega(m);
  Yfit = 0;
  for j = 1:basis.n
    Yfit = Yfit + 1i * om * basis.B(j) / (basis.omega(j)^2 - om^2 + 1i * om * basis.g(j));
  end
  eFit = max(eFit, abs(abs(Yfit) - abs(bo.Y(m))) / abs(bo.Y(m)));
end
ok = eFit < 0.06;
fails = fails + (~ok);
fprintf('modal fit: worst |Y| error over the first four peaks %.1e %s\n', eFit, rv_verdict(ok));

% ---- 3. register selection is emergent ----
% Three fingerings that must each speak their own register at the default
% (tracked) breath: C5 all closed (the easiest fingering to squeak, so this
% is the sharpest test of the jet's instability cutoff), G5, and pinched E6,
% whose fundamental must be suppressed by the vent and not played.
Q = P;
Q.note_dur_s = 0.6;
Q.tail_s = 0.2;
fsHi = P.sample_rate * P.oversample;
centsAbs = zeros(1, 3);
notes3 = [1 5 11];
for i = 1:3
  b = recorder_basis(Q, notes3(i));
  [~, ~, ~, j] = recorder_jet(Q, b, 0.7, 1, fsHi);
  centsAbs(i) = abs(j.cents);
end
ok = max(centsAbs) < 60;
fails = fails + (~ok);
fprintf('registers: C5 %+.0f c, G5 %+.0f c, E6 %+.0f c of equal temperament %s\n', ...
  centsAbs(1), centsAbs(2), centsAbs(3), rv_verdict(ok));

% ---- 4. overblowing is emergent ----
% Quadrupling the breath pressure on the G5 fingering must carry the jet
% transit time out of the fundamental's lock range and onto the second
% resonance: an octave (plus the resonator's inharmonicity) up. Nothing in
% the model prescribes this.
Q = P;
Q.note_dur_s = 0.6;
Q.tail_s = 0.2;
Q.breath_pa = 4 * P.breath_pa;
b = recorder_basis(Q, 5);
[~, ~, ~, j] = recorder_jet(Q, b, 0.7, 1, fsHi);
ratio = j.fPlay / b.f(1);
ok = ratio > 1.9 && ratio < 2.1;
fails = fails + (~ok);
fprintf('overblow at 4x breath: %.1f Hz = %.3f x f1 %s\n', j.fPlay, ratio, rv_verdict(ok));

% ---- 5. the seeded turbulence is exactly reproducible ----
Q = P;
Q.note_dur_s = 0.4;
Q.tail_s = 0.15;
[y1, ~, ~] = recorder_render(Q);
[y2, ~, ~] = recorder_render(Q);
Q.seed = 2;
[y3, ~, ~] = recorder_render(Q);
same = max(abs(y1 - y2)) == 0;
diff = max(abs(y1 - y3)) > 0;
ok = same && diff;
fails = fails + (~ok);
fprintf('seed: identical rerun %d, different seed differs %d %s\n', same, diff, rv_verdict(ok));

% ---- 6. far-field falloff: -6 dB per doubling of distance ----
% RMS over the render, both listeners on one ray from the window so every
% aperture monopole keeps its interference geometry.
Q = P;
Q.normalize = 0;
Q.note_dur_s = 0.5;
Q.tail_s = 0.15;
Q.listener_x_m = 0.5; Q.listener_y_m = 0.4; Q.listener_z_m = 1.0;
[y1, ~, ~] = recorder_render(Q);
Q.listener_x_m = 1.0; Q.listener_y_m = 0.8; Q.listener_z_m = 2.0;
[y2, ~, ~] = recorder_render(Q);
drop = 10 * log10(sum(y2.^2) / sum(y1.^2));
ok = abs(drop + 6.02) < 0.35;
fails = fails + (~ok);
fprintf('falloff 1.2 m -> 2.4 m: %.2f dB rms (target -6.02) %s\n', drop, rv_verdict(ok));

% ---- 7. pinned end-to-end diagnostics at the defaults ----
[y, basis, dinfo] = recorder_render(P);
ok = abs(dinfo.fPlay(1) - 781.88) < 1.5 ...
  && abs(dinfo.peakPa - 1.7982) < 0.05 ...
  && basis.n == 12 ...
  && abs(basis.f(1) - 783.94) < 0.1;
fails = fails + (~ok);
fprintf('pinned render: played %.2f Hz (target %.2f), peak %.4f Pa (target %.4f), %d modes %s\n', ...
  dinfo.fPlay(1), 781.88, dinfo.peakPa, 1.7982, basis.n, rv_verdict(ok));

% ---- 8. the note stops when the breath stops ----
fs = P.sample_rate;
n = length(y);
sus = y(round(0.5 * fs):round(1.0 * fs));
tail = y(n - round(0.2 * fs):n);
dropE = 10 * log10(sum(tail.^2) / length(tail) / (sum(sus.^2) / length(sus)));
ok = dropE < -40;
fails = fails + (~ok);
fprintf('release: last 200 ms %.1f dB re sustain %s\n', dropE, rv_verdict(ok));

if fails > 0
  fprintf('%d CHECK(S) FAILED\n', fails);
else
  fprintf('all checks passed\n');
end

function s = rv_verdict(ok)
if ok
  s = 'ok';
else
  s = 'MISMATCH';
end
end