1// Screen a script for musical chaos: twin growth rate, pitch shift of the
2// strongest low peak vs the lowest linear mode, and spectral flatness.
3// Usage: tsx tests/screen.ts script.m duration
4import { readFileSync, readdirSync } from "fs";
5import { prepareSystem, runScript } from "../src/sim/runScript";
6import { simulate } from "../src/sim/solver";
7import { runTwin } from "../src/sim/twin";
8import { growthRate, spectrum } from "../src/sim/analysis";
9const helpers = readdirSync("matlab/helpers").map(name => ({ name, source: readFileSync("matlab/helpers/" + name, "utf8") }));
10const { sys } = runScript(readFileSync(process.argv[2], "utf8"), helpers, () => {});
11sys.duration = Number(process.argv[3]);
12const notes = prepareSystem(sys);
13if (notes.length) console.log(notes.join(" "));
14const saved: Float64Array[] = [];
15const a = simulate(sys, { onFrame: (_f, _t, x) => saved.push(x.slice()) });
16const { sep } = runTwin(sys, saved, a.dt, 1e-9);
17const g = growthRate(a.frameTimes, sep, 1e-9);
18const fs = a.sampleRate;
19const sp = spectrum(a.audio, Math.round(0.2 * fs), Math.round(0.5 * fs));
20const hz = (b: number) => b * sp.df * fs;
21// flatness over 50 Hz..5 kHz
22let lg = 0, ar = 0, cnt = 0;
23for (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++; }
24const flatness = Math.exp(lg / cnt) / (ar / cnt);
25// nonlinear pitch shift: strongest peak of the same strike scaled down by 1e3, vs the nearest peak of the actual run
26const 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 });
27const spS = spectrum(small.audio, Math.round(0.2 * fs), Math.round(0.5 * fs));
28let bestS = -Infinity, fS = 0;
29for (let b = 1; b < spS.db.length; b++) if (hz(b) > 50 && spS.db[b] > bestS) (bestS = spS.db[b]), (fS = hz(b));
30let best = -Infinity, fb = 0;
31for (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));
32const cents = 1200 * Math.log2(fb / fS);
33const tmax = sep.findIndex(s => s > 0.1);
34console.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`);