% Validate the triangle model. There is no external reference implementation,
% so the checks are (a) closed-form identities the physics must satisfy,
% (b) structural identities of the FEM eigensolve, (c) exact reproducibility
% of the seeded roll process, 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 triangle_validate.m

fails = 0;

% ---- 1. straight free-free rod against closed forms ----
% The same element routines that build the triangle, assembled on a straight
% polyline, must reproduce the textbook free-free bar: bending frequencies
% (beta L)^2 kappa cL / (2 pi L^2), the first axial mode cL / 2L, the first
% torsional mode cT / 2L, and exactly three rigid-body modes per family.
E = 2e11; rho = 7850; nu = 0.30; d = 0.008; L = 0.44;
nodes = [linspace(0, L, 61); zeros(1, 61)];
cL = sqrt(E / rho);
kap = d / 4;
blTab = [4.73004074 7.85320462 10.99560784];
fbClosed = blTab.^2 * kap * cL / (2 * pi * L^2);
faClosed = cL / (2 * L);
ftClosed = sqrt(E / (2 * (1 + nu)) / rho) / (2 * L);
f0 = tv_eigfreqs(nodes, 0, E, rho, nu, d);
f1 = tv_eigfreqs(nodes, 1, E, rho, nu, d);
% the rigid-body zeros come back as numerical noise well below the first
% elastic mode (185 Hz here); 5 Hz separates them cleanly
nRigid0 = sum(f0 < 5);
nRigid1 = sum(f1 < 5);
e0 = f0(f0 > 1);
e1 = f1(f1 > 1);
errB = max(abs(e0(1:3).' - fbClosed) ./ fbClosed);
% rotary inertia lowers the FEM bending values slightly below the closed
% form, which neglects it; 0.5% covers that plus discretization
errA = abs(min(abs(e0 - faClosed)) / faClosed);
errT = abs(min(abs(e1 - ftClosed)) / ftClosed);
ok = nRigid0 == 3 && nRigid1 == 3 && errB < 5e-3 && errA < 2e-3 && errT < 2e-3;
fails = fails + (~ok);
fprintf('straight rod: bend err %.1e, axial %.1e, torsion %.1e, rigid %d/%d %s\n', ...
  errB, errA, errT, nRigid0, nRigid1, tv_verdict(ok));

% ---- 2. eigensolve identities on the triangle itself ----
% The Cholesky reduction must hand back mass-orthonormal shapes that
% diagonalize the stiffness: phi' M phi = I and phi' K phi = diag(om^2).
P = triangle_defaults();
geo = triangle_geometry(P);
[K, M] = triangle_fem(geo.nodes, 0, P.youngs_modulus, P.density, P.poisson_ratio, P.rod_diameter_m);
Rc = chol(M);
Am = (Rc') \ (K / Rc);
Am = (Am + Am') / 2;
[Uv, Sv, ~] = svd(Am);
om2 = flipud(diag(Sv));
Vv = fliplr(Uv);
Phi = Rc \ Vv;
eOrth = max(max(abs(Phi' * M * Phi - eye(length(om2)))));
KP = Phi' * K * 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 %s\n', eOrth, eDiag, tv_verdict(ok));

% ---- 3. Hertz contact-time scaling tau ~ v^(-1/5) ----
% With a beater much lighter than the engaged rod mass the reduced mass is a
% constant and the free-sphere scaling must emerge.
Q = P;
Q.beater_mass_kg = 5e-4;
[~, ~, ~, tau1, ~, ~, ~] = triangle_strikes(Q);
Q.strike_velocity = 2 * P.strike_velocity;
[~, ~, ~, tau2, ~, ~, ~] = triangle_strikes(Q);
ratio = tau2(1) / tau1(1);
ok = abs(ratio - 2^(-1 / 5)) < 0.03;
fails = fails + (~ok);
fprintf('hertz scaling: tau ratio %.4f (target %.4f) %s\n', ratio, 2^(-1 / 5), tv_verdict(ok));

% ---- 4. seeded roll is exactly reproducible ----
Q = P;
Q.strike_count = 12;
[t1r, v1r, j1r, ~, ~, s1r, g1] = triangle_strikes(Q);
[t2r, v2r, j2r, ~, ~, s2r, g2] = triangle_strikes(Q);
same = max(abs(t1r - t2r)) == 0 && max(abs(v1r - v2r)) == 0 && max(abs(s1r - s2r)) == 0;
ok = same ...
  && g1.count == 12 ...
  && abs(g1.meanV - 1.501211020) < 1e-6 ...
  && abs(sum(j1r) - 0.383069278968) < 1e-9 ...
  && abs(t1r(2) - 0.087732781124) < 1e-9;
fails = fails + (~ok);
fprintf('roll (seed %d): %d strikes, mean v %.6f m/s, sum J %.9f N s %s\n', ...
  Q.seed, g1.count, g1.meanV, sum(j1r), tv_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
% origin so every element dipole keeps its directivity factor: off the ray
% the in-plane dipoles fall faster than 1/r, which is physics rather than
% error. Modes are capped at 8 kHz so 2 m is beyond the Rayleigh distance
% of the frame; at 20 kHz it is not, and Fresnel-zone interference shifts
% the ratio by a few tenths of a dB, again physics rather than error.
Q = P;
Q.normalize = 0;
Q.tail_s = 1.0;
Q.max_mode_hz = 8000;
Q.listener_x_m = 0.6; Q.listener_y_m = 0.5; Q.listener_z_m = 1.8;
[y1, ~, ~] = triangle_render(Q);
Q.listener_x_m = 1.2; Q.listener_y_m = 1.0; Q.listener_z_m = 3.6;
[y2, ~, ~] = triangle_render(Q);
drop = 10 * log10(sum(y2.^2) / sum(y1.^2));
ok = abs(drop + 6.02) < 0.2;
fails = fails + (~ok);
fprintf('falloff 2 m -> 4 m: %.2f dB rms (target -6.02) %s\n', drop, tv_verdict(ok));

% ---- 6. pinned end-to-end diagnostics at the defaults ----
[y, basis, dinfo] = triangle_render(P);
ok = abs(dinfo.peakPa - 7.217688) < 0.02 ...
  && abs(basis.f(1) - 163.9200) < 0.1 ...
  && basis.n == 40;
fails = fails + (~ok);
fprintf('pinned render: peak %.4f Pa (target %.4f), mode1 %.2f Hz (target %.2f), %d modes %s\n', ...
  dinfo.peakPa, 7.217688, basis.f(1), 163.9200, basis.n, tv_verdict(ok));

% ---- 7. the ring decays ----
% A triangle rings for seconds, so the tail is not silent; but the energy in
% the last half second must sit well below the first half second's.
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 < -25;
fails = fails + (~ok);
fprintf('ring decay: last/first half second %.1f dB %s\n', dropE, tv_verdict(ok));

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

function f = tv_eigfreqs(nodes, family, E, rho, nu, d)
[K, M] = triangle_fem(nodes, family, E, rho, nu, d);
Rc = chol(M);
Am = (Rc') \ (K / Rc);
Am = (Am + Am') / 2;
[~, Sv, ~] = svd(Am);
f = sort(sqrt(diag(Sv))) / (2 * pi);
end

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