// Screen a script for musical chaos: twin growth rate, pitch shift of the
// strongest low peak vs the lowest linear mode, and spectral flatness.
// Usage: tsx tests/screen.ts script.m duration
import { readFileSync, readdirSync } from "fs";
import { prepareSystem, runScript } from "../src/sim/runScript";
import { simulate } from "../src/sim/solver";
import { runTwin } from "../src/sim/twin";
import { growthRate, spectrum } from "../src/sim/analysis";
const helpers = readdirSync("matlab/helpers").map(name => ({ name, source: readFileSync("matlab/helpers/" + name, "utf8") }));
const { sys } = runScript(readFileSync(process.argv[2], "utf8"), helpers, () => {});
sys.duration = Number(process.argv[3]);
const notes = prepareSystem(sys);
if (notes.length) console.log(notes.join(" "));
const saved: Float64Array[] = [];
const a = simulate(sys, { onFrame: (_f, _t, x) => saved.push(x.slice()) });
const { sep } = runTwin(sys, saved, a.dt, 1e-9);
const g = growthRate(a.frameTimes, sep, 1e-9);
const fs = a.sampleRate;
const sp = spectrum(a.audio, Math.round(0.2 * fs), Math.round(0.5 * fs));
const hz = (b: number) => b * sp.df * fs;
// flatness over 50 Hz..5 kHz
let lg = 0, ar = 0, cnt = 0;
for (let b = 0; b < sp.db.length; b++) { if (hz(b) < 50 || hz(b) > 5000) continue; const p = 10 ** (sp.db[b] / 10); lg += Math.log(p); ar += p; cnt++; }
const flatness = Math.exp(lg / cnt) / (ar / cnt);
// nonlinear pitch shift: strongest peak of the same strike scaled down by 1e3, vs the nearest peak of the actual run
const small = simulate({ ...sys, hits: sys.hits.map(h => ({ ...h, velocity: h.velocity.map(v => v * 1e-3) as any, force: h.force.map(v => v * 1e-3) as any })) }, { maxAnimBytes: 0 });
const spS = spectrum(small.audio, Math.round(0.2 * fs), Math.round(0.5 * fs));
let bestS = -Infinity, fS = 0;
for (let b = 1; b < spS.db.length; b++) if (hz(b) > 50 && spS.db[b] > bestS) (bestS = spS.db[b]), (fS = hz(b));
let best = -Infinity, fb = 0;
for (let b = 1; b < sp.db.length; b++) if (hz(b) > fS * 0.89 && hz(b) < fS * 1.12 && sp.db[b] > best) (best = sp.db[b]), (fb = hz(b));
const cents = 1200 * Math.log2(fb / fS);
const tmax = sep.findIndex(s => s > 0.1);
console.log(`main peak ${fS.toFixed(0)} Hz, nonlinear shift ${cents.toFixed(0)} cents | chaos rate ${g?.toFixed(0) ?? "-"} /s, sep>0.1 at ${tmax >= 0 ? a.frameTimes[tmax].toFixed(2) + " s" : "never"} | flatness ${flatness.toExponential(1)} | drift ${a.energyDrift.toExponential(1)} | ${(a.elapsedMs / 1000).toFixed(1)} s`);