// The modal synthesis: a single mode struck by an impulse gives the closed
// form damped sinusoid; the radiation gain of a uniform piston equals the
// analytic value; the attack and the tail join; the sliced parts sum to the
// whole.
import { test } from "node:test";
import assert from "node:assert/strict";
import { makeGrid, defaultParams, type SimParams, listenerTerms, RHO_AIR, SPEED_OF_SOUND, nodeWeight } from "../src/params.ts";
import { uniformPlate } from "../src/model/reference.ts";
import { assemble } from "../src/model/assemble.ts";
import { computeModes, localRunner } from "../src/modal/modes.ts";
import { constantLoss, LN1000 } from "../src/loss.ts";
import {
  modalModel, modalRates, strikeGains, simulateContact, listenerWeights, radiationGains, tailStart, nodalAccelerationTable, attackPressure, modalTail, movieFrames, synthesize, sliceModel, sliceContact, type ModalModel,
} from "../src/modal/synth.ts";

const steel = { E: 200e9, nu: 0.3, rho: 7850 };

async function smallModel(): Promise<{ model: ModalModel; pl: ReturnType<typeof uniformPlate> }> {
  const g = makeGrid(0.2, 0.15, 0.0125);
  const pl = uniformPlate(g, steel, 1e-3);
  pl.loss = constantLoss(1.0);
  const sys = assemble(pl);
  const modes = await computeModes(sys, 4000, localRunner(sys), { sliceSize: 20 });
  return { model: modalModel(pl, sys, modes), pl };
}

test("a single mode under an impulse follows the closed-form damped sinusoid", async () => {
  const { model, pl } = await smallModel();
  const p: SimParams = { ...defaultParams, fs: 48000, duration: 0.1, fps: 100 };
  // Keep only mode 0 by slicing the model.
  const one = sliceModel(model, 0, 1);
  const sigma = modalRates(one, pl.loss);
  const g = strikeGains(one, pl, p);
  const c = simulateContact(one, sigma, g, p, p.fps);
  assert.ok(c.peakForce > 0 && c.contactTime > 0);
  // After tEnd the mode is q(t) = Re[A exp(s (t - tEnd))]; check the movie frames against it.
  const nFrames = 10;
  const frames = movieFrames(one, sigma, c, p.fps, nFrames);
  const w = 2 * Math.PI * one.freqs[0], s = sigma[0], wd = Math.sqrt(w * w - s * s);
  for (let f = 1; f < nFrames; f++) {
    const t = f / p.fps - c.tEnd;
    if (t < 0) continue;
    const q = Math.exp(-s * t) * (c.aRe[0] * Math.cos(wd * t) - c.aIm[0] * Math.sin(wd * t));
    // Compare at a node away from the edge.
    const k = Math.floor(pl.grid.ny / 2) * pl.grid.nx + Math.floor(pl.grid.nx / 3);
    const phi = one.wShapes[one.wRow[k] * one.count];
    assert.ok(Math.abs(frames[f * one.nGrid + k] - phi * q) < 1e-9 + 1e-6 * Math.abs(phi * q));
  }
  assert.ok(Math.abs(sigma[0] - LN1000 / 1.0) < 1e-9, "decay rate from T60 = 1 s");
});

test("the radiation gain of a uniform piston equals the analytic Rayleigh sum", async () => {
  const { model, pl } = await smallModel();
  const p: SimParams = { ...defaultParams, listenerDist: 1.0, listenerAngle: 0 };
  // A fake model with one mode of uniform shape (mass-normalized to 1/sqrt(M)).
  const uniform: ModalModel = { ...sliceModel(model, 0, 1), freqs: new Float64Array([1000]), padSigma: new Float64Array([0]) };
  let M = 0;
  for (let k = 0; k < uniform.nGrid; k++) M += uniform.masses[k];
  const phi = 1 / Math.sqrt(M);
  for (let r = 0; r < uniform.nRows; r++) uniform.wShapes[r] = phi;
  const { weight, delay } = listenerWeights(pl, p);
  const sigma = new Float64Array([0]);
  const gains = radiationGains(uniform, sigma, weight, delay);
  // Analytic: G = sum_k rho dA / (2 pi r_k) phi exp(-i omega tau_k).
  const lt = listenerTerms(pl, p);
  let dmin = Infinity;
  for (let k = 0; k < lt.weight.length; k++) if (lt.weight[k] > 0) dmin = Math.min(dmin, lt.delay[k]);
  let re = 0, im = 0;
  const w = 2 * Math.PI * 1000;
  for (let k = 0; k < lt.weight.length; k++) {
    if (lt.weight[k] === 0) continue;
    const tau = lt.delay[k] - dmin;
    re += lt.weight[k] * phi * Math.cos(w * tau);
    im -= lt.weight[k] * phi * Math.sin(w * tau);
  }
  // The shapes are float32, so agreement is to about 1e-7.
  assert.ok(Math.abs(gains.re[0] - re) < 1e-6 * Math.abs(re) && Math.abs(gains.im[0] - im) < 1e-6 * Math.max(Math.abs(im), 1e-6));
  // On axis at 1 m the weights sum to rho A / (2 pi r) within the delay spread.
  let area = 0;
  for (let j = 0; j < pl.grid.ny; j++) for (let i = 0; i < pl.grid.nx; i++) area += nodeWeight(pl.grid, i, j) * pl.grid.hx * pl.grid.hy;
  const approx = (RHO_AIR * area) / (2 * Math.PI * 1.0) * phi;
  assert.ok(Math.abs(Math.hypot(re, im) / approx - 1) < 0.05);
  void SPEED_OF_SOUND;
});

test("the attack and the tail join, and the sliced parts sum to the whole", async () => {
  const { model, pl } = await smallModel();
  const p: SimParams = { ...defaultParams, fs: 48000, duration: 0.2, fps: 120, strikeVelocity: 2 };
  const whole = synthesize(model, pl, pl.loss, p, true);
  const nSamples = Math.round(p.duration * p.fs);
  assert.equal(whole.audio.length, nSamples);
  // No discontinuity at the tail start beyond the local sample-to-sample variation.
  const { maxDelay, weight, delay } = listenerWeights(pl, p);
  const n0 = Math.ceil(tailStart(whole.contact, maxDelay) * p.fs - 1e-9);
  let maxStep = 0, jump = Math.abs(whole.audio[n0] - whole.audio[n0 - 1]);
  for (let n = n0 + 2; n < n0 + 200; n++) maxStep = Math.max(maxStep, Math.abs(whole.audio[n] - whole.audio[n - 1]));
  assert.ok(jump < 3 * maxStep, `jump ${jump} vs typical step ${maxStep}`);
  // The tail decays: the last 20 ms is quieter than the first 20 ms.
  const rms = (a: Float32Array) => Math.sqrt(a.reduce((s, v) => s + v * v, 0) / a.length);
  assert.ok(rms(whole.audio.subarray(nSamples - 960)) < rms(whole.audio.subarray(0, 960)));
  // Split the modes into three slices and sum the parts.
  const sigma = modalRates(model, pl.loss);
  const g = strikeGains(model, pl, p);
  const c = simulateContact(model, sigma, g, p, p.fps);
  const tStart = tailStart(c, maxDelay);
  const len = Math.min(c.tableLen, Math.ceil(tStart * c.tableRate) + 2);
  const audio = new Float32Array(nSamples);
  const table = new Float32Array(model.nRows * len);
  const nFrames = Math.round(p.duration * p.fps);
  const frames = new Float32Array(nFrames * model.nGrid);
  const cuts = [0, Math.floor(model.count / 3), Math.floor((2 * model.count) / 3), model.count];
  for (let s = 0; s < 3; s++) {
    const sub = sliceModel(model, cuts[s], cuts[s + 1]);
    const cc = sliceContact(c, cuts[s], cuts[s + 1]);
    const sg = sigma.slice(cuts[s], cuts[s + 1]);
    nodalAccelerationTable(sub, cc, len, table);
    modalTail(sub, sg, cc, radiationGains(sub, sg, weight, delay), p.fs, tStart, nSamples, audio);
    movieFrames(sub, sg, cc, p.fps, nFrames, frames);
  }
  attackPressure(model, table, len, c.tableRate, weight, delay, p.fs, tStart, audio);
  let err = 0, nrm = 0;
  for (let n = 0; n < nSamples; n++) { err += (audio[n] - whole.audio[n]) ** 2; nrm += whole.audio[n] ** 2; }
  assert.ok(Math.sqrt(err / nrm) < 1e-5, `audio split error ${Math.sqrt(err / nrm)}`);
  err = 0; nrm = 0;
  for (let i = 0; i < frames.length; i++) { err += (frames[i] - whole.frames[i]) ** 2; nrm += whole.frames[i] ** 2; }
  assert.ok(Math.sqrt(err / nrm) < 1e-5, `frame split error ${Math.sqrt(err / nrm)}`);
});

test("a lighter hit gives a proportionally smaller signal in the linear model", async () => {
  const { model, pl } = await smallModel();
  const p1: SimParams = { ...defaultParams, fs: 48000, duration: 0.05, strikeVelocity: 0.2 };
  const p2: SimParams = { ...p1, strikeVelocity: 0.1 };
  const a1 = synthesize(model, pl, pl.loss, p1, false).audio, a2 = synthesize(model, pl, pl.loss, p2, false).audio;
  const peak = (a: Float32Array) => a.reduce((m, v) => Math.max(m, Math.abs(v)), 0);
  // The Hertz contact is nonlinear (a faster hit is shorter and brighter), so the ratio exceeds 2 somewhat.
  const ratio = peak(a1) / peak(a2);
  assert.ok(ratio > 1.8 && ratio < 3.5, `ratio ${ratio}`);
});