1// Separation of two runs whose hits differ by a factor (1 + eps).
2// Usage: tsx tests/twin.ts script.m duration [eps]
3import { readFileSync, readdirSync } from "fs";
4import { runScript } from "../src/sim/runScript";
5import { simulate } from "../src/sim/solver";
6const helpers = readdirSync("matlab/helpers").map(name => ({ name, source: readFileSync("matlab/helpers/" + name, "utf8") }));
7const { sys } = runScript(readFileSync(process.argv[2], "utf8"), helpers, s => process.stdout.write(s));
8sys.duration = Number(process.argv[3]);
9const eps = Number(process.argv[4] ?? 1e-9);
10sys.frameRate = 1000;
11const xs: Float64Array[] = [];
12const a = simulate(sys, { onFrame: (_f, _t, x) => xs.push(x.slice()) });
13const sep: number[] = [];
14let scale = 0;
15for (const x of xs) for (let q = 0; q < x.length; q++) scale = Math.max(scale, Math.abs(x[q] - sys.pos0[q]));
16simulate({ ...sys, hits: sys.hits.map(h => ({ ...h, velocity: h.velocity.map(v => v * (1 + eps)) as any, force: h.force.map(v => v * (1 + eps)) as any })) },
17 { onFrame: (f, _t, x) => { let d = 0; for (let q = 0; q < x.length; q++) d = Math.max(d, Math.abs(x[q] - xs[f][q])); sep.push(d / scale); } });
18const out: string[] = [];
19const step = Math.max(1, Math.floor(sep.length / 14));
20for (let i = 0; i < sep.length; i += step) out.push(`${a.frameTimes[i].toFixed(2)}:${Math.log10(sep[i] + 1e-30).toFixed(1)}`);
21console.log(`dt ${a.dt.toExponential(2)} drift ${a.energyDrift.toExponential(1)} ${(a.elapsedMs/1000).toFixed(1)}s | log10 sep: ${out.join(" ")}`);