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

fails = 0;

% ---- 1. spherical Hankel functions against closed forms ----
x = 3.7;
[h, hp] = maraca_hankel(4, x);
h0c = sin(x) / x - 1i * (-cos(x) / x);
h1c = (sin(x) / x^2 - cos(x) / x) - 1i * (-cos(x) / x^2 - sin(x) / x);
e1 = max(abs(h(1) - h0c), abs(h(2) - h1c));
d = 1e-6;
[ha, ~] = maraca_hankel(4, x + d);
[hb, ~] = maraca_hankel(4, x - d);
e2 = max(abs(hp - (ha - hb) / (2 * d)));
xs = 0.5;                                % small argument: Miller recursion
[hs, ~] = maraca_hankel(2, xs);
j2c = (3 / xs^3 - 1 / xs) * sin(xs) - 3 / xs^2 * cos(xs);
e3 = abs(real(hs(3)) - j2c);
ok = e1 < 1e-12 && e2 < 1e-6 && e3 < 1e-12;
fails = fails + (~ok);
fprintf('hankel closed forms: errs %.1e %.1e %.1e %s\n', e1, e2, e3, verdict(ok));

% ---- 2. breathing mode: membrane closed form + cavity stiffness ----
P = maraca_defaults();
basis = maraca_basis(P);
nu = P.poisson_ratio;
cp = sqrt(P.youngs_modulus / (P.density * (1 - nu^2)));
fdry = sqrt(2 * (1 + nu)) * cp / (2 * pi * P.shell_radius_m) * P.stiffness_trim;
cav = P.cavity_coupling * 3 * P.air_density * P.speed_of_sound^2 ...
  / (P.density * P.shell_thickness_m * P.shell_radius_m);
fwet = sqrt((2 * pi * fdry)^2 + cav) / (2 * pi);
ok = abs(basis.breathingHz - fwet) < 1e-6 * fwet && all(basis.zeta > 0);
fails = fails + (~ok);
fprintf('breathing mode: %.1f Hz (closed form %.1f) %s\n', basis.breathingHz, fwet, verdict(ok));

% ---- 3. seeded impact process is exactly reproducible ----
[timp, jimp, tauimp, fimp, ximp, ginfo] = maraca_impacts(P);
ok = ginfo.count == 3628 ...
  && abs(ginfo.meanV - 2.462374) < 1e-5 ...
  && abs(sum(jimp) - 3.586588617e-01) < 1e-9 ...
  && abs(timp(1) - 0.016172179) < 1e-8;
fails = fails + (~ok);
fprintf('impact process (seed %d): %d impacts, mean v %.4f m/s, sum J %.6f N s %s\n', ...
  P.seed, ginfo.count, ginfo.meanV, sum(jimp), verdict(ok));

% ---- 4. Hertz contact-time scaling tau ~ v^(-1/5) ----
Q = P;
Q.grain_count = 150;
[~, ~, tau1, ~, ~, g1] = maraca_impacts(Q);
Q.shake_velocity = 2 * P.shake_velocity;
[~, ~, tau2, ~, ~, g2] = maraca_impacts(Q);
ratio = (sum(tau2) / length(tau2)) / (sum(tau1) / length(tau1));
ok = abs(ratio - 2^(-1 / 5)) < 0.03;     % same seed, so same draws
fails = fails + (~ok);
fprintf('hertz scaling: mean tau ratio %.4f (target %.4f) %s\n', ratio, 2^(-1 / 5), verdict(ok));

% ---- 5. far-field falloff: -6 dB per doubling of distance ----
% RMS, not peak: the per-mode radiation gains scale exactly as 1/r, but the
% near-field phase terms of h_l shift the interference at the peak sample by
% a few tenths of a dB between distances, which is physics rather than error.
Q = P;
Q.normalize = 0;
Q.listener_x_m = 0; Q.listener_y_m = 0; Q.listener_z_m = 1.0;
[y1, ~, ~] = maraca_render(Q);
Q.listener_z_m = 2.0;
[y2, ~, ~] = maraca_render(Q);
drop = 10 * log10(sum(y2.^2) / sum(y1.^2));
ok = abs(drop + 6.02) < 0.15;
fails = fails + (~ok);
fprintf('falloff 1 m -> 2 m: %.2f dB rms (target -6.02) %s\n', drop, verdict(ok));

% ---- 6. pinned end-to-end diagnostics at the defaults ----
[y, basis, dinfo] = maraca_render(P);
ok = abs(dinfo.peakPa - 14.935221) < 0.02 ...
  && abs(basis.f(1) - 6445.5174) < 0.1 ...
  && basis.n == 11;
fails = fails + (~ok);
fprintf('pinned render: peak %.3f Pa (target 14.935), mode1 %.1f Hz (target 6445.5), %d modes %s\n', ...
  dinfo.peakPa, basis.f(1), basis.n, verdict(ok));

% ---- 7. the tail decays to silence ----
fs = P.sample_rate;
n = length(y);
tail = y(n - round(0.08 * fs):n);
tailRms = sqrt(sum(tail.^2) / length(tail));
ok = tailRms < 0.02 * max(abs(y));
fails = fails + (~ok);
fprintf('tail decay: rms %.2e of peak %s\n', tailRms / max(abs(y)), verdict(ok));

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

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