% Validate the bass model. There is no external reference implementation,
% so the checks are (a) closed-form identities the physics must satisfy,
% (b) structural identities of the coupled eigensolve, (c) exact
% reproducibility of the seeded walking line, and (d) 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 bass_validate.m

fails = 0;

% ---- 1. decoupled stiff string against the closed form ----
% With the bridge mobility driven to zero the body cannot move, and the
% coupled solve must hand back the textbook pinned-pinned stiff string,
% f_k = sqrt((k pi/Lv)^2 T/mu + (k pi/Lv)^4 EI/mu) / 2 pi, including the
% inharmonic stretch. Checked on the open A with an exaggerated B.
P = bass_defaults();
Q = P;
Q.bridge_mobility = 1e-9;
Q.inharmonicity = 2e-4;
basis = bass_basis(Q, 3, 0);
[~, ~, omS, ~, meta] = bass_assemble(Q, 3, 0);
errS = 0;
for k = 1:10
  fClosed = omS(k) / (2 * pi);
  fGot = bv_nearest(basis.f, fClosed);
  errS = max(errS, abs(fGot - fClosed) / fClosed);
end
stretch = (omS(10) / (10 * omS(1))) - 1;
ok = errS < 1e-6 && stretch > 0.005;
fails = fails + (~ok);
fprintf('stiff string: err %.1e vs closed form, partial-10 stretch %.3f %s\n', ...
  errS, stretch, bv_verdict(ok));

% ---- 2. eigensolve identities on the full coupled system ----
% The Cholesky reduction must hand back mass-orthonormal shapes that
% diagonalize the stiffness: Phi' M Phi = I and Phi' K Phi = diag(om^2).
% Checked on the open E, the largest system of the default chart.
[M, Kmat, ~, ~, meta] = bass_assemble(P, 4, 0);
Rc = chol(M);
Am = (Rc') \ (Kmat / Rc);
Am = (Am + Am') / 2;
[Uv, Sv, ~] = svd(Am);
om2 = flipud(diag(Sv));
Phi = Rc \ fliplr(Uv);
eOrth = max(max(abs(Phi' * M * Phi - eye(meta.n))));
KP = Phi' * Kmat * Phi;
eDiag = max(max(abs(KP - diag(om2)))) / max(om2);
ok = eOrth < 1e-8 && eDiag < 1e-10;
fails = fails + (~ok);
fprintf('eigensolve: orthonormality %.1e, diagonalization %.1e (n = %d) %s\n', ...
  eOrth, eDiag, meta.n, bv_verdict(ok));

% ---- 3. a midpoint pluck silences the even partials ----
% b_k ~ sin(k pi xp/Lv): plucking at the middle of the vibrating length
% puts every even string mode at a node of the initial shape, so the
% second partial must all but vanish from the rendered pressure, while an
% ordinary pluck position leaves it strong. This exercises the projection,
% the synthesis, and the radiation end to end.
Q = P;
Q.normalize = 0;
Q.tail_s = 3.0;
Q.pluck_pos = 0.5;
[y1, b1, d1] = bass_render(Q);
fs = Q.sample_rate;
f1 = d1.fundHz;
f2 = bv_nearest(b1.f, 2 * f1);
relMid = 20 * log10(bv_level(y1, fs, f2) / bv_level(y1, fs, f1) + 1e-15);
Q.pluck_pos = 0.30;
[y2, b2, d2] = bass_render(Q);
f2b = bv_nearest(b2.f, 2 * d2.fundHz);
relOrd = 20 * log10(bv_level(y2, fs, f2b) / bv_level(y2, fs, d2.fundHz) + 1e-15);
ok = relMid < -35 && relOrd > -25;
fails = fails + (~ok);
fprintf('midpoint pluck: partial 2 at %.1f dB (ordinary pluck %.1f dB) %s\n', ...
  relMid, relOrd, bv_verdict(ok));

% ---- 4. two-way coupling: veering at a string-body coincidence ----
% Tuning the strongest corpus mode onto the open G fundamental must split
% it into a mixed pair (this is the wolf-note mechanism), and the splitting
% must grow as the square root of the bridge mobility, which is what a
% two-way coupled pair does and a one-way filter does not.
d05 = bv_split(P, 0.5);
d20 = bv_split(P, 2.0);
ratio = d20.split / d05.split;
ok = abs(ratio - 2) < 0.25 && d20.mix > 0.15;
fails = fails + (~ok);
fprintf('veering: split %.2f -> %.2f Hz for 4x mobility (ratio %.2f, mix %.2f) %s\n', ...
  d05.split, d20.split, ratio, d20.mix, bv_verdict(ok));

% ---- 5. far-field falloff: -6 dB per doubling of distance ----
% RMS over the render, with the two listeners on the same ray from the
% bridge so the dipole directivity is held fixed. At the fundamental the
% dipole near-field term still contributes a few hundredths of a dB at
% 4 m, which is physics rather than error; the tolerance covers it.
Q = P;
Q.normalize = 0;
Q.tail_s = 2.5;
Q.listener_x_m = 1.2; Q.listener_y_m = 0.75; Q.listener_z_m = 3.0;
[ya, ~, ~] = bass_render(Q);
Q.listener_x_m = 2.4; Q.listener_y_m = 1.5; Q.listener_z_m = 6.0;
[yb, ~, ~] = bass_render(Q);
drop = 10 * log10(sum(yb.^2) / sum(ya.^2));
ok = abs(drop + 6.02) < 0.35;
fails = fails + (~ok);
fprintf('falloff 3.3 m -> 6.6 m: %.2f dB rms (target -6.02) %s\n', drop, bv_verdict(ok));

% ---- 6. the seeded walking line is exactly reproducible ----
Q = P;
Q.note_count = 8;
[t1, m1, s1, st1, fo1, g1] = bass_notes(Q);
[t2, m2, s2, st2, fo2, g2] = bass_notes(Q);
same = max(abs(t1 - t2)) == 0 && max(abs(m1 - m2)) == 0 ...
  && max(abs(s1 - s2)) == 0 && max(abs(fo1 - fo2)) == 0;
ok = same ...
  && g1.count == 8 ...
  && sum(m1) == 293 ...
  && abs(g1.meanForce - 14.014948614) < 1e-6 ...
  && abs(t1(2) - 0.623242473696) < 1e-9;
fails = fails + (~ok);
fprintf('walk (seed %d): %d notes, sum midi %d, mean force %.6f N %s\n', ...
  Q.seed, g1.count, sum(m1), g1.meanForce, bv_verdict(ok));

% ---- 7. pinned end-to-end diagnostics at the defaults ----
[y, basis, dinfo] = bass_render(P);
ok = abs(dinfo.peakPa - 0.270351) < 0.001 ...
  && abs(dinfo.fundHz - 54.4515) < 0.005 ...
  && basis.n == 61;
fails = fails + (~ok);
fprintf('pinned render: peak %.6f Pa (target %.6f), fund %.4f Hz (target %.4f), %d modes %s\n', ...
  dinfo.peakPa, 0.270351, dinfo.fundHz, 54.4515, basis.n, bv_verdict(ok));

% ---- 8. the note decays ----
% An open-string pizzicato rings for seconds, so the tail is not silent;
% but the energy in the last half second must sit well below the opening
% half second, which carries the attack.
fs = P.sample_rate;
n = length(y);
head = y(1:round(0.5 * fs));
tail = y(n - round(0.5 * fs):n);
dropE = 10 * log10(sum(tail.^2) / sum(head.^2));
ok = dropE < -15;
fails = fails + (~ok);
fprintf('note decay: last/first half second %.1f dB %s\n', dropE, bv_verdict(ok));

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

function fv = bv_nearest(fvec, target)
% The entry of fvec nearest target.
fv = fvec(1);
best = abs(fvec(1) - target);
for i = 2:length(fvec)
  d = abs(fvec(i) - target);
  if d < best
    best = d;
    fv = fvec(i);
  end
end
end

function a = bv_level(y, fs, f)
% Amplitude of y at frequency f over the first two seconds (plain
% quadrature dot products; the window is long against 1/f).
n = min(length(y), floor(2.0 * fs));
t = (0:n - 1) / fs;
seg = y(1:n);
a = sqrt((seg * (cos(2 * pi * f * t))')^2 + (seg * (sin(2 * pi * f * t))')^2);
end

function r = bv_split(P, mobility)
% Frequency splitting of the coupled pair when the strongest corpus mode
% is tuned onto the open G fundamental, plus how mixed the pair is.
Q = P;
Q.bridge_mobility = mobility;
Q.body_trim = 97.9989 / 97.0;
basis = bass_basis(Q, 1, 0);
f0 = 97.9989;
% the two coupled modes nearest the coincidence
d = abs(basis.f - f0);
[~, order] = sort(d);
i1 = order(1);
i2 = order(2);
r.split = abs(basis.f(i1) - basis.f(i2));
r.mix = min(min(basis.bodyFrac(i1), 1 - basis.bodyFrac(i1)), ...
  min(basis.bodyFrac(i2), 1 - basis.bodyFrac(i2)));
end

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