// The crash in miniature: a Fermi-Pasta-Ulam chain, masses on springs whose
// force carries a quadratic (geometric) correction, started in its lowest
// mode. At small amplitude the energy stays where it was put; at large
// amplitude the nonlinearity couples the modes and energy migrates up the
// ladder. The shell's bending-stretching coupling does to its nodal-diameter
// modes what this chain does to its sine modes.

import { useMemo, useState } from 'react';
import { useCanvas, drawFrame, PLOT_COLORS, colormap } from '../../components/plots/canvasUtils';
import { Figure, FigControls, FigSlider, usePlayer } from '../figparts';

const HEIGHT = 260;
const N = 24;          // interior masses
const ALPHA = 0.25;    // quadratic coupling
const TSTEPS = 30000;
const DT = 0.06;
const FRAMES = 240;    // energy snapshots across the run

function simulate(amp: number) {
  const u = new Float64Array(N + 2);
  const v = new Float64Array(N + 2);
  for (let i = 1; i <= N; i++) u[i] = amp * Math.sin((Math.PI * i) / (N + 1));

  const wk = new Float64Array(N);
  for (let k = 1; k <= N; k++) wk[k - 1] = 2 * Math.sin((Math.PI * k) / (2 * (N + 1)));

  const energies = new Float32Array(FRAMES * N);
  const a = new Float64Array(N + 2);
  const accel = (uu: Float64Array) => {
    for (let i = 1; i <= N; i++) {
      const dr = uu[i + 1] - uu[i];
      const dl = uu[i] - uu[i - 1];
      a[i] = dr - dl + ALPHA * (dr * dr - dl * dl);
    }
  };

  const snap = (frame: number) => {
    for (let k = 1; k <= N; k++) {
      let q = 0, p = 0;
      for (let i = 1; i <= N; i++) {
        const s = Math.sin((Math.PI * k * i) / (N + 1));
        q += u[i] * s;
        p += v[i] * s;
      }
      const c = Math.sqrt(2 / (N + 1));
      q *= c;
      p *= c;
      energies[frame * N + (k - 1)] = 0.5 * (p * p + wk[k - 1] * wk[k - 1] * q * q);
    }
  };

  const perFrame = Math.floor(TSTEPS / FRAMES);
  // leapfrog: cheap and symplectic, so the ladder is not an artifact of decay
  accel(u);
  for (let f = 0; f < FRAMES; f++) {
    snap(f);
    for (let s = 0; s < perFrame; s++) {
      for (let i = 1; i <= N; i++) v[i] += 0.5 * DT * a[i];
      for (let i = 1; i <= N; i++) u[i] += DT * v[i];
      accel(u);
      for (let i = 1; i <= N; i++) v[i] += 0.5 * DT * a[i];
    }
  }
  return energies;
}

// Sonification: the same chain, retuned so mode 1 sits at 220 Hz (its top
// mode then lands near 3.5 kHz), lightly damped, and integrated at audio
// rate; the output is the acceleration of one off-center mass. At small
// amplitude it is a plain 220 Hz tone; at large amplitude the cascade is
// audible as a brightening after the onset.
const AUDIO_SR = 24000;
const AUDIO_F1 = 220;

export function renderChainAudio(amp: number): Float32Array<ArrayBuffer> {
  const dur = 2.5;
  const n = Math.round(AUDIO_SR * dur);
  const u = new Float64Array(N + 2);
  const v = new Float64Array(N + 2);
  const a = new Float64Array(N + 2);
  for (let i = 1; i <= N; i++) u[i] = amp * Math.sin((Math.PI * i) / (N + 1));
  const w1 = 2 * Math.sin(Math.PI / (2 * (N + 1)));
  const tauPerSec = (2 * Math.PI * AUDIO_F1) / w1; // natural time per real second
  const sub = 2;
  const dtau = tauPerSec / AUDIO_SR / sub;
  const cdamp = (2 * 2.76) / tauPerSec; // amplitude T60 near 2.5 s

  const accel = () => {
    for (let i = 1; i <= N; i++) {
      const dr = u[i + 1] - u[i];
      const dl = u[i] - u[i - 1];
      a[i] = dr - dl + ALPHA * (dr * dr - dl * dl) - cdamp * v[i];
    }
  };
  const iOut = 5;
  const y = new Float32Array(n);
  accel();
  for (let k = 0; k < n; k++) {
    for (let s = 0; s < sub; s++) {
      for (let i = 1; i <= N; i++) v[i] += 0.5 * dtau * a[i];
      for (let i = 1; i <= N; i++) u[i] += dtau * v[i];
      accel();
      for (let i = 1; i <= N; i++) v[i] += 0.5 * dtau * a[i];
    }
    y[k] = a[iOut];
  }
  let pk = 0;
  for (let k = 0; k < n; k++) pk = Math.max(pk, Math.abs(y[k]));
  const fadeIn = Math.round(0.004 * AUDIO_SR);
  const fadeOut = Math.round(0.12 * AUDIO_SR);
  for (let k = 0; k < n; k++) {
    let g = 0.45 / (pk || 1);
    if (k < fadeIn) g *= k / fadeIn;
    if (k > n - fadeOut) g *= (n - k) / fadeOut;
    y[k] *= g;
  }
  return y;
}

export default function CascadeFig() {
  const [amp, setAmp] = useState(4);
  const energies = useMemo(() => simulate(amp), [amp]);
  const play = usePlayer();

  const ref = useCanvas(
    (ctx, w, h) => {
      drawFrame(ctx, w, h);
      const padL = 70;
      const padR = 14;
      const padT = 10;
      const padB = 28;
      const pw = w - padL - padR;
      const ph = h - padT - padB;

      let e0 = 0;
      for (let k = 0; k < N; k++) e0 += energies[k];
      const floorDb = -55;

      // paint at data resolution, then stretch (putImageData ignores the
      // dpr transform, so go through an offscreen canvas as the plots do)
      const off = document.createElement('canvas');
      off.width = FRAMES;
      off.height = N;
      const octx = off.getContext('2d')!;
      const img = octx.createImageData(FRAMES, N);
      for (let mode = 0; mode < N; mode++) {
        const row = N - 1 - mode; // mode 1 at the bottom
        for (let frame = 0; frame < FRAMES; frame++) {
          const e = energies[frame * N + mode] / (e0 || 1);
          const db = 10 * Math.log10(Math.max(1e-9, e));
          const t = Math.min(1, Math.max(0, (db - floorDb) / -floorDb));
          const [r, g, b] = colormap(t);
          const o = 4 * (row * FRAMES + frame);
          img.data[o] = r;
          img.data[o + 1] = g;
          img.data[o + 2] = b;
          img.data[o + 3] = 255;
        }
      }
      octx.putImageData(img, 0, 0);
      ctx.imageSmoothingEnabled = false;
      ctx.drawImage(off, padL, padT, pw, ph);

      ctx.fillStyle = PLOT_COLORS.axis;
      ctx.font = '10px system-ui';
      ctx.textAlign = 'right';
      for (const k of [1, 6, 12, 18, 24]) {
        const yy = padT + ph - ((k - 0.5) / N) * ph;
        ctx.fillText(`${k}`, padL - 6, yy + 3);
      }
      ctx.save();
      ctx.translate(14, padT + ph / 2);
      ctx.rotate(-Math.PI / 2);
      ctx.textAlign = 'center';
      ctx.fillText('mode number (energy, dB re start)', 0, 0);
      ctx.restore();
      ctx.textAlign = 'center';
      ctx.fillText('time →', padL + pw / 2, h - 8);
      ctx.textAlign = 'left';
      ctx.fillText('all energy starts in mode 1 (bottom row)', padL + 6, h - 8);
    },
    [energies],
    HEIGHT,
  );

  return (
    <Figure
      caption={
        <>
          Modal energy (color, log scale) against time for a chain of {N} masses whose spring
          forces carry a quadratic geometric term, started with all of its energy in mode 1.
          Below an amplitude near 1 the system is effectively linear and the energy stays in
          the bottom row for the whole run. Push the amplitude up and the modes cease to be
          independent: energy climbs the ladder, and at the top of the range the spectrum
          fills in. This, plus damping, is the crash: the strike puts energy in at the bottom,
          and the nonlinearity spends the next several hundred milliseconds carrying it
          upward, which is why a crash brightens after the attack rather than at it. The play
          button integrates the same chain retuned so its lowest mode sits at 220 Hz and
          lightly damped: quiet strikes give a plain tone, and somewhere past an amplitude of
          2 the brightening-after-onset becomes plainly audible.
        </>
      }
    >
      <canvas ref={ref} />
      <FigControls>
        <FigSlider
          label="strike amplitude"
          value={amp}
          min={0.2}
          max={8}
          step={0.1}
          format={(v) => v.toFixed(1)}
          onChange={setAmp}
        />
        <button onClick={() => play(renderChainAudio(amp), AUDIO_SR)}>
          ▶ Hear this chain (mode 1 tuned to 220 Hz)
        </button>
      </FigControls>
    </Figure>
  );
}