/**
 * Every object in closed form: exact modes of the discrete schemes, and the
 * mic trace synthesized from them against precomputed room kernels.
 *
 * This is the mathematics behind the fast render path (src/studio/room.ts
 * records the kernels it consumes; src/studio/fast.ts strings the two
 * together). The facts it rests on, verified numerically in
 * scripts/modal-poc.ts and test/checks.ts:
 *
 *  - Every object's discrete operator has exact sine eigenmodes. The pinned
 *    5-point Laplacian's are product sines; the second difference's are
 *    sines with lam_m = (4/h²) sin²(mπ/2(n−1)); and the fourth difference
 *    with antisymmetric ghosts is exactly the square of the pinned second
 *    difference (the README's own point about it), so it shares those sines
 *    with eigenvalue lam². The damping terms are built from the same
 *    operators. Each modal amplitude therefore obeys a scalar two-term
 *    recursion — the leapfrog restricted to one mode — whose solution is a
 *    closed-form geometric series. No object is ever timestepped.
 *  - The coupling into the air is linear and one-way, so every source the
 *    air sees is a sum of *fixed spatial patterns* times scalar series in
 *    those modal amplitudes: an object's rad route is one pattern per mode
 *    (prof .* rad(phi_m)), and a string's bridge route is a single pattern
 *    (bprof) times the slope its modes arrive at the far end with.
 *  - The whole system is LTI, so the mic trace is the sum over (pattern,
 *    series) pairs of the series convolved with that pattern's impulse
 *    response of the room — the kernel, which depends only on geometry.
 *
 *  Each series being exactly Re(c·z^k) per characteristic root, its
 *  convolution reduces to a complex one-pole recursion over the kernel: the
 *  whole synthesis is O(poles × samples) scalar work on the CPU, with the
 *  closed-form exponential tail falling out for free once the kernel ends.
 *
 * Everything here matches the compiled f32 schemes to their own rounding
 * noise; the approximations are the membrane mode set kept (the caller's
 * eps/fmax — line objects have few modes and always keep all) and the
 * kernel length.
 */
import type { AirGrid, ObjectGrid } from '../grids.ts';
import type { ObjectSpec, Params } from './types.ts';
import { suffixed } from './types.ts';

/** One root of a mode's recursion, carrying the displacement series:
 *  alpha[k] += Re(a · z^k). A real root is the zIm = 0 case. */
export interface BasePole {
  aRe: number;
  aIm: number;
  zRe: number;
  zIm: number;
}

/** One term of a source series, ready to filter a kernel: Re(c · z^k). */
export interface SeriesPole {
  cRe: number;
  cIm: number;
  zRe: number;
  zIm: number;
}

/** One mode of one object: which sine it is (n = 0 on a line grid), its
 *  discrete frequency and decay under the current parameters, and its
 *  displacement series as poles. */
export interface ModalMode {
  /** Sine indices, 1-based; a line mode has n = 0. */
  m: number;
  n: number;
  /** Eigenvalue of minus the discrete second difference / Laplacian. */
  lam: number;
  /** Discrete frequency (Hz) and 1/e decay time (s) under these params. */
  freq: number;
  tau: number;
  /** Peak source strength |acc coefficient|, for truncation decisions. */
  amp: number;
  poles: BasePole[];
}

/** sin(mπi/(n−1)) for m = 1..n−2, i = 0..n−1, row-major by m. */
export function sineTable(n: number): Float64Array {
  const t = new Float64Array((n - 2) * n);
  for (let m = 1; m <= n - 2; m++)
    for (let i = 0; i < n; i++) t[(m - 1) * n + i] = Math.sin((m * Math.PI * i) / (n - 1));
  return t;
}

/* ----------------------------------------------------- the recursions --- */

/**
 * Solve one mode's recursion in closed form.
 *
 * Restricted to an eigenvalue, every object's step is
 *   a[k+1] = cA·a[k] + cB·a[k−1],
 * and the excitation is the pair (a[0], a[−1]): a strike is a velocity
 * (a[0] = 0, a[−1] = b), a pluck a held displacement (a[0] = a[−1] = b).
 * The characteristic roots give a[k] as a geometric series; underdamped
 * modes (all the audible ones) as one conjugate pair written Re(a·z^k),
 * modes pushed past critical damping near the grid's Nyquist as two real
 * roots.
 */
export function solveBasePoles(
  cA: number,
  cB: number,
  a0: number,
  am1: number,
): BasePole[] {
  const disc = cA * cA + 4 * cB;
  if (disc < 0) {
    const zRe = cA / 2;
    const zIm = Math.sqrt(-disc) / 2;
    const zAbs2 = zRe * zRe + zIm * zIm;
    // alpha[0] = Re(a), alpha[-1] = Re(a·conj(z))/|z|².
    const aRe = a0;
    const aIm = zIm !== 0 ? (am1 * zAbs2 - aRe * zRe) / zIm : 0;
    return [{ aRe, aIm, zRe, zIm }];
  }
  const sq = Math.sqrt(Math.max(disc, 1e-300));
  const r1 = (cA + sq) / 2;
  const r2 = (cA - sq) / 2;
  if (r1 === r2 || r1 === 0 || r2 === 0) return [];
  // a1 + a2 = a0;  a1/r1 + a2/r2 = am1.
  const det = 1 / r1 - 1 / r2;
  const a1 = (am1 - a0 / r2) / det;
  const a2 = a0 - a1;
  return [
    { aRe: a1, aIm: 0, zRe: r1, zIm: 0 },
    { aRe: a2, aIm: 0, zRe: r2, zIm: 0 },
  ];
}

const freqOf = (poles: BasePole[], dt: number): { freq: number; tau: number } => {
  if (poles.length === 0) return { freq: 0, tau: 0 };
  const p = poles[0];
  if (p.zIm !== 0) {
    const zAbs = Math.hypot(p.zRe, p.zIm);
    return {
      freq: Math.abs(Math.atan2(p.zIm, p.zRe)) / (2 * Math.PI * dt),
      tau: zAbs < 1 ? -dt / Math.log(zAbs) : Infinity,
    };
  }
  const rMax = Math.max(...poles.map((q) => Math.abs(q.zRe)));
  return {
    freq: poles.some((q) => q.zRe < 0) ? 1 / (2 * dt) : 0,
    tau: rMax > 0 && rMax < 1 ? -dt / Math.log(rMax) : Infinity,
  };
};

/** Complex-multiply a base pole's coefficient by f(z). */
const scalePole = (p: BasePole, fRe: number, fIm: number): SeriesPole => ({
  cRe: p.aRe * fRe - p.aIm * fIm,
  cIm: p.aRe * fIm + p.aIm * fRe,
  zRe: p.zRe,
  zIm: p.zIm,
});

/** The acceleration series (un − 2u + um)/dt² of a mode, as poles:
 *  coefficient scaled by (z − 2 + 1/z)/dt² per root. */
export function accPoles(mode: ModalMode, dt: number, gain = 1): SeriesPole[] {
  return mode.poles.map((p) => {
    const zAbs2 = p.zRe * p.zRe + p.zIm * p.zIm;
    const fRe = (gain * (p.zRe - 2 + p.zRe / zAbs2)) / (dt * dt);
    const fIm = (gain * (p.zIm - p.zIm / zAbs2)) / (dt * dt);
    return scalePole(p, fRe, fIm);
  });
}

/** The series alpha[k+1] (the object's *new* displacement at step k, which
 *  is what a string's bridge slope is read from): coefficient scaled by z. */
export function nextPoles(mode: ModalMode, gain = 1): SeriesPole[] {
  return mode.poles.map((p) => scalePole(p, gain * p.zRe, gain * p.zIm));
}

/* ------------------------------------------------------ the mode sets --- */

export interface ModeSetOptions {
  /** Membranes only: keep modes whose source strength is at least eps × the
   *  strongest. 0 keeps everything the strike touches. Line objects have
   *  few modes and always keep all. */
  eps?: number;
  /** Membranes only: keep modes below this discrete frequency, Hz. */
  fmax?: number;
}

const ampOf = (poles: BasePole[], dt: number): number => {
  let amp = 0;
  for (const p of accPoles({ poles } as ModalMode, dt)) amp += Math.hypot(p.cRe, p.cIm);
  return amp;
};

/** Eigenvalue of minus the pinned second difference for sine mode m. */
export const lineLam = (g: ObjectGrid, m: number): number => {
  const s = Math.sin((m * Math.PI) / (2 * (g.n1 - 1)));
  return (4 / (g.h1 * g.h1)) * s * s;
};

/** Eigenvalue of minus the pinned 5-point Laplacian for sine mode (m, n). */
export const membraneLam = (g: ObjectGrid, m: number, n: number): number => {
  const s1 = Math.sin((m * Math.PI) / (2 * (g.n1 - 1)));
  const s2 = Math.sin((n * Math.PI) / (2 * (g.n2 - 1)));
  return (4 / (g.h1 * g.h1)) * s1 * s1 + (4 / (g.h2 * g.h2)) * s2 * s2;
};

/** Project a line field onto every interior sine mode, normalized. */
export function projectLine(g: ObjectGrid, v: Float64Array): Float64Array {
  const n = g.n1;
  const t = sineTable(n);
  const out = new Float64Array(n - 2);
  for (let m = 0; m < n - 2; m++) {
    let s = 0;
    for (let i = 0; i < n; i++) s += t[m * n + i] * v[i];
    out[m] = s / ((n - 1) / 2);
  }
  return out;
}

/** Project a membrane field onto every interior sine mode: two separable
 *  contractions, out[(m−1) + (n1−2)(n−1)], normalized. */
export function projectModes(g: ObjectGrid, v: Float64Array): Float64Array {
  const { n1, n2 } = g;
  const M1 = n1 - 2;
  const M2 = n2 - 2;
  const s1 = sineTable(n1);
  const s2 = sineTable(n2);
  const tmp = new Float64Array(M1 * n2);
  for (let m = 0; m < M1; m++) {
    for (let j = 0; j < n2; j++) {
      let s = 0;
      for (let i = 0; i < n1; i++) s += s1[m * n1 + i] * v[i + n1 * j];
      tmp[m + M1 * j] = s;
    }
  }
  const out = new Float64Array(M1 * M2);
  const norm = ((n1 - 1) / 2) * ((n2 - 1) / 2);
  for (let m = 0; m < M1; m++) {
    for (let n = 0; n < M2; n++) {
      let s = 0;
      for (let j = 0; j < n2; j++) s += s2[n * n2 + j] * tmp[m + M1 * j];
      out[m + M1 * n] = s / norm;
    }
  }
  return out;
}

/** One object's parameter, unsuffixed, defaulting hard on absence. */
const par = (params: Params, id: string, key: string): number => {
  const v = params[suffixed(key, id)];
  if (v === undefined) throw new Error(`fast path: parameter '${suffixed(key, id)}' missing`);
  return v;
};

/**
 * The mode set for one object under its current parameters, sorted
 * strongest first: everything the synthesis needs. For strings and bars
 * every interior mode is kept (there are a few dozen); for membranes
 * `opts` bounds the set, which with the kernel length is the fast path's
 * only approximation.
 *
 * The excitations are exactly the generated init's: a string's triangular
 * pluck held then released (a[0] = a[−1]), a bar's and membrane's
 * mallet-shaped initial velocity (a[0] = 0).
 */
export function objectModeSet(
  obj: ObjectSpec,
  g: ObjectGrid,
  dt: number,
  params: Params,
  opts: ModeSetOptions = {},
): ModalMode[] {
  const id = obj.id;
  const t60 = par(params, id, 't60');
  const sig1 = par(params, id, 'sig1');
  const d = (6.91 / t60) * dt;
  /** cA/cB from the stiffness's squared angular frequency at eigenvalue lam. */
  const coeffs = (w2: number, lam: number): { cA: number; cB: number } => ({
    cA: (2 - dt * dt * w2 - 2 * sig1 * dt * lam) / (1 + d),
    cB: (-(1 - d) + 2 * sig1 * dt * lam) / (1 + d),
  });

  if (obj.type === 'membrane') {
    const eps = opts.eps ?? 1e-6;
    const fmax = opts.fmax ?? Infinity;
    const f0 = par(params, id, 'f0');
    const cm = (2 * f0) / Math.sqrt(1 / (g.L1 * g.L1) + 1 / (g.L2 * g.L2));
    const um0 = new Float64Array(g.npts);
    const sx0 = par(params, id, 'strikex') * g.L1;
    const sy0 = par(params, id, 'strikey') * g.L2;
    const sr = par(params, id, 'striker');
    const vel = par(params, id, 'vel');
    for (let k = 0; k < g.npts; k++) {
      const dx = g.xs64[k] - sx0;
      const dy = g.ys64[k] - sy0;
      um0[k] = -(dt * vel) * g.pin[k] * Math.exp(-(dx * dx + dy * dy) / (sr * sr));
    }
    const b = projectModes(g, um0);
    const M1 = g.n1 - 2;
    const all: ModalMode[] = [];
    for (let n = 1; n <= g.n2 - 2; n++) {
      for (let m = 1; m <= M1; m++) {
        const bmn = b[m - 1 + M1 * (n - 1)];
        if (bmn === 0) continue;
        const lam = membraneLam(g, m, n);
        const { cA, cB } = coeffs(cm * cm * lam, lam);
        const poles = solveBasePoles(cA, cB, 0, bmn);
        const { freq, tau } = freqOf(poles, dt);
        if (poles.length === 0 || freq > fmax) continue;
        all.push({ m, n, lam, freq, tau, amp: ampOf(poles, dt), poles });
      }
    }
    all.sort((a, bb) => bb.amp - a.amp);
    const ampMax = all.length ? all[0].amp : 0;
    return all.filter((md) => md.amp >= eps * ampMax);
  }

  // A line object: pluck (string) or strike (bar), then the same recursion
  // with the type's own stiffness.
  const n = g.n1;
  const u0 = new Float64Array(n);
  let w2: (lam: number) => number;
  if (obj.type === 'string') {
    const f0 = par(params, id, 'f0');
    const B = par(params, id, 'B');
    const cs = 2 * g.L1 * f0;
    const kap2 = (B * (cs * g.L1) ** 2) / (Math.PI * Math.PI);
    w2 = (lam) => cs * cs * lam + kap2 * lam * lam;
    const pp = par(params, id, 'pluckpos') * g.L1;
    const amp = par(params, id, 'amp');
    for (let i = 0; i < n; i++) {
      const a = g.xs64[i] / pp;
      const bb = (g.L1 - g.xs64[i]) / (g.L1 - pp);
      u0[i] = amp * 0.5 * (a + bb - Math.abs(a - bb));
    }
  } else {
    const f1 = par(params, id, 'f1');
    const kap = (2 * f1 * g.L1 * g.L1) / Math.PI;
    w2 = (lam) => kap * kap * lam * lam;
    const sp = par(params, id, 'strikepos') * g.L1;
    const sr = par(params, id, 'striker');
    const vel = par(params, id, 'vel');
    for (let i = 0; i < n; i++) {
      const r = (g.xs64[i] - sp) / sr;
      u0[i] = -(dt * vel) * g.pin[i] * Math.exp(-r * r);
    }
  }
  const b = projectLine(g, u0);
  const out: ModalMode[] = [];
  for (let m = 1; m <= n - 2; m++) {
    if (b[m - 1] === 0) continue;
    const lam = lineLam(g, m);
    const { cA, cB } = coeffs(w2(lam), lam);
    // A pluck starts held (a[0] = a[−1]); a strike starts moving (a[0] = 0).
    const poles =
      obj.type === 'string'
        ? solveBasePoles(cA, cB, b[m - 1], b[m - 1])
        : solveBasePoles(cA, cB, 0, b[m - 1]);
    if (poles.length === 0) continue;
    const { freq, tau } = freqOf(poles, dt);
    out.push({ m, n: 0, lam, freq, tau, amp: ampOf(poles, dt), poles });
  }
  return out;
}

/* -------------------------------------------------------- air tables ---- */

/** The bilinearly interpolated sine tables the coupling ops effectively
 *  apply, per air column: rad of a product of sines factorizes per axis,
 *  and a line object's rad is the x factor alone. */
export function radSineTables(
  g: ObjectGrid,
  air: AirGrid,
): { sx: Float32Array; sy: Float32Array; M1: number; M2: number } {
  const { n1, n2, h1, h2, L1, L2 } = g;
  const M1 = n1 - 2;
  const s1 = sineTable(n1);
  const sx = new Float32Array(M1 * air.nx);
  for (let ix = 0; ix < air.nx; ix++) {
    const xw = -air.Lx / 2 + (ix + 0.5) * air.h;
    const u = (xw + L1 / 2) / h1;
    const j = Math.floor(u);
    if (j < 0 || j >= n1 - 1) continue;
    const f = u - j;
    for (let m = 0; m < M1; m++)
      sx[m * air.nx + ix] = (1 - f) * s1[m * n1 + j] + f * s1[m * n1 + j + 1];
  }
  if (g.kind === 'line') {
    // No second axis: the y factor is 1 everywhere (the tube profile in the
    // scene's prof carries the y and z structure).
    return { sx, sy: new Float32Array(air.ny).fill(1), M1, M2: 1 };
  }
  const M2 = n2 - 2;
  const s2 = sineTable(n2);
  const sy = new Float32Array(M2 * air.ny);
  for (let iy = 0; iy < air.ny; iy++) {
    const yw = -air.Ly / 2 + (iy + 0.5) * air.h;
    const v = (yw + L2 / 2) / h2;
    const j = Math.floor(v);
    if (j < 0 || j >= n2 - 1) continue;
    const f = v - j;
    for (let n = 0; n < M2; n++)
      sy[n * air.ny + iy] = (1 - f) * s2[n * n2 + j] + f * s2[n * n2 + j + 1];
  }
  return { sx, sy, M1, M2 };
}

/* --------------------------------------------------------- synthesis ---- */

/** One (series, kernel) pair: the poles' series convolved with the kernel. */
export interface Contribution {
  poles: SeriesPole[];
  /** One kernel row, `kernelSteps` samples. */
  kernel: Float32Array;
}

/**
 * Synthesize the mic trace: every contribution's kernel filtered by its
 * poles, summed. The trace may be longer than the kernels: a pole keeps
 * ringing the room's response down after the kernel ends, which is exactly
 * the closed-form tail.
 */
export function synthesize(
  contribs: Contribution[],
  kernelSteps: number,
  totalSteps: number,
): Float32Array {
  const y = new Float64Array(totalSteps);
  for (const c of contribs) {
    const h = c.kernel;
    for (const pl of c.poles) {
      const { cRe, cIm, zRe, zIm } = pl;
      let wRe = 0;
      let wIm = 0;
      let wPk = 0;
      const lim = Math.min(kernelSteps, totalSteps);
      for (let t = 0; t < lim; t++) {
        const nRe = zRe * wRe - zIm * wIm + h[t];
        const nIm = zRe * wIm + zIm * wRe;
        wRe = nRe;
        wIm = nIm;
        const m2 = wRe * wRe + wIm * wIm;
        if (m2 > wPk) wPk = m2;
        y[t] += cRe * wRe - cIm * wIm;
      }
      // The closed-form tail: the pole rings the kernel's response down with
      // no further input. Cut it 240 dB under its own peak — far below both
      // f32 and the .wav — which also keeps the recursion out of subnormal
      // floats, where every multiply is slow.
      const cut = Math.max(wPk * 1e-24, 1e-280);
      for (let t = lim; t < totalSteps; t++) {
        const nRe = zRe * wRe - zIm * wIm;
        const nIm = zRe * wIm + zIm * wRe;
        wRe = nRe;
        wIm = nIm;
        if (wRe * wRe + wIm * wIm < cut) break;
        y[t] += cRe * wRe - cIm * wIm;
      }
    }
  }
  const out = new Float32Array(totalSteps);
  for (let t = 0; t < totalSteps; t++) out[t] = y[t];
  return out;
}