// The viscoelastic felt: its chain of Maxwell arms has the loss factor it
// was built for across the audio band, the storage stiffness it was
// normalized to, and an exact time step.
import { test } from "node:test";
import assert from "node:assert/strict";
import { feltChain, feltModulus, feltStepConstants, FELT_ARMS } from "../src/felt.ts";

test("the loss factor of the chain is flat across the band and the storage is one at the reference", () => {
  for (const eta of [0.05, 0.3, 0.6]) {
    for (const fRef of [15, 96, 400]) {
      const ch = feltChain(eta, fRef);
      assert.equal(ch.tau.length, FELT_ARMS);
      const ref = feltModulus(ch, 2 * Math.PI * fRef);
      assert.ok(Math.abs(ref.storage - 1) < 1e-12);
      assert.ok(Math.abs(ref.loss / ref.storage / eta - 1) < 0.02, `eta ${eta} at the reference ${fRef} Hz: ${ref.loss / ref.storage}`);
      let worst = 0, prev = 0;
      for (let f = 5; f <= 10000; f *= 1.1) {
        const m = feltModulus(ch, 2 * Math.PI * f);
        worst = Math.max(worst, Math.abs(m.loss / m.storage / eta - 1));
        assert.ok(m.storage > prev, "the storage stiffness rises with frequency");
        prev = m.storage;
      }
      const tol = eta > 0.5 ? 0.1 : 0.06;
      assert.ok(worst < tol, `eta ${eta}, reference ${fRef} Hz: loss factor off by up to ${(100 * worst).toFixed(1)}% between 5 Hz and 10 kHz`);
      // Relaxed and glassy stiffness bracket the reference.
      const glassy = ch.kInf + ch.k.reduce((a, b) => a + b, 0);
      assert.ok(ch.kInf < 1 && glassy > 1);
    }
  }
  const none = feltChain(0, 100);
  assert.equal(none.tau.length, 0);
  assert.equal(none.kInf, 1);
  assert.deepEqual(feltModulus(none, 1000), { storage: 1, loss: 0 });
});

test("the exact step of an arm reproduces its response to a ramp and to a sinusoid", () => {
  const ch = feltChain(0.3, 50);
  const a = 6; // an arm in the middle of the band
  const tau = ch.tau[a];
  // A ramp w = s t from a relaxed arm: p(t) = s tau (1 - exp(-t / tau)), whatever the step.
  for (const dt of [tau / 50, tau / 3, 2 * tau]) {
    const { decay, growth } = feltStepConstants(ch, dt);
    const s = 1e-3;
    let p = 0;
    const steps = 20;
    for (let i = 0; i < steps; i++) p = p * decay[a] + s * dt * growth[a];
    const exact = s * tau * (1 - Math.exp(-(steps * dt) / tau));
    assert.ok(Math.abs(p / exact - 1) < 1e-12, `dt = ${dt / tau} tau: ${p} vs ${exact}`);
  }
  // A sinusoid w = A sin(omega t), stepped finely: the work per cycle against
  // the arm's force -k p is pi k A^2 omega tau / (1 + omega^2 tau^2).
  const omega = 1 / tau, A = 1e-4, steps = 2000;
  const dt = (2 * Math.PI) / omega / steps;
  const { decay, growth } = feltStepConstants(ch, dt);
  let p = 0, work = 0;
  for (let i = 0; i < 6 * steps; i++) {
    const w0 = A * Math.sin(omega * i * dt), w1 = A * Math.sin(omega * (i + 1) * dt);
    const f0 = ch.k[a] * p;
    p = p * decay[a] + (w1 - w0) * growth[a];
    const f1 = ch.k[a] * p;
    if (i >= 5 * steps) work += 0.5 * (f0 + f1) * (w1 - w0);
  }
  const expected = (Math.PI * ch.k[a] * A * A * (omega * tau)) / (1 + omega * omega * tau * tau);
  assert.ok(Math.abs(work / expected - 1) < 1e-3, `work per cycle ${work} vs ${expected}`);
});