// Perturbed twin: rerun with every hit scaled by (1 + eps) and measure, at
// each frame, the largest difference in any bead coordinate, relative to the
// largest displacement of the original run over the preceding 50 ms (so that
// the measure stays meaningful while the motion decays).

import { simulate } from "./solver";
import type { BeadSystem } from "./types";

export function runTwin(
  sys: BeadSystem,
  saved: Float64Array[], // positions of the original run at each frame
  dt: number,
  eps: number,
  onProgress?: (t: number) => void
): { audio: Float32Array; sep: Float64Array } {
  const disp = saved.map(x => {
    let d = 0;
    for (let q = 0; q < x.length; q++) d = Math.max(d, Math.abs(x[q] - sys.pos0[q]));
    return d;
  });
  const win = Math.max(1, Math.round(0.05 * sys.frameRate));
  const scale = disp.map((_, f) => Math.max(...disp.slice(Math.max(0, f - win), f + 1)));
  const sep = new Float64Array(saved.length);
  const f = 1 + eps;
  const r = simulate(
    {
      ...sys,
      dt,
      hits: sys.hits.map(h => ({
        ...h,
        velocity: [h.velocity[0] * f, h.velocity[1] * f, h.velocity[2] * f],
        force: [h.force[0] * f, h.force[1] * f, h.force[2] * f],
      })),
    },
    {
      maxAnimBytes: 0,
      onProgress,
      onFrame: (fr, _t, x) => {
        const y = saved[fr];
        if (!y) return;
        let d = 0;
        for (let q = 0; q < x.length; q++) d = Math.max(d, Math.abs(x[q] - y[q]));
        sep[fr] = scale[fr] > 0 ? d / scale[fr] : 0;
      },
    }
  );
  return { audio: r.audio, sep };
}