function [acc, dinfo, fhist] = cymbal2_sim(P, S, MB)
% CYMBAL2_SIM Time-domain integration of the full system: von Karman
% shallow spherical shell (implicit in the stiff linear terms, explicit in
% the nonlinearity), washer support built into the boundary rows, mallet
% head + stick coupled through the Hertz contact. Returns the normal
% acceleration of the shell at the pickup node at the simulation rate,
% diagnostics, and the contact force history.

G = S.G;
nrt = G.nrt;
Nr = G.Nr;
nth = G.nth;
ht = S.ht;
rhoh = S.rhoh;
Rs = S.Rs;
EhRs = S.Eh / Rs;
Eh2 = S.Eh / 2;
ls = S.ls; js = S.js; lp = S.lp; jp = S.jp;
fsc = S.forceScale;
nlr0 = S.eomRow0;

K = max(1, round(P.strike_count));

% Stroke schedule; per-stroke peak hand speeds with a seeded spread. Each
% stroke is launched early by its predicted hand travel time (acceleration
% plus coast over the rest clearance), so the contacts land near the
% nominal beats instead of drifting with velocity. A nonpositive
% strike_gap_m asks for the clearance to be sized to the stroke: contact
% arrives just after the hand reaches full speed, which keeps the tip's
% whip phase, and with it the contact speed, consistent across strokes.
tstr = P.pre_delay_ms * 1e-3 + (0:K - 1) / max(1e-3, P.strike_rate_hz);
vstr = zeros(1, K);
state = P.seed * 2654435761;
for i = 1:K
  state = mod(state * 1664525 + 1013904223, 4294967296);
  x = state / 4294967296;
  vstr(i) = max(0.05, P.strike_velocity * (1 + P.velocity_spread * (2 * x - 1)));
end
gap = P.strike_gap_m;
if gap <= 0
  gap = 1.05 * max(vstr)^2 / (2 * P.strike_accel) + 0.005;
end
travel = gap ./ vstr + vstr / (2 * P.strike_accel);
tlaunch = tstr - travel;
off = -min(0, min(tlaunch));
tlaunch = tlaunch + off;

duration = off + (K - 1) / max(1e-3, P.strike_rate_hz) + P.tail_s;
nsteps = floor(duration * P.sample_rate * P.oversample);

KM = P.head_stiffness;
alM = P.head_alpha;
tipi = MB.tipi;
ps = MB.ps;

% Hand state machine: the hand accelerates downward, coasts at the stroke
% speed until the head makes contact (or a miss depth), holds briefly,
% then retracts at half speed and brakes so it arrives back at rest,
% ready for the next scheduled stroke. This replaces the paper's
% constant-acceleration hand, which never reverses. The machine runs one
% step ahead because the beam coupling needs u_H at k + 1.
aH = P.strike_accel;
aHr = P.retract_accel;
holdSteps = max(1, round(P.strike_hold_ms * 1e-3 / ht));
missDepth = gap + 0.06;
uHm = 0; uH0 = 0; uHp = 0;
vH = 0;
phase = 0;
strokeIdx = 1;
holdLeft = 0;
vtar = 0;

D1p = G.D1; D2p = G.D2;
T1t = G.T1t; T2 = G.T2;
F = G.F; Gi = G.Gi;
Rd = G.Rd; R2d = G.R2d;
AiB = S.AiB; AiC = S.AiC; Ai = S.Ai; Fi = S.Fi;
sNL = ht^2 / rhoh;

U = zeros(nrt, nth);
U1 = zeros(nrt, nth);
Uh3 = zeros(nrt, 1, nth);
Uh13 = zeros(nrt, 1, nth);
d = zeros(MB.nfree, 1);
d1 = zeros(MB.nfree, 1);

acc = zeros(1, nsteps);
fhist = zeros(1, nsteps);
shellActive = 0;
peakDisp = 0;
peakXi = 0;
inContact = 0;
cvel = zeros(1, K);

for k = 1:nsteps
  % ---- mallet: contact force from the state at step k ----
  xi = U(ls, js) - d(tipi) - gap;
  if xi > 0
    fM1 = KM * xi^alM;
    shellActive = 1;
    if inContact == 0 && strokeIdx <= K && cvel(strokeIdx) == 0
      cvel(strokeIdx) = (d(tipi) - d1(tipi)) / ht;   % tip speed at contact
    end
    inContact = 1;
    if xi > peakXi
      peakXi = xi;
    end
  else
    fM1 = 0;
    inContact = 0;
  end
  fhist(k) = fM1;
  hA = -(uHp - 2 * uH0 + uHm);
  hB = -(ht / 2) * (uHp - uHm);
  hC = -ht^2 * (ps * uH0 + ((1 - ps) / 2) * (uHp + uHm));
  dnew = MB.ASiB * d + MB.ASiC * d1 + MB.vM * hA + MB.vR * hB + MB.vK * hC ...
    + (ht^2 * fM1) * MB.tipcol;

  % ---- advance the hand one step ----
  if phase == 0
    if strokeIdx <= K && k * ht >= tlaunch(strokeIdx)
      phase = 1;
      vtar = vstr(strokeIdx);
    end
  elseif phase == 1
    vH = max(vH - aH * ht, -vtar);
    if inContact ~= 0
      phase = 3;
      holdLeft = holdSteps;
    elseif vH <= -vtar
      phase = 2;
    end
  elseif phase == 2
    if inContact ~= 0
      phase = 3;
      holdLeft = holdSteps;
    elseif uHp <= -missDepth
      phase = 4;
    end
  elseif phase == 3
    holdLeft = holdLeft - 1;
    if holdLeft <= 0
      phase = 4;
    end
  elseif phase == 4
    vH = min(vH + aHr * ht, 0.5 * vtar);
    if vH > 0 && uHp >= -vH^2 / (2 * aHr)
      phase = 5;
    end
  else
    vH = vH - aHr * ht;
    if vH <= 0 || uHp >= 0
      vH = 0;
      phase = 0;
      strokeIdx = strokeIdx + 1;
    end
  end
  uHm = uH0;
  uH0 = uHp;
  uHp = uHp + vH * ht;

  if shellActive ~= 0
    % ---- shell derivatives and von Karman brackets at step k ----
    Ur = D1p * U;
    Urr = D2p * U;
    v1 = Rd .* Ur + R2d .* (U * T2);
    v2 = Rd .* (Ur * T1t) - R2d .* (U * T1t);
    LapU = Urr + v1;
    rhsA = EhRs * LapU - Eh2 * 2 * (Urr .* v1 - v2 .* v2);
    rhsA(1:3, :) = 0;
    rhsA(Nr:nrt, :) = 0;
    Ph3 = pagemtimes(Fi, reshape(rhsA * F, nrt, 1, nth));
    Phi = reshape(Ph3, nrt, nth) * Gi;
    Pr = D1p * Phi;
    Prr = D2p * Phi;
    w1 = Rd .* Pr + R2d .* (Phi * T2);
    w2 = Rd .* (Pr * T1t) - R2d .* (Phi * T1t);
    NL = Prr .* v1 + Urr .* w1 - 2 * w2 .* v2 - (Prr + w1) / Rs;
    NL(ls, js) = NL(ls, js) - fM1 * fsc;
    NL(1:nlr0 - 1, :) = 0;
    NL(Nr + 1:nrt, :) = 0;
    NLh = pagemtimes(Ai, reshape((sNL * NL) * F, nrt, 1, nth));
    Uh3n = pagemtimes(AiB, Uh3) + pagemtimes(AiC, Uh13) + NLh;
    Unew = reshape(Uh3n, nrt, nth) * Gi;

    acc(k) = (Unew(lp, jp) - 2 * U(lp, jp) + U1(lp, jp)) / ht^2;
    U1 = U;
    U = Unew;
    Uh13 = Uh3;
    Uh3 = Uh3n;
    if mod(k, 512) == 0
      pd = max(abs(U(:, js)));
      if pd > peakDisp
        peakDisp = pd;
      end
      if ~isfinite(pd)
        error('cymbal2: simulation diverged at t = %.4f s (reduce strike_velocity, or raise n_theta / oversample)', k * ht);
      end
    end
  end

  d1 = d;
  d = dnew;
end

% Contact bookkeeping from the force history
inC = fhist > 0;
dinfo.peakForceN = max(fhist);
dinfo.contactSteps = sum(inC);
dinfo.contactMs = 1000 * ht * dinfo.contactSteps;
dinfo.contactEvents = sum(diff([0, inC]) > 0);
dinfo.peakXiMm = peakXi * 1000;
dinfo.peakDispM = peakDisp;
dinfo.peakDispThickness = peakDisp / P.thickness_m;
dinfo.peakAccel = max(abs(acc));
dinfo.nsteps = nsteps;
dinfo.simRate = P.sample_rate * P.oversample;
dinfo.tstr = tstr;
dinfo.vstr = vstr;
dinfo.meanV = mean(vstr);
dinfo.contactVel = cvel;   % tip speed at first contact of each stroke
dinfo.strokesPlayed = min(strokeIdx - 1 + (phase > 0), K);
end