import { describe, expect, it } from 'vitest';
import { buildCymbal, buildCymbalExcitation, buildMic, deriveCymbal, type CymbalParams } from '../src/instruments/bead-cymbal/model';
import { CYMBAL_PRESETS } from '../src/instruments/bead-cymbal/presets';
import { bendingEnergy, simulate } from '../src/sim/solver';
import { spectrum } from '../src/ui/fft';

// the flat plate (the tests compare with flat-plate theory)
const flat = CYMBAL_PRESETS.find((p) => p.name.startsWith('Flat 8" toy, stick'))!.params;
const base = (over: Partial<CymbalParams> = {}): CymbalParams => ({ ...flat, ...over });

describe('cymbal lattice', () => {
  it('bending energy matches the plate (nu = 1/3) for bending, doming and twisting', () => {
    const p = base({ holeRadius: 0, rings: 12 });
    const d = deriveCymbal(p);
    const sys = buildCymbal(d);
    const k = 0.5; // 1/m, small: k * a << 1
    const R = p.diameter / 2;
    const nu = 1 / 3;
    const inner = (t: number) => d.mesh.tris[t].every((i) => Math.hypot(sys.rest[3 * i], sys.rest[3 * i + 2]) < 0.6 * R);
    let area = 0;
    d.mesh.tris.forEach((tr, t) => {
      if (!inner(t)) return;
      const [u, v, w] = tr.map((i) => [sys.rest[3 * i], sys.rest[3 * i + 2]]);
      area += 0.5 * Math.abs((v[0] - u[0]) * (w[1] - u[1]) - (w[0] - u[0]) * (v[1] - u[1]));
    });
    const shapes: [string, (x: number, z: number) => number, number][] = [
      ['cylinder', (x) => 0.5 * k * x * x, 0.5 * k * k],
      ['sphere', (x, z) => 0.5 * k * (x * x + z * z), 0.5 * (4 - 2 * (1 - nu)) * k * k],
      ['saddle', (x, z) => 0.5 * k * (x * x - z * z), (1 - nu) * k * k],
      ['twist', (x, z) => k * x * z, (1 - nu) * k * k],
    ];
    for (const [name, f, density] of shapes) {
      const X = new Float64Array(3 * (sys.n + sys.anchors.length));
      X.set(sys.rest);
      for (let i = 0; i < sys.n; i++) X[3 * i + 1] = f(sys.rest[3 * i], sys.rest[3 * i + 2]);
      const ratio = bendingEnergy(sys.bending!, X, inner) / (d.D * density * area);
      console.log(name, 'bending energy / plate', ratio.toFixed(4));
      expect(ratio).toBeCloseTo(1, 2);
    }
  });

  it('a free flat plate rings near the (2,0) mode of a free circular plate', () => {
    // no felts, no hole, no damping; initial shape r^2 cos(2 phi), small
    const p = base({ holeRadius: 0, rings: 9, t60Low: 0, t60High: 0, duration: 0.3, frameRate: 50 });
    const d = deriveCymbal(p);
    const sys = buildCymbal(d);
    const exc = buildCymbalExcitation({ ...p, excitation: 'strike', strikeVel: 0 }, d);
    for (let i = 0; i < sys.n; i++) {
      const x = sys.rest[3 * i], z = sys.rest[3 * i + 2];
      exc.u0[3 * i + 1] = 1e-7 * (x * x - z * z) / 0.01;
    }
    const pick = new Float64Array(3 * sys.n);
    let edge = 0;
    for (let i = 0; i < sys.n; i++) if (sys.rest[3 * i] > sys.rest[3 * edge]) edge = i;
    pick[3 * edge + 1] = 1;
    const r = simulate(sys, exc, { duration: p.duration, fs: p.fs, cfl: p.cfl, frameRate: p.frameRate, pickup: pick });
    const s = spectrum(r.audio, r.fs);
    let best = 1;
    for (let k = 1; k < s.freq.length && s.freq[k] < 300; k++) if (s.db[k] > s.db[best]) best = k;
    const R = p.diameter / 2;
    const ideal = (5.253 / (2 * Math.PI * R * R)) * d.bendingCoef; // Leissa, nu ~ 1/3
    console.log('free (2,0) measured', s.freq[best].toFixed(1), 'Hz, ideal', ideal.toFixed(1), 'Hz');
    expect(Math.abs(s.freq[best] / ideal - 1)).toBeLessThan(0.05);
  }, 60000);

  it('energy balances with a glancing stick hit, friction, damping and rigid motion removed', () => {
    const p = base({ rings: 7, hitAngle: 60, stickSpeed: 6, bowHeight: 0.01, bellHeight: 0.01, duration: 0.03, frameRate: 2000 });
    const d = deriveCymbal(p);
    const mic = buildMic(p, d);
    const r = simulate(buildCymbal(d), buildCymbalExcitation(p, d), { duration: p.duration, fs: p.fs, cfl: p.cfl, frameRate: p.frameRate, pickup: mic.weights, pickupDelay: mic.delay });
    const emax = Math.max(...r.energy);
    const bal = Math.max(...r.balance.map(Math.abs));
    const cf = r.contactForce!;
    console.log('contact samples', cf.reduce((c, v) => c + (v > 0 ? 1 : 0), 0), 'peak N', Math.max(...cf).toFixed(0), 'balance', (bal / emax).toExponential(2));
    expect(bal / emax).toBeLessThan(1e-2);
    expect(cf[cf.length - 1]).toBe(0);
    expect(r.audio.some((v) => v !== 0)).toBe(true);
  });

  it('removing rigid motion leaves the plate with no net momentum and keeps the hit point vibrating', () => {
    const p = base({ rings: 7, duration: 0.05, frameRate: 2000 });
    const d = deriveCymbal(p);
    const sys = buildCymbal(d);
    const mic = buildMic(p, d);
    const r = simulate(sys, buildCymbalExcitation(p, d), { duration: p.duration, fs: p.fs, cfl: p.cfl, frameRate: p.frameRate, pickup: mic.weights, pickupDelay: mic.delay });
    // mean vertical displacement (rigid translation) stays near zero compared with the vibration
    const n = sys.n;
    let meanMax = 0, ampMax = 0;
    for (let f = 0; f < r.frameTimes.length; f++) {
      let m = 0, mt = 0;
      for (let i = 0; i < n; i++) {
        m += sys.mass[i] * r.frames[f * 3 * n + 3 * i + 1];
        mt += sys.mass[i];
        ampMax = Math.max(ampMax, Math.abs(r.frames[f * 3 * n + 3 * i + 1]));
      }
      meanMax = Math.max(meanMax, Math.abs(m / mt));
    }
    console.log('mean y', (meanMax * 1000).toFixed(4), 'mm, max |y|', (ampMax * 1000).toFixed(3), 'mm');
    expect(meanMax / ampMax).toBeLessThan(1e-2);
    expect(Math.max(...r.balance.map(Math.abs)) / Math.max(...r.energy)).toBeLessThan(1e-2);
  }, 60000);
});