aa03146Maraca shake simulator: MATLAB model and interactive web appmagland 1function [timp, jimp, tauimp, fimp, ximp, ginfo] = maraca_impacts(P)
2% MARACA_IMPACTS Stochastic grain-shell impact process for a shake pattern.
3%
4% The shell holds grain_count grains. A stroke throws them at one wall, where
5% each makes bounces_per_grain contacts on average (one arrival plus whatever
6% rattle it has before the next reversal), so the expected number of impacts
7% in stroke k is grain_count * bounces_per_grain * sqrt(accent_k). The grains
8% are not tracked individually: at a few hundred grains in a shell barely
9% wider than a few grain diameters, the arrival process is what is audible,
10% not which grain arrived. That is the central approximation of this model,
11% and it is what makes a few thousand impacts cost milliseconds.
12%
13% Each stroke k of the gesture produces a burst of independent grain impacts:
14% count ~ Poisson with that mean (Gaussian approximation, exact enough
15% at hundreds per stroke);
16% times = stroke onset + slosh_lag + Gamma(3, burst_width/3) delay, the
17% collective slosh of the grain mass arriving at the far wall;
18% speeds ~ Rayleigh with mean shake_velocity * sqrt(accent_k);
19% points ~ von Mises-Fisher about the pole the grains are thrown at, which
20% alternates between +x and -x on successive strokes;
21% size ~ lognormal spread of the grain radius (grain_size_spread).
22%
23% Each impact is a Hertz contact of a small sphere on the shell wall:
24% k_h = (4/3) E* sqrt(r_g), 1/E* = (1-nu_g^2)/E_g + (1-nu_s^2)/E_s
25% tau = 3.218 (m^2 / (k_h^2 v))^{1/5} (Hertz contact time)
26% J = (1 + restitution) m v (momentum transferred)
27% Fpk = pi J / (2 tau) (half-sine pulse peak)
28% The force pulse is modeled as a half sine of duration tau carrying impulse J,
29% so J, tau, and Fpk are mutually consistent by construction.
30%
31% All randomness comes from maraca_rng seeded with P.seed, so a given seed
32% reproduces the identical process in the CLI, Node, and the browser.
33%
34% Returns row vectors sorted by time: timp (s), jimp (N s), tauimp (s),
35% fimp (peak N), and ximp (3 x I unit impact positions on the shell).
36% ginfo carries the expected-rate envelope and summary statistics for display.
38state = mod(floor(abs(P.seed)), 4294967296);
39K = P.stroke_count;
40lag = P.slosh_lag_ms * 1e-3;
41theta = (P.burst_width_ms * 1e-3) / 3; % Gamma(3, theta): mean = width
42spreadRad = P.impact_spread_deg * pi / 180;
43kappa = 1 / max(1e-4, spreadRad^2);
44duration = P.stroke_count / P.stroke_rate_hz + P.tail_s;
46% Hertz constants shared by all impacts (grain Poisson ratio fixed at 0.35)
47nuG = 0.35;
48Estar = 1 / ((1 - nuG^2) / P.grain_youngs + (1 - P.poisson_ratio^2) / P.youngs_modulus);
49rg0 = P.grain_radius_m;
50mg0 = (4 / 3) * pi * rg0^3 * P.grain_density;
52maxI = ceil(K * P.grain_count * P.bounces_per_grain * sqrt(P.accent_gain) * 1.5) + 64;
53timp = zeros(1, maxI);
54jimp = zeros(1, maxI);
55tauimp = zeros(1, maxI);
56fimp = zeros(1, maxI);
57vimp = zeros(1, maxI);
58ximp = zeros(3, maxI);
59I = 0;
61counts = zeros(K, 1);
62accents = zeros(K, 1);
63for k = 1:K
64 tk = (k - 1) / P.stroke_rate_hz;
65 if mod(k - 1, max(1, round(P.accent_period))) == 0
66 a = P.accent_gain;
67 else
68 a = 1.0;
69 end
70 accents(k) = a;
71 meanCount = P.grain_count * P.bounces_per_grain * sqrt(a);
72 [u, state] = maraca_rng(state, 2);
73 z = sqrt(-2 * log(u(1))) * cos(2 * pi * u(2));
74 nk = max(0, round(meanCount + sqrt(meanCount) * z));
75 counts(k) = nk;
76 if nk == 0
77 continue;
78 end
79 [u, state] = maraca_rng(state, nk * 8);
80 U = reshape(u, nk, 8);
82 % times: Gamma(3, theta) as the sum of three exponentials
83 tt = tk + lag - theta * log(U(:, 1) .* U(:, 2) .* U(:, 3));
84 % speeds: Rayleigh, mean shake_velocity * sqrt(a), capped at 2.5x the mean
85 % (a grain cannot exceed a small multiple of the shell speed)
86 sigV = P.shake_velocity * sqrt(a) / sqrt(pi / 2);
87 vv = sigV * sqrt(-2 * log(U(:, 4)));
88 vv = min(vv, 2.5 * P.shake_velocity * sqrt(a));
89 % positions: von Mises-Fisher about the active pole (+x odd, -x even strokes)
90 w = 1 + log(U(:, 5) * (1 - exp(-2 * kappa)) + exp(-2 * kappa)) / kappa;
91 st = sqrt(max(0, 1 - w.^2));
92 phi = 2 * pi * U(:, 6);
93 pole = 1 - 2 * mod(k - 1, 2); % +1, -1, +1, ...
94 xs = [pole * w.'; (st .* cos(phi)).'; (st .* sin(phi)).'];
95 % grain size: lognormal scale on the radius
96 zs = sqrt(-2 * log(U(:, 7))) .* cos(2 * pi * U(:, 8));
97 s = exp(P.grain_size_spread * zs);
98 mg = mg0 * s.^3;
99 kh = (4 / 3) * Estar * sqrt(rg0 * s);
101 tau = 3.218 * (mg.^2 ./ (kh.^2 .* vv)).^(1 / 5);
102 J = (1 + P.restitution) * mg .* vv;
103 Fpk = pi * J ./ (2 * tau);
105 keep = find(tt < duration - 2e-3);
106 nkeep = length(keep);
107 idx = I + 1:I + nkeep;
108 timp(idx) = tt(keep).';
109 jimp(idx) = J(keep).';
110 tauimp(idx) = tau(keep).';
111 fimp(idx) = Fpk(keep).';
112 vimp(idx) = vv(keep).';
113 ximp(:, idx) = xs(:, keep);
114 I = I + nkeep;
115end
116if I == 0
117 error('The impact process produced no impacts -- check the gesture parameters.');
118end
120[timp, order] = sort(timp(1, 1:I));
121timp = reshape(timp, 1, I);
122jimp = reshape(jimp(order), 1, I);
123tauimp = reshape(tauimp(order), 1, I);
124fimp = reshape(fimp(order), 1, I);
125vimp = reshape(vimp(order), 1, I);
126xs = ximp(:, 1:I);
127ximp = xs(:, order);
129% Expected-rate envelope on a 1 ms grid, for plotting and validation
130tg = 0:1e-3:duration;
131rate = zeros(size(tg));
132for k = 1:K
133 tk = (k - 1) / P.stroke_rate_hz + lag;
134 tau0 = max(0, tg - tk);
135 g = (tau0.^2 / (2 * theta^3)) .* exp(-tau0 / theta); % Gamma(3,theta) pdf
136 rate = rate + P.grain_count * P.bounces_per_grain * sqrt(accents(k)) * g;
137end
139ginfo.count = I;
140ginfo.counts = counts;
141ginfo.accents = accents;
142ginfo.meanV = sum(vimp) / I;
143ginfo.maxV = max(vimp);
144ginfo.vimp = vimp;
145ginfo.meanTauUs = sum(tauimp) / I * 1e6;
146ginfo.rateT = tg;
147ginfo.rateHz = rate;
148ginfo.durationS = duration;
149end