% Validate the cymbal2 model. There is no external reference
% implementation, so the checks are (a) closed-form limits the
% discretization must reproduce (the free circular plate, rigid-body
% motion on the washer, cantilever statics of the stick), (b) agreement
% with the linear-FEM eigenfrequencies published in Samejima (2021) for
% the modes that model resolves the same way, (c) structural identities
% of the transform and of the assembled update, (d) exact reproducibility,
% and (e) pinned end-to-end diagnostics from a reduced configuration,
% recorded when the model was built and treated as a regression reference.
% Run with: numbl run cymbal2_validate.m

fails = 0;
P = cymbal2_defaults();

% ---- 1. flat limit: the completely free circular plate ----
% With the curvature removed and the washer stiffness zeroed, the shell
% must reproduce the free plate's frequency equation (the 2 x 2 J_n / I_n
% determinant under zero radial moment and zero Kelvin-Kirchhoff shear).
% m = 0 is excluded: its limit as the free hole closes is logarithmic, so
% a finite hole shifts it far more than the m >= 1 modes.
Q = P;
Q.rise_m = 1e-9;
Q.washer_stiffness = 0;
Q.washer_radius_m = 0.010;
Q.n_radial = 96;
Dp = Q.youngs_modulus * Q.thickness_m^3 / (12 * (1 - Q.poisson_ratio^2));
scal = sqrt(Dp / (Q.density * Q.thickness_m)) / Q.radius_m^2 / (2 * pi);
errP = 0;
for m = 1:4
  f = cymbal2_modes(Q, m);
  fe = f(f > 0.5);
  lam = c2v_plate_roots(m, Q.poisson_ratio);
  fc = lam.^2 * scal;
  nc = min(3, min(length(fe), length(fc)));
  errP = max(errP, max(abs(fe(1:nc) - fc(1:nc)') ./ fc(1:nc)'));
end
ok = errP < 0.04;
fails = fails + (~ok);
fprintf('flat plate : max freq err %.1e over m = 1..4 %s\n', errP, c2v_verdict(ok));

% ---- 2. rigid limit: bounce and tilt on the washer ----
% With the plate stiffened a hundredfold the lowest m = 0 and m = 1 modes
% must approach the closed-form rigid-body values on the line spring:
% bounce f^2 = KW 2 pi c / (rho h pi (a^2 - c^2)) / (2 pi)^2 and tilt
% f^2 = 4 KW c^3 / (rho h (a^4 - c^4)) / (2 pi)^2.
Q = P;
Q.rise_m = 1e-9;
Q.youngs_modulus = P.youngs_modulus * 100;
Q.n_radial = 96;
a = Q.radius_m; c = Q.washer_radius_m;
rhoh = Q.density * Q.thickness_m;
KW = Q.washer_stiffness;
f0r = sqrt(KW * 2 * pi * c / (rhoh * pi * (a^2 - c^2))) / (2 * pi);
f1r = sqrt(4 * KW * c^3 / (rhoh * (a^4 - c^4))) / (2 * pi);
fa = cymbal2_modes(Q, 0); fa = fa(fa > 0.01);
fb = cymbal2_modes(Q, 1); fb = fb(fb > 0.01);
e0 = abs(fa(1) / f0r - 1);
e1 = abs(fb(1) / f1r - 1);
ok = e0 < 0.03 && e1 < 0.03;
fails = fails + (~ok);
fprintf('washer     : bounce err %.1e, tilt err %.1e vs rigid closed form %s\n', e0, e1, c2v_verdict(ok));

% ---- 3. the paper's linear FEM reference, where comparable ----
% Samejima's Fig. 9 gives linear-FEM eigenfrequencies of the same shell.
% For the clamped center, the (2,1), (3,1), (4,1) modes must agree to a
% few percent ((1,1) is excluded: it is controlled by the support model,
% where this shell's moment-free hole and their 2-D elements differ).
Q = P;
Q.support_clamped = 1;
ref = [20.1, 43.3, 77.7];
got = zeros(1, 3);
for m = 2:4
  f = cymbal2_modes(Q, m);
  f = f(f > 1);
  got(m - 1) = f(1);
end
errF = max(abs(got - ref) ./ ref);
ok = errF < 0.06;
fails = fails + (~ok);
fprintf('paper FEM  : clamped (2,1)(3,1)(4,1) = %.1f %.1f %.1f Hz (ref %.1f %.1f %.1f), err %.1e %s\n', ...
  got(1), got(2), got(3), ref(1), ref(2), ref(3), errF, c2v_verdict(ok));

% ---- 4. transform and update-structure identities ----
% The azimuthal transform must be exactly invertible, and the physical
% derivative pipeline must agree with the per-wavenumber operators: for a
% field cos(m theta) g(r), the physical-space Laplacian must equal the
% per-m radial Laplacian applied to g.
G = cymbal2_grid(P);
eT = max(max(abs(G.F * G.Gi - eye(G.nth))));
gr = exp(-((G.r - 0.12) / 0.05).^2);
m = 7;
Uf = gr * cos(m * G.th');
LapPhys = G.D2 * Uf + G.Rd .* (G.D1 * Uf) + G.R2d .* (Uf * G.T2);
lam2 = (2 * cos(m * G.hth) - 2) / G.hth^2;
op = cymbal2_perm(P, G, lam2);
LapM = (op.Lap * gr) * cos(m * G.th');
r2 = 2:G.nrt - 1;
eL = max(max(abs(LapPhys(r2, :) - LapM(r2, :)))) / max(max(abs(LapM(r2, :))));
ok = eT < 1e-10 && eL < 1e-10;
fails = fails + (~ok);
fprintf('operators  : transform inverse %.1e, physical vs per-m Laplacian %.1e %s\n', ...
  eT, eL, c2v_verdict(ok));

% ---- 5. the assembled update is not amplifying ----
% Every azimuthal column's two-step update must have spectral radius at
% most 1 (the damping makes physical modes strictly contractive; boundary
% constraints contribute zeros). The system builder truncates columns
% that fail; at the defaults none may be.
S = cymbal2_system(P);
nrt = S.G.nrt;
srMax = 0;
for col = [1, 2, S.G.nth / 2, S.G.nth]
  T = [S.AiB(:, :, col), S.AiC(:, :, col); eye(nrt), zeros(nrt)];
  srMax = max(srMax, max(abs(eig(T))));
end
ok = S.nDropped == 0 && srMax <= 1 + 1e-9;
fails = fails + (~ok);
fprintf('stability  : %d truncated columns, sampled spectral radius %.9f %s\n', ...
  S.nDropped, srMax, c2v_verdict(ok));

% ---- 6. stick statics ----
% A static tip load on the clamped-hand stick must deflect by
% F L^3 / (3 E I); the rotational-spring hand adds F L^2 / K_R times L.
MBc = cymbal2_mallet(setfield(P, 'grip_clamped', 1));
MBs = cymbal2_mallet(P);
EI = P.stick_youngs * P.stick_inertia_m4;
L = P.stick_length_m;
Fl = 1.0;
% static balance of the update: at a fixed point, (I - A^-1 B - A^-1 C) d
% equals A^-1 ht^2 K d, so the deflection under a constant tip load F is
% (I - ASiB - ASiC) \ (tipcol * ht^2 * F)
htv = MBc.ht;
dc = ((eye(size(MBc.ASiB, 1)) - MBc.ASiB - MBc.ASiC) \ (MBc.tipcol * htv^2)) * Fl;
ds = ((eye(size(MBs.ASiB, 1)) - MBs.ASiB - MBs.ASiC) \ (MBs.tipcol * htv^2)) * Fl;
wc = dc(MBc.tipi);
ws = ds(MBs.tipi);
wcRef = Fl * L^3 / (3 * EI);
wsRef = wcRef + (Fl * L / P.grip_spring) * L;   % grip rotation F L / K_R times the arm
ec = abs(wc / wcRef - 1);
es = abs(ws / wsRef - 1);
ok = ec < 1e-6 && es < 1e-6;
fails = fails + (~ok);
fprintf('stick      : clamped-hand tip %.3e (ref %.3e), spring-hand %.3e (ref %.3e) %s\n', ...
  wc, wcRef, ws, wsRef, c2v_verdict(ok));

% ---- 7. reproducibility and pinned diagnostics ----
% Two identical renders of a reduced configuration must agree exactly,
% and its diagnostics are pinned. The trajectory of a nonlinear shell is
% sensitive, so the acoustic pins carry percent-level tolerances while
% the contact pins are tight.
Q = P;
Q.sample_rate = 44100;
Q.oversample = 1;
Q.tail_s = 1.5;
[y1, d1] = cymbal2_render(Q);
[y2, d2] = cymbal2_render(Q);
rms1 = sqrt(mean(y1.^2));
rms2 = sqrt(mean(y2.^2));
same = max(abs(y1 - y2)) == 0;
ok = same ...
  && d1.contactEvents == d2.contactEvents ...
  && abs(d1.peakForceN - 311.14) < 4 ...
  && abs(d1.contactVel(1) + 8.223) < 0.08 ...
  && abs(rms1 - 0.0973) < 0.010;
fails = fails + (~ok);
fprintf('pinned     : repro %d, peakF %.2f N (ref 311.14), cvel %.3f (ref -8.223), rms %.4f (ref 0.0973) %s\n', ...
  same, d1.peakForceN, d1.contactVel(1), rms1, c2v_verdict(ok));

% ---- 8. the crash decays ----
n = length(y1);
fs = Q.sample_rate;
head = y1(1:round(0.2 * fs));
tail = y1(n - round(0.2 * fs):n);
dropE = 10 * log10(sum(tail.^2) / sum(head.^2));
ok = dropE < -3 && dropE > -50;
fails = fails + (~ok);
fprintf('decay      : last/first 0.2 s %.1f dB %s\n', dropE, c2v_verdict(ok));

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

function lam = c2v_plate_roots(n, nu)
% Roots of the free-edge frequency equation of a circular plate
lam = [];
prev = c2v_plate_det(n, nu, 0.30);
for x = 0.30:0.01:20
  cur = c2v_plate_det(n, nu, x + 0.01);
  if prev * cur < 0
    lo = x; hi = x + 0.01;
    for it = 1:60
      mid = (lo + hi) / 2;
      if c2v_plate_det(n, nu, lo) * c2v_plate_det(n, nu, mid) < 0
        hi = mid;
      else
        lo = mid;
      end
    end
    lam = [lam, (lo + hi) / 2];
  end
  prev = cur;
end
end

function d = c2v_plate_det(n, nu, x)
J = besselj(n, x); I = besseli(n, x);
Jp = besselj(n - 1, x) - (n / x) * J;
Ip = besseli(n - 1, x) - (n / x) * I;
Jpp = -Jp / x + (n^2 / x^2 - 1) * J;
Ipp = -Ip / x + (n^2 / x^2 + 1) * I;
a11 = Jpp + nu * (Jp / x - n^2 * J / x^2);
a12 = Ipp + nu * (Ip / x - n^2 * I / x^2);
a21 = -Jp - (1 - nu) * (n^2 / x^2) * (Jp - J / x);
a22 = Ip - (1 - nu) * (n^2 / x^2) * (Ip - I / x);
d = a11 * a22 - a12 * a21;
end

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