function [tstr, vstr, jstr, taustr, fstr, sstr, sinfo] = cymbal_strikes(P, basis)
% CYMBAL_STRIKES Stick contacts: a single strike, or a seeded repeated figure.
%
% strike_count = 1 gives one deterministic strike of strike_velocity at
% strike_r at t = 0; no random numbers are drawn. strike_count > 1 repeats it
% at strike_rate_hz, alternating between strike_r and roll_r2 (a ride pattern
% moving between the bow and the edge), with Gaussian timing jitter
% (timing_jitter_ms) and lognormal velocity spread (velocity_spread). The
% first hit stays at t = 0 so the onset is fixed. All randomness comes from
% cymbal_rng seeded with P.seed, so a given seed reproduces the identical
% pattern in the CLI, in Node, and in the browser. Every contact sits on the
% meridian theta = 0; see cymbal_basis.m for why.
%
% Each contact is a Hertz impact of the stick tip on the plate:
%   E*   = 1 / ((1 - nu_t^2) / E_t + (1 - nu_p^2) / E_p)
%   R_eff= 1 / sqrt((1/r_tip + 1/R_s)(1/r_tip + 1/R_theta))   (sphere on a
%          doubly curved surface, at the local principal curvatures)
%   k_h  = (4/3) E* sqrt(R_eff)
%   tau  = 3.218 (m*^2 / (k_h^2 v))^{1/5}                     (Hertz contact time)
%   J    = (1 + restitution) m* v
%   Fpk  = pi J / (2 tau)                                     (half-sine peak)
%
% The reduced mass m* pairs the stick with the plate mass engaged during the
% contact. Here the sibling triangle's flexural-group-velocity estimate is
% replaced by something exact: the driving-point impedance of an infinite
% plate is a pure resistance, Z = 8 sqrt(D rho t), independent of frequency,
% so the mass the stick sees during a pulse of duration tau is simply Z tau. A
% short fixed-point iteration closes the loop, since tau depends on m* which
% depends on tau. What this captures is that a heavy stick, a soft tip or a
% slow stroke all lengthen the contact and darken the strike; what it leaves
% out is the millisecond detail of the contact itself, which is prescribed
% rather than simulated.
%
% Returns row vectors sorted by time: tstr (s), vstr (m/s), jstr (N s),
% taustr (s), fstr (peak N), sstr (arclength fraction of the contact), and
% sinfo with summary statistics.

K = max(1, round(P.strike_count));
rate = max(1e-3, P.strike_rate_hz);
duration = (K - 1) / rate + P.tail_s;

E = P.youngs_modulus;
nu = P.poisson_ratio;
rho = P.density;
nut = 0.35;                          % stick tip Poisson ratio, fixed
Estar = 1 / ((1 - nut^2) / P.stick_youngs + (1 - nu^2) / E);
rt = max(1e-4, P.tip_radius_m);

state = mod(floor(abs(P.seed)), 4294967296);
z1 = zeros(K, 1);
z2 = zeros(K, 1);
if K > 1
  [u, state] = cymbal_rng(state, 4 * K);
  U = reshape(u, K, 4);
  z1 = sqrt(-2 * log(U(:, 1))) .* cos(2 * pi * U(:, 2));
  z2 = sqrt(-2 * log(U(:, 3))) .* cos(2 * pi * U(:, 4));
end

tstr = zeros(1, K);
vstr = zeros(1, K);
jstr = zeros(1, K);
taustr = zeros(1, K);
fstr = zeros(1, K);
sstr = zeros(1, K);
for i = 1:K
  t = (i - 1) / rate;
  if i > 1
    t = max(0, t + P.timing_jitter_ms * 1e-3 * z1(i));
  end
  v = P.strike_velocity;
  if K > 1
    v = v * exp(P.velocity_spread * z2(i));
  end
  v = max(0.01, v);
  if mod(i, 2) == 1
    sfrac = P.strike_r;
  else
    sfrac = P.roll_r2;
  end
  sfrac = min(0.999, max(0.001, sfrac));

  % local plate properties at the contact
  [~, ~, thk, ~, ~, cs, ct] = cymbal_at(basis, sfrac);
  Rrel = 1 / sqrt(max(1e-9, (1 / rt + cs) * (1 / rt + ct)));
  kh = (4 / 3) * Estar * sqrt(Rrel);
  D = E * thk^3 / (12 * (1 - nu^2));
  Zp = 8 * sqrt(D * rho * thk);      % driving-point resistance of the plate

  tau = 1e-4;
  mstar = P.stick_mass_kg;
  for it = 1:6
    mloc = Zp * tau;
    mstar = P.stick_mass_kg * mloc / (P.stick_mass_kg + mloc);
    tau = 3.218 * (mstar^2 / (kh^2 * v))^(1 / 5);
  end

  tstr(i) = t;
  vstr(i) = v;
  jstr(i) = (1 + P.restitution) * mstar * v;
  taustr(i) = tau;
  fstr(i) = pi * jstr(i) / (2 * tau);
  sstr(i) = sfrac;
end

[tstr, order] = sort(tstr);
tstr = reshape(tstr, 1, K);
vstr = reshape(vstr(order), 1, K);
jstr = reshape(jstr(order), 1, K);
taustr = reshape(taustr(order), 1, K);
fstr = reshape(fstr(order), 1, K);
sstr = reshape(sstr(order), 1, K);

sinfo.count = K;
sinfo.meanV = sum(vstr) / K;
sinfo.maxV = max(vstr);
sinfo.meanTauUs = sum(taustr) / K * 1e6;
sinfo.peakForceN = max(fstr);
sinfo.durationS = duration;
end