% Validate the cymbal model. There is no external reference implementation, so
% the checks are (a) closed-form identities the physics must satisfy, (b)
% structural identities of the shell eigensolve and of the radiation integral,
% (c) exact reproducibility of the seeded strike 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 cymbal_validate.m

fails = 0;

% ---- 1. flat free circular plate against the closed form ----
% The same shell routine that builds the cymbal, run on a flat meridian of
% constant thickness with a pinhole at the centre, must reproduce the
% completely free circular plate: the roots of the 2 x 2 determinant built
% from J_n and I_n under vanishing radial moment and Kelvin-Kirchhoff shear,
% and exactly two rigid-body modes at n = 0 (axial translation and rotation
% about the axis), two at n = 1 (translation and tilt) and none above.
P = cymbal_defaults();
Q = P;
Q.bow_height_m = 0; Q.bell_height_m = 0;
Q.thickness_bell_m = 0.001; Q.thickness_rim_m = 0.001;
Q.hole_radius_m = 0.002; Q.radius_m = 0.2; Q.n_elements = 40;
Gf = cymbal_profile(Q);
Dp = Q.youngs_modulus * 0.001^3 / (12 * (1 - Q.poisson_ratio^2));
scal = sqrt(Dp / (Q.density * 0.001)) / Q.radius_m^2 / (2 * pi);
Lr = Gf.S / (Gf.nn - 1);
errP = 0;
rigid = zeros(1, 5);
for n = 0:4
  f = cv_freqs(Gf, n, Q, Lr);
  rigid(n + 1) = sum(f < 3);
  fe = f(f > 3);
  lam = cv_plate_roots(n, Q.poisson_ratio);
  fc = lam.^2 * scal;
  nc = min(4, min(length(fe), length(fc)));
  errP = max(errP, max(abs(fe(1:nc).' - fc(1:nc)) ./ fc(1:nc)));
end
ok = errP < 3e-3 && rigid(1) == 2 && rigid(2) == 2 && sum(rigid(3:5)) == 0;
fails = fails + (~ok);
fprintf('flat plate : max freq err %.1e, rigid modes %d %d %d %d %d %s\n', ...
  errP, rigid(1), rigid(2), rigid(3), rigid(4), rigid(5), cv_verdict(ok));

% ---- 2. eigensolve identities on the cymbal 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).
G = cymbal_profile(P);
Lref = G.S / (G.nn - 1);
[K, M] = cymbal_shell(G, 5, P.youngs_modulus, P.poisson_ratio, P.density, Lref);
Rc = chol(M);
Am = (Rc') \ (K / 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(length(om2)))));
eDiag = max(max(abs(Phi' * K * Phi - diag(om2)))) / max(om2);
ok = eOrth < 1e-8 && eDiag < 1e-10;
fails = fails + (~ok);
fprintf('eigensolve : orthonormality %.1e, diagonalization %.1e %s\n', eOrth, eDiag, cv_verdict(ok));

% ---- 3. the Bessel far field against a direct sum over the dipole sheet ----
% cymbal_basis reduces the azimuthal integral of the far-field kernel to
% Bessel functions in order to get the radiated power; that reduction must
% agree with brute force summation over a finely discretized dipole sheet.
Q = P;
Q.n_diameters = 4; Q.max_mode_hz = 1200;
b4 = cymbal_basis(Q);
errF = 0;
for m = [1 3 5 9]
  [pd, pb] = cv_farfield(Q, b4, m, 200, 720, 0.6, 0.4);
  errF = max(errF, abs(pb / pd - 1));
end
ok = errF < 0.01;
fails = fails + (~ok);
fprintf('far field  : Bessel vs direct dipole sum, max err %.2e %s\n', errF, cv_verdict(ok));

% ---- 4. on the axis, every mode with a nodal diameter cancels exactly ----
% A cos(n theta) source on a surface of revolution has no monopole, dipole or
% any lower multipole than 2^n, and on the symmetry axis all of them vanish.
% This is the sharpest available check that the azimuthal sum in the listener
% transfer is assembled correctly. It is a statement about the IDEAL plate, so
% the hammering is switched off: hybridizing the modes is precisely what gives
% the high-order ones a way out, and with it on they radiate on the axis like
% anything else.
Q = P;
Q.hammer_mix = 0; Q.hammer_detune = 0;
Q.n_diameters = 6; Q.max_mode_hz = 1500;
Q.listener_x_m = 0; Q.listener_y_m = 0; Q.listener_z_m = 0.6;
b6 = cymbal_basis(Q);
[tstr, ~, jstr, taustr, fstr, sstr, sinf] = cymbal_strikes(Q, b6);
[env, tfr, DV, kf] = cymbal_cascade(Q, b6, tstr, jstr, taustr, sstr, 0.05);
[~, ri] = cymbal_radiate(Q, b6, env, tfr, DV, kf, tstr, taustr, fstr, sstr, Q.sample_rate);
g0 = max(ri.gainAbs(b6.ndia == 0));
g1 = max(ri.gainAbs(b6.ndia >= 1));
ok = g1 / g0 < 1e-10;
fails = fails + (~ok);
fprintf('on axis    : |G| of n>=1 modes is %.1e of the n=0 modes %s\n', g1 / g0, cv_verdict(ok));

% ---- 4b. the hammering mixing is exactly orthogonal ----
% The hybridization is a product of Givens rotations, so it must be orthogonal
% to round-off. That is what makes it energy conserving: it redistributes the
% stroke among the modes and the radiation among them again, without inventing
% or destroying any.
bh = cymbal_basis(P);
eMix = max(max(abs(bh.mix * bh.mix.' - eye(bh.n))));
spread = sum(bh.mix.^2 > 0.01, 2);
ok = eMix < 1e-12 && median(spread) > 2;
fails = fails + (~ok);
fprintf('hammering  : orthogonality %.1e, median %d ideal modes per hybrid %s\n', ...
  eMix, round(median(spread)), cv_verdict(ok));

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

% ---- 6. the cascade conserves energy, and switches off cleanly ----
% With the damping removed, the top of the ladder reflecting and the
% fluctuation off, the transfer scheme is a Markov generator on the modal
% energies and must conserve their sum to round-off. With cascade_gain = 0
% nothing is exchanged, so the fluctuation freezes at zero on its own and
% every envelope must collapse onto the closed-form exp(-zeta om t) of an
% uncoupled damped oscillator, fluctuation still switched on.
Q = P;
Q.n_diameters = 6; Q.max_mode_hz = 1500; Q.cascade_leak = 0; Q.turbulence = 0;
bz = b6;
bz.zeta = zeros(bz.n, 1);
[tstr, ~, jstr, taustr, ~, sstr] = cymbal_strikes(Q, bz);
[envC, ~, DVc] = cymbal_cascade(Q, bz, tstr, jstr, taustr, sstr, 0.4);
Etot = sum(envC.^2 / 2, 1);
eCons = max(abs(Etot / Etot(1) - 1));
Q.cascade_gain = 0;
Q.turbulence = 0.9;   % explicitly on, so the freeze is actually tested
[envL, tfrL] = cymbal_cascade(Q, b6, tstr, jstr, taustr, sstr, 0.4);
pred = envL(:, 1) * exp(-(b6.zeta .* b6.omega).' * 0 + 0);
pred = (envL(:, 1) * ones(1, length(tfrL))) .* exp(-(b6.zeta .* b6.omega) * tfrL);
sel = envL(:, 1) > 1e-12;
eLin = max(max(abs(envL(sel, :) ./ pred(sel, :) - 1)));
ok = eCons < 1e-10 && eLin < 1e-9;
fails = fails + (~ok);
fprintf('cascade    : energy drift %.1e, linear limit err %.1e %s\n', eCons, eLin, cv_verdict(ok));

% ---- 7. 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.
Q = P;
Q.normalize = 0; Q.n_diameters = 12; Q.max_mode_hz = 3000;
Q.tail_s = 0.4; Q.oversample = 1; Q.n_azimuth = 96;
Q.listener_x_m = 0.8; Q.listener_y_m = 0.5; Q.listener_z_m = 1.8;
y1 = cymbal_render(Q);
Q.listener_x_m = 1.6; Q.listener_y_m = 1.0; Q.listener_z_m = 3.6;
y2 = cymbal_render(Q);
drop = 10 * log10(sum(y2.^2) / sum(y1.^2));
ok = abs(drop + 6.02) < 0.25;
fails = fails + (~ok);
fprintf('falloff    : 2 m -> 4 m %.2f dB rms (target -6.02) %s\n', drop, cv_verdict(ok));

% ---- 8. the seeded ride figure is exactly reproducible ----
Q = P;
Q.strike_count = 8;
[t1, v1, j1, ~, ~, s1, g1] = cymbal_strikes(Q, b6);
[t2, v2, ~, ~, ~, s2] = cymbal_strikes(Q, b6);
same = max(abs(t1 - t2)) == 0 && max(abs(v1 - v2)) == 0 && max(abs(s1 - s2)) == 0;
ok = same ...
  && g1.count == 8 ...
  && abs(g1.meanV - 3.111518886) < 2e-6 ...
  && abs(sum(j1) - 0.150818144) < 1e-9 ...
  && abs(t1(2) - 0.159427523) < 1e-9;
fails = fails + (~ok);
fprintf('ride (seed %d): %d strikes, mean v %.6f m/s, sum J %.9f N s %s\n', ...
  Q.seed, g1.count, g1.meanV, sum(j1), cv_verdict(ok));

% ---- 9. pinned end-to-end diagnostics at the defaults ----
% The rms is pinned alongside the peak deliberately. The peak of a sum of
% several hundred oscillators is a fragile thing to compare across machines,
% where the rms is not, so a mismatch on the peak alone says the phases have
% moved and a mismatch on both says the physics has.
[y, basis, dinfo] = cymbal_render(P);
yrms = sqrt(sum(y.^2) / length(y));
ok = abs(dinfo.peakPa - 13.6155) < 0.02 ...
  && abs(yrms - 0.04652826) < 2e-6 ...
  && abs(basis.f(1) - 49.4130) < 0.02 ...
  && basis.n == 378;
fails = fails + (~ok);
fprintf('pinned     : peak %.4f Pa (target %.4f), rms %.8f (target %.8f), mode1 %.3f Hz, %d modes %s\n', ...
  dinfo.peakPa, 13.6155, yrms, 0.04652826, basis.f(1), basis.n, cv_verdict(ok));

% ---- 10. the wash decays ----
% A cymbal rings for seconds, so the tail is not silent, but the energy in the
% last half second must sit well below that of the opening half second.
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 < -20 && dropE > -80;
fails = fails + (~ok);
fprintf('decay      : last/first half second %.1f dB %s\n', dropE, cv_verdict(ok));

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

function f = cv_freqs(G, n, P, Lref)
[K, M] = cymbal_shell(G, n, P.youngs_modulus, P.poisson_ratio, P.density, Lref);
Rc = chol(M);
Am = (Rc') \ (K / Rc);
Am = (Am + Am') / 2;
[~, Sv, ~] = svd(Am);
f = sort(sqrt(diag(Sv))) / (2 * pi);
end

function lam = cv_plate_roots(n, nu)
% Roots of the free-edge frequency equation of a circular plate, bracketed on
% a coarse scan and bisected.
lam = [];
prev = cv_plate_det(n, nu, 0.30);
for x = 0.30:0.01:26
  cur = cv_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 cv_plate_det(n, nu, lo) * cv_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 = cv_plate_det(n, nu, x)
% W = A J_n(lam r/a) + B I_n(lam r/a); the free edge sets the radial moment
% and the Kelvin-Kirchhoff shear to zero, and the 2 x 2 determinant vanishes.
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 [pd, pb] = cv_farfield(P, b, m, R, nth, tho, pho)
% Far-field pressure of one mode, per unit modal velocity, computed by direct
% summation over a theta-discretized dipole sheet and by the Bessel reduction.
c = P.speed_of_sound;
rho0 = P.air_density;
om = b.omega(m);
k = om / c;
n = b.ndia(m);
ld = b.ldip(m);
rh = [sin(tho) * cos(pho), sin(tho) * sin(pho), cos(tho)];
th = ((0:nth - 1) + 0.5) / nth * 2 * pi;
dth = 2 * pi / nth;
acc = 0;
for e = 1:length(b.surfR)
  x = b.surfR(e) * cos(th);
  y = b.surfR(e) * sin(th);
  z = b.surfZ(e) * ones(1, nth);
  nx = -b.surfDz(e) * cos(th);
  ny = -b.surfDz(e) * sin(th);
  nz = b.surfDr(e) * ones(1, nth);
  dx = R * rh(1) - x;
  dy = R * rh(2) - y;
  dz = R * rh(3) - z;
  rr = sqrt(dx.^2 + dy.^2 + dz.^2);
  proj = (nx .* dx + ny .* dy + nz .* dz) ./ rr * b.SW(e, m) .* cos(n * th);
  V = ld * b.surfR(e) * b.surfLen(e) * dth;
  acc = acc + sum((1i * om * rho0 / (4 * pi)) * V * proj .* (1i * k + 1 ./ rr) .* exp(-1i * k * rr) ./ rr);
end
pd = abs(acc);
mu = cos(tho);
sth = sin(tho);
al = k * b.surfR * sth;
Jn = besselj(n, al);
Jp1 = besselj(n + 1, al);
if n == 0
  Jm1 = -Jp1;
else
  Jm1 = besselj(n - 1, al);
end
ip = 1i^n;
Tk = ip * (2 * pi * Jn .* (b.surfDr * mu) - 1i * pi * (Jp1 - Jm1) .* (b.surfDz * sth));
An = sum(ld * b.surfR .* b.surfLen .* b.SW(:, m) .* Tk .* exp(1i * (k * b.surfZ) * mu));
pb = abs(-(om^2 * rho0 / (4 * pi * c * R)) * exp(-1i * k * R) * An * cos(n * pho));
end

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