function [env, tframe, DV, kf, cinfo] = cymbal_cascade(P, basis, tstr, jstr, taustr, sstr, duration)
% CYMBAL_CASCADE Modal excitation and the nonlinear flow of energy up the
% frequency ladder.
%
% This is the stage that makes the instrument a cymbal rather than a gong or a
% bell. A cymbal struck hard is not a linear resonator: the transverse
% displacement is comparable to the plate thickness, so the von Karman
% stretching of the mid-surface couples the modes, and energy injected in the
% low modes by the stick migrates upward over tens of milliseconds. That is
% what the ear hears as a crash swelling into shimmer after its attack, and it
% is why the same cymbal struck softly just pings: the coupling strength goes
% as the square of the amplitude, so a quiet stroke stays linear.
%
% Integrating the von Karman equations is out of reach here (the cubic
% coupling is a dense N^3 tensor, and the sibling repositories render in
% seconds). What is integrated instead is the ENERGY of each mode. Each mode
% is given an up-target, the lowest retained mode at least cascade_ratio times
% its frequency, and a down-target the same ratio below, and energy performs a
% biased random walk between them:
%   dE_m/dt = -(1 + beta) Gam_m E_m + (arrivals from below and above)
%           - 2 zeta_m om_m E_m
%   Gam_m   = cascade_gain om_m (w_m / t_mean)^2 = cascade_gain E_m / (om_m M t^2)
%   w_m^2   = E_m / (om_m^2 M)
% with beta = cascade_back < 1 the downward-to-upward ratio, so the net flux is
% toward high frequency. Steps are taken in log frequency rather than in mode
% index on purpose: four-wave interactions move energy by a fraction of an
% octave at a time, not by one mode, and a ladder defined on the index would
% make the swell depend on how densely the band happens to be resolved. The
% rate is proportional to om_m because the higher a mode sits the faster its
% interactions run, and to the square of that mode's own
% amplitude-to-thickness ratio because that is the order of the geometric
% nonlinearity. Taking the amplitude locally rather than as a plate-wide rms
% is what keeps the cascade honest at the top of the band: the flux out of a
% mode goes as the SQUARE of its energy, so a nearly empty mode passes nothing
% on, the front of the cascade advances rather than the whole band draining at
% once, and the swell stops on its own as the energy thins out.
%
% The scheme conserves energy exactly when the damping is off and the top of
% the ladder reflects (cascade_leak = 0), which the validation checks; with
% cascade_leak = 1 the flux out of the top is lost, standing in for the
% truncated modes above max_mode_hz. cascade_gain is the one fitted number in
% the model: it is a dimensionless coupling constant, set so that a default
% crash takes about the 40 ms a real one takes to open up.
%
% What this buys is a cymbal that responds to dynamics the way the instrument
% does, at the cost of one control-rate loop rather than a full nonlinear
% integration. What it gives up is phase: energy moves between modes but the
% intermodulation products of a real cascade, and the split-pair beating of an
% out-of-round cymbal, are not reproduced. cascade_gain = 0 turns the whole
% mechanism off and leaves an exactly linear plate.
%
% Each contact adds ½ (dv_m)^2 to the energy of every mode, where
%   dv_m = J S(om_m tau) (phi_m(x_c) . dhat)
% with S(x) = cos(x/2) / (1 - (x/pi)^2) the spectrum of the half-sine contact
% pulse and dhat the stroke direction (strike_angle_deg away from the surface
% normal, in the meridian plane). Energies add incoherently, which is right
% for a plate already ringing and exact for the first stroke into a silent one.
%
% Note that only the MAGNITUDE of the excitation survives into the energy
% state. The signed jumps DV and the frame each contact lands in are returned
% alongside, because cymbal_radiate.m needs them for the oscillator phases:
% amplitude comes from the cascade, phase from the linear impulse response.
%
% Returns env (N x nf+1, modal VELOCITY amplitude at each frame boundary),
% tframe (1 x nf+1), DV (N x number of contacts), kf (the frame index of each
% contact) and cinfo diagnostics.

N = basis.n;
omega = basis.omega;
zeta = basis.zeta;
h = max(1e-4, P.cascade_step_ms * 1e-3);
nf = max(1, ceil(duration / h));
tframe = (0:nf) * h;

% ---- modal velocity jump of each contact ----
I = length(tstr);
DV = zeros(N, I);
ang = P.strike_angle_deg * pi / 180;
for i = 1:I
  [uh, wh] = cymbal_at(basis, sstr(i));
  proj = -cos(ang) * wh + sin(ang) * uh;      % 1 x N, along the stroke
  x = omega * taustr(i);
  den = 1 - (x / pi).^2;
  S = cos(x / 2) ./ den;
  sing = abs(den) < 1e-6;
  S(sing) = pi / 4;                            % limit at x -> pi
  DV(:, i) = jstr(i) * S .* reshape(proj, N, 1);
end
kf = min(nf, max(1, floor(tstr / h) + 1));     % frame each contact lands in

decay = exp(-2 * zeta .* omega * h);
gmax = 0.5 / (h * (1 + P.cascade_back));
Mtot = basis.mass;
teff = basis.thkMean;

% ---- the log-frequency ladder ----
% upIdx(m) is the lowest mode at least cascade_ratio times f(m); modes with no
% such partner point at N + 1, the sink standing in for the truncated band.
% dnIdx(m) is the highest mode at most f(m) / cascade_ratio, or m itself at the
% bottom, which reflects.
ratio = max(1.02, P.cascade_ratio);
f = basis.f;
upIdx = zeros(N, 1);
dnIdx = zeros(N, 1);
for m = 1:N
  ju = find(f > f(m) * ratio);
  if length(ju) == 0
    upIdx(m) = N + 1;
  else
    upIdx(m) = ju(1);
  end
  jd = find(f < f(m) / ratio);
  if length(jd) == 0
    dnIdx(m) = m;
  else
    dnIdx(m) = jd(length(jd));
  end
end

env = zeros(N, nf + 1);
Ee = zeros(N, 1);
peakStrength = 0;
peakWrms = 0;
lost = 0;
for f = 1:nf
  hit = find(kf == f);
  for q = 1:length(hit)
    Ee = Ee + 0.5 * DV(:, hit(q)).^2;
  end
  env(:, f) = sqrt(2 * Ee);

  wrms2 = sum(Ee ./ omega.^2) / Mtot;
  strength = wrms2 / teff^2;
  if strength > peakStrength
    peakStrength = strength;
    peakWrms = sqrt(wrms2);
  end
  gam = min(gmax, P.cascade_gain * Ee ./ (omega * Mtot * teff^2));
  up = gam .* Ee * h;
  dn = P.cascade_back * gam .* Ee * h;
  arrU = accumarray(upIdx, up, [N + 1, 1]);
  arrD = accumarray(dnIdx, dn, [N, 1]);
  leak = arrU(N + 1);
  En = Ee - up - dn + arrU(1:N) + arrD;
  En(N) = En(N) + (1 - P.cascade_leak) * leak;
  lost = lost + P.cascade_leak * leak;
  Ee = En .* decay;
end
env(:, nf + 1) = sqrt(2 * Ee);

cinfo.frames = nf;
cinfo.stepS = h;
cinfo.peakStrength = peakStrength;
cinfo.peakWrmsUm = peakWrms * 1e6;
cinfo.cascadeLostJ = lost;
cinfo.injectedJ = sum(sum(0.5 * DV.^2));
end