// Pair-invariant comparison of two modes files: the transfer function from a point force to the
// listener pressure, |H(w)| = |sum_m G_m phi_m(p) (-w^2) / (w0^2 - w^2 + 2 i sigma w)|, on a
// frequency grid, and the profile of the frequency differences along the spectrum.
//   node tests/transfer.ts a.bin b.bin
import { readFileSync } from "node:fs";
import { locate, modeValuesAt, parseModes, type ModesData } from "../src/modes.ts";

const sigmaOf = (d: ModesData) => {
  const mat = (d.meta.config as any).material;
  return Float64Array.from(d.dampQ, (q) => mat.sigma0 + mat.sigma1 * q);
};

const A = parseModes(readFileSync(process.argv[2]).buffer as ArrayBuffer);
const B = parseModes(readFileSync(process.argv[3]).buffer as ArrayBuffer);
const n = Math.min(A.nModes, B.nModes);
for (const [lo, hi] of [[0, 100], [100, 300], [300, 600], [600, 1000], [1000, n]]) {
  let worst = 0, mean = 0, cnt = 0;
  for (let m = lo; m < Math.min(hi, n); m++) { const r = Math.abs(A.freq[m] - B.freq[m]) / B.freq[m]; worst = Math.max(worst, r); mean += r; cnt++; }
  console.log(`modes ${lo}-${Math.min(hi, n)}: freq rel diff mean ${(mean / cnt).toExponential(2)} max ${worst.toExponential(2)}`);
}
function transfer(d: ModesData, px: number, py: number, freqs: number[]): number[] {
  const sig = sigmaOf(d);
  const loc = locate(d, px, py);
  if (!loc) throw new Error("point outside mesh");
  const phi = new Float32Array(d.nModes);
  modeValuesAt(d, loc.nodes, loc.weights, phi, 0);
  return freqs.map((f) => {
    const w = 2 * Math.PI * f;
    let re = 0, im = 0;
    for (let m = 0; m < d.nModes; m++) {
      const w0 = 2 * Math.PI * d.freq[m], s = sig[m];
      const dr = w0 * w0 - w * w, di = 2 * s * w;
      const den = dr * dr + di * di;
      const hr = dr / den, hi = -di / den;            // 1 / (dr + i di)
      const g = -w * w * phi[m];
      const gr = (d.gainFront[2 * m] - d.gainBack[2 * m]) * g, gi = (d.gainFront[2 * m + 1] - d.gainBack[2 * m + 1]) * g;
      re += gr * hr - gi * hi; im += gr * hi + gi * hr;
    }
    return Math.hypot(re, im);
  });
}
const freqs = Array.from({ length: 60 }, (_, i) => 80 * Math.pow(15000 / 80, i / 59));
for (const [px, py] of [[0.137 * Math.cos(0.3), 0.137 * Math.sin(0.3)], [0.19 * Math.cos(2), 0.19 * Math.sin(2)]]) {
  const ha = transfer(A, px, py, freqs), hb = transfer(B, px, py, freqs);
  let worst = 0, mean = 0;
  for (let i = 0; i < freqs.length; i++) { const r = Math.abs(ha[i] - hb[i]) / hb[i]; worst = Math.max(worst, r); mean += r / freqs.length; }
  console.log(`|H| at (${px.toFixed(3)}, ${py.toFixed(3)}): mean rel diff ${mean.toExponential(2)}, max ${worst.toExponential(2)}`);
  console.log("  sample |H| (A|B) at 100, 1k, 5k, 12k Hz:", [0, 22, 38, 54].map((i) => `${ha[i].toExponential(3)}|${hb[i].toExponential(3)} @${freqs[i].toFixed(0)}`).join("  "));
}