// Why the stiff terms are integrated implicitly. The test problem is one
// undamped mode, u'' = -w^2 u, stepped by the paper's scheme
//   (1 + g z) u+ = (2 - (1 - 2g) z) u - (1 + g z) u-,   z = (w ht)^2,
// which is the explicit central difference at g = 0. The slider pushes the
// mode frequency past the explicit ceiling w ht = 2 while the implicit
// scheme stays on the unit circle.

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

const HEIGHT = 250;
const G = 0.2498;
const NSTEPS = 90;

function stepScheme(z: number, g: number): { u: number[]; rho: number } {
  const A = 1 + g * z;
  const B = 2 - (1 - 2 * g) * z;
  // exact start: u(-dt) = u(0) cos(w dt) with w dt = sqrt(z)
  const u: number[] = [1, Math.cos(Math.sqrt(z))];
  // recurrence runs forward from (u[1], u[0])
  let um = u[0];
  let uc = u[1];
  for (let k = 2; k < NSTEPS; k++) {
    const up = (B * uc - A * um) / A;
    u.push(up);
    um = uc;
    uc = up;
  }
  // spectral radius of the two-step recurrence
  const disc = B * B - 4 * A * A;
  let rho: number;
  if (disc <= 0) rho = 1;
  else rho = (Math.abs(B) + Math.sqrt(disc)) / (2 * A);
  return { u, rho };
}

// For the audio comparison: the frequency the implicit scheme actually
// produces for this mode, from the phase angle of its recurrence roots.
function schemeFreq(z: number, g: number, stepRate: number): number {
  const A = 1 + g * z;
  const B = 2 - (1 - 2 * g) * z;
  const theta = Math.atan2(Math.sqrt(Math.max(0, 4 * A * A - B * B)), B);
  return (theta / (2 * Math.PI)) * stepRate;
}

const AUDIO_SR = 24000;
// A pretend step rate low enough that the slider's whole range is audible.
const DEMO_STEP_RATE = 8000;

function renderDetuneAudio(fTrue: number, fNum: number): Float32Array<ArrayBuffer> {
  const seg = 0.4;
  const gap = 0.1;
  const n = Math.round(AUDIO_SR * (2 * seg + gap));
  const y = new Float32Array(n);
  const tone = (f: number, start: number) => {
    const n0 = Math.round(start * AUDIO_SR);
    const n1 = n0 + Math.round(seg * AUDIO_SR);
    const fade = Math.round(0.012 * AUDIO_SR);
    for (let k = n0; k < n1; k++) {
      let g2 = 0.35;
      if (k - n0 < fade) g2 *= (k - n0) / fade;
      if (n1 - k < fade) g2 *= (n1 - k) / fade;
      y[k] = g2 * Math.sin((2 * Math.PI * f * (k - n0)) / AUDIO_SR);
    }
  };
  tone(fTrue, 0);
  tone(fNum, seg + gap);
  return y;
}

export default function StabilityFig() {
  const [wht, setWht] = useState(2.2);
  const play = usePlayer();

  const z = wht * wht;
  const ex = stepScheme(z, 0);
  const im = stepScheme(z, G);
  const fTrue = (wht / (2 * Math.PI)) * DEMO_STEP_RATE;
  const fNum = schemeFreq(z, G, DEMO_STEP_RATE);

  const ref = useCanvas(
    (ctx, w, h) => {
      drawFrame(ctx, w, h);
      const panels: { u: number[]; rho: number; name: string; color: string }[] = [
        { ...ex, name: 'explicit (g = 0)', color: PLOT_COLORS.warn },
        { ...im, name: `implicit (g = ${G})`, color: PLOT_COLORS.trace },
      ];
      const pw = (w - 40) / 2;
      panels.forEach((p, pi) => {
        const x0 = 14 + pi * (pw + 12);
        const y0 = 12;
        const ph = h - 50;
        ctx.strokeStyle = PLOT_COLORS.grid;
        ctx.strokeRect(x0, y0, pw, ph);
        ctx.beginPath();
        ctx.moveTo(x0, y0 + ph / 2);
        ctx.lineTo(x0 + pw, y0 + ph / 2);
        ctx.stroke();

        // clip to the panel so divergence does not spray over the figure
        ctx.save();
        ctx.beginPath();
        ctx.rect(x0, y0, pw, ph);
        ctx.clip();
        ctx.strokeStyle = p.color;
        ctx.lineWidth = 1.6;
        ctx.beginPath();
        for (let k = 0; k < NSTEPS; k++) {
          const x = x0 + (k / (NSTEPS - 1)) * pw;
          const y = y0 + ph / 2 - Math.max(-1e6, Math.min(1e6, p.u[k])) * (ph / 2.4);
          if (k === 0) ctx.moveTo(x, y);
          else ctx.lineTo(x, y);
        }
        ctx.stroke();
        ctx.lineWidth = 1;
        ctx.restore();

        ctx.fillStyle = p.color;
        ctx.font = '11px system-ui';
        ctx.textAlign = 'left';
        const growth =
          p.rho <= 1 + 1e-12 ? 'bounded (|amplification| = 1)' : `grows ×${p.rho.toFixed(3)} per step`;
        ctx.fillText(`${p.name}: ${growth}`, x0 + 6, y0 + ph + 16);
      });

      ctx.fillStyle = PLOT_COLORS.axis;
      ctx.textAlign = 'center';
      ctx.fillText(
        `one mode at ω·ht = ${wht.toFixed(2)}  (explicit ceiling: 2)`,
        w / 2,
        h - 8,
      );
    },
    [wht, z],
    HEIGHT,
  );

  return (
    <Figure
      caption={
        <>
          Ninety steps of one undamped mode under both schemes, at the frequency set by the
          slider. The explicit central difference is exact in cost and cheap in memory, but
          the moment &omega;&middot;h<sub>t</sub> passes 2 its recurrence has a root outside
          the unit circle and the mode explodes; a plate's biharmonic stiffness puts its
          highest grid modes far past that ceiling at any audio time step. With the implicit
          weighting at g &ge; 1/4 the roots stay on the unit circle at every frequency (the
          paper's g = 0.2498 sits just below 1/4, and is stable to &omega;&middot;h
          <sub>t</sub> &asymp; 71, far beyond anything the grid carries). The price of the
          implicit form is a linear solve per step, which Section 7 makes affordable, and a
          detuning that grows with frequency: the play button compares the mode's true pitch
          against the pitch the scheme actually produces, scaled into the audible range by
          pretending the step rate were 8 kHz. Near the left end of the slider the two tones
          are indistinguishable, which is the regime the audible modes live in; the real
          simulation steps at 96 kHz, so a mode this badly flattened would sit far above the
          spectrum that matters.
        </>
      }
    >
      <canvas ref={ref} />
      <FigControls>
        <FigSlider
          label="ω·ht"
          value={wht}
          min={0.3}
          max={4.5}
          step={0.01}
          format={(v) => v.toFixed(2)}
          onChange={setWht}
        />
        <button onClick={() => play(renderDetuneAudio(fTrue, fNum), AUDIO_SR)}>
          ▶ Hear the detuning: {fTrue.toFixed(0)} Hz true, then {fNum.toFixed(0)} Hz as the
          scheme plays it
        </button>
      </FigControls>
    </Figure>
  );
}