aa03146Maraca shake simulator: MATLAB model and interactive web appmagland 1% Validate the maraca model. There is no external reference implementation,
2% so the checks are (a) closed-form identities the physics must satisfy,
3% (b) exact reproducibility of the seeded impact process, and (c) pinned
4% end-to-end diagnostics from the default configuration, recorded when the
5% model was built and treated as a regression reference since.
6% Run with: numbl run maraca_validate.m
8fails = 0;
10% ---- 1. spherical Hankel functions against closed forms ----
11x = 3.7;
12[h, hp] = maraca_hankel(4, x);
13h0c = sin(x) / x - 1i * (-cos(x) / x);
14h1c = (sin(x) / x^2 - cos(x) / x) - 1i * (-cos(x) / x^2 - sin(x) / x);
15e1 = max(abs(h(1) - h0c), abs(h(2) - h1c));
16d = 1e-6;
17[ha, ~] = maraca_hankel(4, x + d);
18[hb, ~] = maraca_hankel(4, x - d);
19e2 = max(abs(hp - (ha - hb) / (2 * d)));
20xs = 0.5; % small argument: Miller recursion
21[hs, ~] = maraca_hankel(2, xs);
22j2c = (3 / xs^3 - 1 / xs) * sin(xs) - 3 / xs^2 * cos(xs);
23e3 = abs(real(hs(3)) - j2c);
24ok = e1 < 1e-12 && e2 < 1e-6 && e3 < 1e-12;
25fails = fails + (~ok);
26fprintf('hankel closed forms: errs %.1e %.1e %.1e %s\n', e1, e2, e3, verdict(ok));
28% ---- 2. breathing mode: membrane closed form + cavity stiffness ----
29P = maraca_defaults();
30basis = maraca_basis(P);
31nu = P.poisson_ratio;
32cp = sqrt(P.youngs_modulus / (P.density * (1 - nu^2)));
33fdry = sqrt(2 * (1 + nu)) * cp / (2 * pi * P.shell_radius_m) * P.stiffness_trim;
34cav = P.cavity_coupling * 3 * P.air_density * P.speed_of_sound^2 ...
35 / (P.density * P.shell_thickness_m * P.shell_radius_m);
36fwet = sqrt((2 * pi * fdry)^2 + cav) / (2 * pi);
37ok = abs(basis.breathingHz - fwet) < 1e-6 * fwet && all(basis.zeta > 0);
38fails = fails + (~ok);
39fprintf('breathing mode: %.1f Hz (closed form %.1f) %s\n', basis.breathingHz, fwet, verdict(ok));
41% ---- 3. seeded impact process is exactly reproducible ----
42[timp, jimp, tauimp, fimp, ximp, ginfo] = maraca_impacts(P);
43ok = ginfo.count == 3628 ...
44 && abs(ginfo.meanV - 2.462374) < 1e-5 ...
45 && abs(sum(jimp) - 3.586588617e-01) < 1e-9 ...
46 && abs(timp(1) - 0.016172179) < 1e-8;
47fails = fails + (~ok);
48fprintf('impact process (seed %d): %d impacts, mean v %.4f m/s, sum J %.6f N s %s\n', ...
49 P.seed, ginfo.count, ginfo.meanV, sum(jimp), verdict(ok));
51% ---- 4. Hertz contact-time scaling tau ~ v^(-1/5) ----
52Q = P;
53Q.grain_count = 150;
54[~, ~, tau1, ~, ~, g1] = maraca_impacts(Q);
55Q.shake_velocity = 2 * P.shake_velocity;
56[~, ~, tau2, ~, ~, g2] = maraca_impacts(Q);
57ratio = (sum(tau2) / length(tau2)) / (sum(tau1) / length(tau1));
58ok = abs(ratio - 2^(-1 / 5)) < 0.03; % same seed, so same draws
59fails = fails + (~ok);
60fprintf('hertz scaling: mean tau ratio %.4f (target %.4f) %s\n', ratio, 2^(-1 / 5), verdict(ok));
62% ---- 5. far-field falloff: -6 dB per doubling of distance ----
63% RMS, not peak: the per-mode radiation gains scale exactly as 1/r, but the
64% near-field phase terms of h_l shift the interference at the peak sample by
65% a few tenths of a dB between distances, which is physics rather than error.
66Q = P;
67Q.normalize = 0;
68Q.listener_x_m = 0; Q.listener_y_m = 0; Q.listener_z_m = 1.0;
69[y1, ~, ~] = maraca_render(Q);
70Q.listener_z_m = 2.0;
71[y2, ~, ~] = maraca_render(Q);
72drop = 10 * log10(sum(y2.^2) / sum(y1.^2));
73ok = abs(drop + 6.02) < 0.15;
74fails = fails + (~ok);
75fprintf('falloff 1 m -> 2 m: %.2f dB rms (target -6.02) %s\n', drop, verdict(ok));
77% ---- 6. pinned end-to-end diagnostics at the defaults ----
78[y, basis, dinfo] = maraca_render(P);
79ok = abs(dinfo.peakPa - 14.935221) < 0.02 ...
80 && abs(basis.f(1) - 6445.5174) < 0.1 ...
81 && basis.n == 11;
82fails = fails + (~ok);
83fprintf('pinned render: peak %.3f Pa (target 14.935), mode1 %.1f Hz (target 6445.5), %d modes %s\n', ...
84 dinfo.peakPa, basis.f(1), basis.n, verdict(ok));
86% ---- 7. the tail decays to silence ----
87fs = P.sample_rate;
88n = length(y);
89tail = y(n - round(0.08 * fs):n);
90tailRms = sqrt(sum(tail.^2) / length(tail));
91ok = tailRms < 0.02 * max(abs(y));
92fails = fails + (~ok);
93fprintf('tail decay: rms %.2e of peak %s\n', tailRms / max(abs(y)), verdict(ok));
95if fails > 0
96 fprintf('%d CHECK(S) FAILED\n', fails);
97else
98 fprintf('all checks passed\n');
99end
101function s = verdict(ok)
102if ok
103 s = 'ok';
104else
105 s = 'MISMATCH';
106end
107end