function [timp, jimp, tauimp, fimp, ximp, ginfo] = maraca_impacts(P)
% MARACA_IMPACTS Stochastic grain-shell impact process for a shake pattern.
%
% The shell holds grain_count grains. A stroke throws them at one wall, where
% each makes bounces_per_grain contacts on average (one arrival plus whatever
% rattle it has before the next reversal), so the expected number of impacts
% in stroke k is grain_count * bounces_per_grain * sqrt(accent_k). The grains
% are not tracked individually: at a few hundred grains in a shell barely
% wider than a few grain diameters, the arrival process is what is audible,
% not which grain arrived. That is the central approximation of this model,
% and it is what makes a few thousand impacts cost milliseconds.
%
% Each stroke k of the gesture produces a burst of independent grain impacts:
%   count   ~ Poisson with that mean (Gaussian approximation, exact enough
%             at hundreds per stroke);
%   times   = stroke onset + slosh_lag + Gamma(3, burst_width/3) delay, the
%             collective slosh of the grain mass arriving at the far wall;
%   speeds  ~ Rayleigh with mean shake_velocity * sqrt(accent_k);
%   points  ~ von Mises-Fisher about the pole the grains are thrown at, which
%             alternates between +x and -x on successive strokes;
%   size    ~ lognormal spread of the grain radius (grain_size_spread).
%
% Each impact is a Hertz contact of a small sphere on the shell wall:
%   k_h  = (4/3) E* sqrt(r_g),  1/E* = (1-nu_g^2)/E_g + (1-nu_s^2)/E_s
%   tau  = 3.218 (m^2 / (k_h^2 v))^{1/5}          (Hertz contact time)
%   J    = (1 + restitution) m v                  (momentum transferred)
%   Fpk  = pi J / (2 tau)                         (half-sine pulse peak)
% The force pulse is modeled as a half sine of duration tau carrying impulse J,
% so J, tau, and Fpk are mutually consistent by construction.
%
% All randomness comes from maraca_rng seeded with P.seed, so a given seed
% reproduces the identical process in the CLI, Node, and the browser.
%
% Returns row vectors sorted by time: timp (s), jimp (N s), tauimp (s),
% fimp (peak N), and ximp (3 x I unit impact positions on the shell).
% ginfo carries the expected-rate envelope and summary statistics for display.

state = mod(floor(abs(P.seed)), 4294967296);
K = P.stroke_count;
lag = P.slosh_lag_ms * 1e-3;
theta = (P.burst_width_ms * 1e-3) / 3;          % Gamma(3, theta): mean = width
spreadRad = P.impact_spread_deg * pi / 180;
kappa = 1 / max(1e-4, spreadRad^2);
duration = P.stroke_count / P.stroke_rate_hz + P.tail_s;

% Hertz constants shared by all impacts (grain Poisson ratio fixed at 0.35)
nuG = 0.35;
Estar = 1 / ((1 - nuG^2) / P.grain_youngs + (1 - P.poisson_ratio^2) / P.youngs_modulus);
rg0 = P.grain_radius_m;
mg0 = (4 / 3) * pi * rg0^3 * P.grain_density;

maxI = ceil(K * P.grain_count * P.bounces_per_grain * sqrt(P.accent_gain) * 1.5) + 64;
timp = zeros(1, maxI);
jimp = zeros(1, maxI);
tauimp = zeros(1, maxI);
fimp = zeros(1, maxI);
vimp = zeros(1, maxI);
ximp = zeros(3, maxI);
I = 0;

counts = zeros(K, 1);
accents = zeros(K, 1);
for k = 1:K
  tk = (k - 1) / P.stroke_rate_hz;
  if mod(k - 1, max(1, round(P.accent_period))) == 0
    a = P.accent_gain;
  else
    a = 1.0;
  end
  accents(k) = a;
  meanCount = P.grain_count * P.bounces_per_grain * sqrt(a);
  [u, state] = maraca_rng(state, 2);
  z = sqrt(-2 * log(u(1))) * cos(2 * pi * u(2));
  nk = max(0, round(meanCount + sqrt(meanCount) * z));
  counts(k) = nk;
  if nk == 0
    continue;
  end
  [u, state] = maraca_rng(state, nk * 8);
  U = reshape(u, nk, 8);

  % times: Gamma(3, theta) as the sum of three exponentials
  tt = tk + lag - theta * log(U(:, 1) .* U(:, 2) .* U(:, 3));
  % speeds: Rayleigh, mean shake_velocity * sqrt(a), capped at 2.5x the mean
  % (a grain cannot exceed a small multiple of the shell speed)
  sigV = P.shake_velocity * sqrt(a) / sqrt(pi / 2);
  vv = sigV * sqrt(-2 * log(U(:, 4)));
  vv = min(vv, 2.5 * P.shake_velocity * sqrt(a));
  % positions: von Mises-Fisher about the active pole (+x odd, -x even strokes)
  w = 1 + log(U(:, 5) * (1 - exp(-2 * kappa)) + exp(-2 * kappa)) / kappa;
  st = sqrt(max(0, 1 - w.^2));
  phi = 2 * pi * U(:, 6);
  pole = 1 - 2 * mod(k - 1, 2);                  % +1, -1, +1, ...
  xs = [pole * w.'; (st .* cos(phi)).'; (st .* sin(phi)).'];
  % grain size: lognormal scale on the radius
  zs = sqrt(-2 * log(U(:, 7))) .* cos(2 * pi * U(:, 8));
  s = exp(P.grain_size_spread * zs);
  mg = mg0 * s.^3;
  kh = (4 / 3) * Estar * sqrt(rg0 * s);

  tau = 3.218 * (mg.^2 ./ (kh.^2 .* vv)).^(1 / 5);
  J = (1 + P.restitution) * mg .* vv;
  Fpk = pi * J ./ (2 * tau);

  keep = find(tt < duration - 2e-3);
  nkeep = length(keep);
  idx = I + 1:I + nkeep;
  timp(idx) = tt(keep).';
  jimp(idx) = J(keep).';
  tauimp(idx) = tau(keep).';
  fimp(idx) = Fpk(keep).';
  vimp(idx) = vv(keep).';
  ximp(:, idx) = xs(:, keep);
  I = I + nkeep;
end
if I == 0
  error('The impact process produced no impacts -- check the gesture parameters.');
end

[timp, order] = sort(timp(1, 1:I));
timp = reshape(timp, 1, I);
jimp = reshape(jimp(order), 1, I);
tauimp = reshape(tauimp(order), 1, I);
fimp = reshape(fimp(order), 1, I);
vimp = reshape(vimp(order), 1, I);
xs = ximp(:, 1:I);
ximp = xs(:, order);

% Expected-rate envelope on a 1 ms grid, for plotting and validation
tg = 0:1e-3:duration;
rate = zeros(size(tg));
for k = 1:K
  tk = (k - 1) / P.stroke_rate_hz + lag;
  tau0 = max(0, tg - tk);
  g = (tau0.^2 / (2 * theta^3)) .* exp(-tau0 / theta);   % Gamma(3,theta) pdf
  rate = rate + P.grain_count * P.bounces_per_grain * sqrt(accents(k)) * g;
end

ginfo.count = I;
ginfo.counts = counts;
ginfo.accents = accents;
ginfo.meanV = sum(vimp) / I;
ginfo.maxV = max(vimp);
ginfo.vimp = vimp;
ginfo.meanTauUs = sum(tauimp) / I * 1e6;
ginfo.rateT = tg;
ginfo.rateHz = rate;
ginfo.durationS = duration;
end