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.
%
% FLUCTUATION. What the equations above track is the MEAN energy of each mode,
% and the mean is not what a listener hears. A mode sitting in a turbulent
% field does not hold a smooth amplitude; it is continually being filled and
% drained by its neighbours, and its amplitude wanders about that mean on the
% timescale of the exchange. This matters more than it sounds: a few hundred
% modes with smooth envelopes are a few hundred audible partials, spaced tens
% of hertz apart with bandwidths of a couple of hertz, which is a chime. Give
% each one a wandering amplitude with a correlation time of a few milliseconds
% and each partial broadens by a couple of hundred hertz, the modal overlap
% passes one, and the comb closes into the continuum a wash actually is.
%
% So each mode carries a log-amplitude that follows an Ornstein-Uhlenbeck
% process discretized at the control rate,
%   rho_m = exp(-(energy exchanged this step) / E_m)
%   g_m  <- rho_m g_m + sqrt(1 - rho_m^2) turbulence z_m
% with the correlation rate taken as the mode's own fractional energy turnover
% rather than as a fitted time: a mode being fed while nearly empty
% decorrelates completely every step, a mode merely ringing down does not
% decorrelate at all. A running variance keeps exp(g) unbiased, so the
% fluctuation redistributes energy in time without creating or destroying any,
% and the process is seeded, so a render is reproducible. Where the cascade is
% off the turnover is zero, rho is one, g stays where it started, and the
% envelopes are exactly the closed-form exponentials of a linear plate.
% turbulence = 0 removes the fluctuation and leaves the mean-energy solve.
%
% 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). The projection is taken through the
% hammering mixing of cymbal_basis.m, since the modes the stick actually
% drives are the hybrids, not the ideal cos(n theta) shapes. 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
  proj = (basis.mix * reshape(proj, N, 1)).'; % hammered modes are hybrids
  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;

% one independent LCG stream per mode, for the amplitude fluctuation
sig = P.turbulence;
st1 = mod(floor(abs(P.seed)) * 2654435761 + (1:N).' * 40503, 4294967296);
st2 = mod(st1 * 69069 + 12345, 4294967296);
gl = zeros(N, 1);
gv = zeros(N, 1);

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) .* exp(gl - gv / 2);

  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;

  if sig > 0
    % fractional energy turnover this step sets the decorrelation
    moved = up + dn + arrU(1:N) + arrD;
    rho = exp(-min(30, moved ./ max(Ee, 1e-300)));
    [u1, st1] = cymbal_rngv(st1);
    [u2, st2] = cymbal_rngv(st2);
    zg = sqrt(-2 * log(u1)) .* cos(2 * pi * u2);
    w = 1 - rho.^2;
    gl = rho .* gl + sqrt(w) * sig .* zg;
    gv = rho.^2 .* gv + w * sig^2;
  end

  Ee = En .* decay;
end
env(:, nf + 1) = sqrt(2 * Ee) .* exp(gl - gv / 2);

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