// WebGPU version of simulate(): same scheme and result, single precision.
// See gpu.wgsl for the kernel. The CPU prepares the rest geometry in double
// precision, chooses the step (prepare.ts), and computes the energies and
// mode energies from the positions and velocities the GPU returns at each
// frame.

import shader from "./gpu.wgsl?raw";
import { hingeAngle } from "./hinge";
import { ModeTracker, prepare, staticForces } from "./prepare";
import { frequencyErrorCents, stepCoefficients } from "./solver";
import type { SimOptions } from "./solver";
import type { BeadSystem, SimResult } from "./types";

const HDR_WORDS = 44;
const TARGET_MS = 120; // aim for dispatches of about this length (avoids GPU watchdog resets)

export async function gpuAvailable(): Promise<boolean> {
  try {
    return !!(typeof navigator !== "undefined" && navigator.gpu && (await navigator.gpu.requestAdapter()));
  } catch {
    return false;
  }
}

// incidence lists: per bead, (slot, coefficient) pairs
function csr(n: number, entries: [number, number, number][]): { rowptr: Uint32Array; data: Uint32Array } {
  const count = new Uint32Array(n + 1);
  for (const [b] of entries) count[b + 1]++;
  for (let i = 0; i < n; i++) count[i + 1] += count[i];
  const pos = count.slice(0, n);
  const data = new Uint32Array(2 * Math.max(1, entries.length));
  const f = new Float32Array(1);
  const u = new Uint32Array(f.buffer);
  for (const [b, slot, coef] of entries) {
    const q = pos[b]++;
    data[2 * q] = slot;
    f[0] = coef;
    data[2 * q + 1] = u[0];
  }
  return { rowptr: count, data };
}

export async function simulateGPU(sys: BeadSystem, opts: SimOptions = {}): Promise<SimResult> {
  const t0 = performance.now();
  const adapter = await navigator.gpu.requestAdapter();
  if (!adapter) throw new Error("WebGPU: no adapter");
  const device = await adapter.requestDevice();
  try {
    return await run(device, sys, opts, t0);
  } finally {
    device.destroy();
  }
}

async function run(device: GPUDevice, sys: BeadSystem, opts: SimOptions, t0: number): Promise<SimResult> {
  const { n, pos0, mass, pinned } = sys;
  const prep = prepare(sys);
  const warnings = [...prep.warnings];
  const { substeps, h, nSamples, peRest } = prep;
  const fs = sys.sampleRate;
  const Delta = 1 / fs;
  const cs = stepCoefficients(sys.method);
  const ns = sys.si.length;
  const nb = sys.bi.length;
  const nh = sys.hi.length;
  const keepD = Array.from(sys.dc.keys()).filter(d => sys.dc[d] !== 0 || sys.dca[d] !== 0);
  const nd = keepD.length;

  // ---- parameters (f32) ------------------------------------------------------
  const pulses = sys.hits.filter(hh => hh.duration > 0);
  const oInvm = 0;
  const oMass = n;
  const oFric = 2 * n;
  const oS = 3 * n;
  const oB = oS + 8 * ns;
  const oH = oB + 8 * nb;
  const oD = oH + 16 * nh;
  const oPulse = oD + 8 * nd;
  const oCs = oPulse + 8 * pulses.length;
  const P = new Float32Array(oCs + cs.length + 1);
  for (let i = 0; i < n; i++) {
    P[oInvm + i] = pinned[i] ? 0 : 1 / mass[i];
    P[oMass + i] = pinned[i] ? 0 : mass[i];
    P[oFric + i] = sys.friction[i];
  }
  const rel = (a: number, b: number, c: number) => pos0[3 * b + c] - pos0[3 * a + c]; // x_b - x_a at rest
  for (let s = 0; s < ns; s++) {
    const i = sys.si[s];
    const j = sys.sj[s];
    const r = [rel(i, j, 0), rel(i, j, 1), rel(i, j, 2)];
    const L = Math.hypot(r[0], r[1], r[2]);
    P.set([r[0], r[1], r[2], L, L - sys.sL0[s], sys.sk[s], sys.sk3[s]], oS + 8 * s);
  }
  for (let t = 0; t < nb; t++) {
    const [i, j, k] = [sys.bi[t], sys.bj[t], sys.bk[t]];
    const D = [0, 1, 2].map(c => pos0[3 * i + c] - 2 * pos0[3 * j + c] + pos0[3 * k + c]);
    P.set([D[0], D[1], D[2], sys.bkb[t], sys.bcb[t]], oB + 8 * t);
  }
  for (let q = 0; q < nh; q++) {
    const [i, j, k, l] = [sys.hi[q], sys.hj[q], sys.hk[q], sys.hl[q]];
    const R = [i, k, l].flatMap(b => [rel(j, b, 0), rel(j, b, 1), rel(j, b, 2)]);
    P.set([...R, sys.hkh[q], sys.hch[q]], oH + 16 * q);
  }
  keepD.forEach((d, q) => {
    const [i, j] = [sys.di[d], sys.dj[d]];
    P.set([rel(i, j, 0), rel(i, j, 1), rel(i, j, 2), sys.dc[d], sys.dca[d]], oD + 8 * q);
  });
  pulses.forEach((p, q) => P.set([p.t, p.duration, ...p.force], oPulse + 8 * q));
  P.set(cs, oCs);

  // ---- indices (u32) -----------------------------------------------------------
  const eFb = 3 * ns;
  const eFh = eFb + 3 * nb;
  const eDb = 3 * nd;
  const eDh = eDb + 3 * nb;
  const forceEntries: [number, number, number][] = [];
  const dampEntries: [number, number, number][] = [];
  for (let s = 0; s < ns; s++) {
    forceEntries.push([sys.si[s], 3 * s, 1], [sys.sj[s], 3 * s, -1]);
  }
  for (let t = 0; t < nb; t++) {
    forceEntries.push([sys.bi[t], eFb + 3 * t, -1], [sys.bj[t], eFb + 3 * t, 2], [sys.bk[t], eFb + 3 * t, -1]);
    if (sys.bcb[t] !== 0) dampEntries.push([sys.bi[t], eDb + 3 * t, 1], [sys.bj[t], eDb + 3 * t, -2], [sys.bk[t], eDb + 3 * t, 1]);
  }
  for (let q = 0; q < nh; q++) {
    [sys.hi[q], sys.hj[q], sys.hk[q], sys.hl[q]].forEach((b, p) => {
      forceEntries.push([b, eFh + 12 * q + 3 * p, -1]);
      if (sys.hch[q] !== 0) dampEntries.push([b, eDh + 12 * q + 3 * p, 1]);
    });
  }
  keepD.forEach((d, q) => dampEntries.push([sys.di[d], 3 * q, 1], [sys.dj[d], 3 * q, -1]));
  const fcsr = csr(n, forceEntries);
  const dcsr = csr(n, dampEntries);
  const pulseInc: number[][] = Array.from({ length: n }, () => []);
  pulses.forEach((p, q) => p.beads.forEach(b => pulseInc[b].push(q)));
  const prow = new Uint32Array(n + 1);
  for (let i = 0; i < n; i++) prow[i + 1] = prow[i] + pulseInc[i].length;
  const pdat = Uint32Array.from(pulseInc.flat());
  // impulses, grouped by audio sample
  const impulses = new Map<number, [number, number, number, number][]>();
  for (const hit of sys.hits) {
    if (hit.duration > 0) continue;
    const s = Math.round(hit.t * fs);
    if (s >= nSamples) continue;
    if (!impulses.has(s)) impulses.set(s, []);
    for (const b of hit.beads) impulses.get(s)!.push([b, ...hit.velocity]);
  }
  const impSamples = [...impulses.keys()].sort((a, b) => a - b);
  const impList: [number, number, number, number][] = [];
  const impStart = new Map<number, number>();
  for (const s of impSamples) {
    impStart.set(s, impList.length);
    impList.push(...impulses.get(s)!);
  }

  const parts: Uint32Array[] = [];
  let off = 0;
  const place = (a: Uint32Array) => {
    const o = off;
    parts.push(a);
    off += a.length;
    return o;
  };
  const iS = place(new Uint32Array(Array.from({ length: ns }, (_, s) => [sys.si[s], sys.sj[s]]).flat()));
  const iB = place(new Uint32Array(Array.from({ length: nb }, (_, t) => [sys.bi[t], sys.bj[t], sys.bk[t]]).flat()));
  const iH = place(new Uint32Array(Array.from({ length: nh }, (_, q) => [sys.hi[q], sys.hj[q], sys.hk[q], sys.hl[q]]).flat()));
  const iD = place(new Uint32Array(keepD.flatMap(d => [sys.di[d], sys.dj[d]])));
  const iFr = place(fcsr.rowptr);
  const iFe = place(fcsr.data);
  const iDr = place(dcsr.rowptr);
  const iDe = place(dcsr.data);
  const iPr = place(prow);
  const iPe = place(pdat.length ? pdat : new Uint32Array(1));
  const iL = place(Uint32Array.from(sys.listenBeads));
  const impWords = new Uint32Array(4 * Math.max(1, impList.length));
  {
    const f = new Float32Array(impWords.buffer);
    impList.forEach(([b, vx, vy, vz], q) => {
      impWords[4 * q] = b;
      f[4 * q + 1] = vx;
      f[4 * q + 2] = vy;
      f[4 * q + 3] = vz;
    });
  }
  const iImp = place(impWords);
  const I = new Uint32Array(off + 1);
  {
    let o = 0;
    for (const a of parts) {
      I.set(a, o);
      o += a.length;
    }
  }

  // ---- damping: number of Taylor terms ------------------------------------------
  let normA = 0;
  {
    const row = new Float64Array(n);
    for (let i = 0; i < n; i++) row[i] = sys.friction[i];
    keepD.forEach(d => {
      const c = 2 * (sys.dc[d] + sys.dca[d]);
      row[sys.di[d]] += c;
      row[sys.dj[d]] += c;
    });
    for (let t = 0; t < nb; t++) {
      row[sys.bi[t]] += 4 * sys.bcb[t];
      row[sys.bj[t]] += 8 * sys.bcb[t];
      row[sys.bk[t]] += 4 * sys.bcb[t];
    }
    const g = new Float64Array(12);
    for (let q = 0; q < nh; q++) {
      if (sys.hch[q] === 0) continue;
      hingeAngle(pos0, sys.hi[q], sys.hj[q], sys.hk[q], sys.hl[q], g);
      const total = g.reduce((a, x) => a + Math.abs(x), 0);
      // 1.5x margin: the gradients change as the shell moves
      [sys.hi[q], sys.hj[q], sys.hk[q], sys.hl[q]].forEach((b, p) => {
        row[b] += 1.5 * sys.hch[q] * total * (Math.abs(g[3 * p]) + Math.abs(g[3 * p + 1]) + Math.abs(g[3 * p + 2]));
      });
    }
    for (let i = 0; i < n; i++) if (!pinned[i]) normA = Math.max(normA, row[i] / mass[i]);
  }
  const hasDamp = normA > 0;
  const nSplit = Math.max(1, Math.ceil((Delta * normA) / 0.5));
  let K = 1;
  {
    const z = (Delta / nSplit) * normA;
    let term = z;
    while (term > 1e-9 && K < 30) {
      K++;
      term *= z / K;
    }
  }

  // ---- buffers ---------------------------------------------------------------------
  const frameStride = Math.max(1, Math.round(fs / sys.frameRate));
  const nFrames = Math.floor((nSamples - 1) / frameStride) + 1;
  const frameSize = 7 * n; // u, v, dissipated energy per bead
  const maxChunk = Math.min(nSamples, 8192);
  const framesPerChunk = Math.floor(maxChunk / frameStride) + 2;
  const oFrames = maxChunk;
  const outWords = oFrames + framesPerChunk * frameSize;
  const usage = GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST | GPUBufferUsage.COPY_SRC;
  const mk = (words: number, u = usage) => device.createBuffer({ size: Math.max(16, 4 * words), usage: u });
  const bufH = device.createBuffer({ size: 4 * HDR_WORDS, usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST });
  const bufP = mk(P.length);
  const bufI = mk(I.length);
  const bufS = mk(17 * n + nh + 1);
  const bufE = mk(Math.max(3 * ns + 3 * nb + 12 * nh, 3 * nd + 3 * nb + 12 * nh, 1));
  const bufO = mk(outWords);
  const bufC = mk(4);
  const staging = mk(outWords, GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST);
  device.queue.writeBuffer(bufP, 0, P);
  device.queue.writeBuffer(bufI, 0, I);

  // per-phase constants, selected with a dynamic offset
  const align = device.limits.minUniformBufferOffsetAlignment;
  const phases: { kick: number; drift: number; tOff: number; gather: number; c: number; saveAv: number; last: number }[] = [];
  const addPhase = (p: Partial<(typeof phases)[number]>) => {
    phases.push({ kick: 0, drift: 0, tOff: 0, gather: 0, c: 0, saveAv: 0, last: 0, ...p });
    return (phases.length - 1) * align;
  };
  // During the sub-steps of sample s the sample index on the GPU is already
  // s + 1 (the record phase advanced it), so t = (s + 1) Delta + tOff.
  const m = substeps * cs.length;
  const csAt = (j: number) => cs[j % cs.length];
  const phStart = addPhase({ kick: 0.5 * csAt(0) * h, drift: csAt(0) * h });
  const phMid: number[] = [];
  let tAcc = csAt(0) * h;
  for (let j = 1; j < m; j++) {
    phMid.push(addPhase({ gather: 1, kick: 0.5 * (csAt(j - 1) + csAt(j)) * h, drift: csAt(j) * h, tOff: tAcc - Delta }));
    tAcc += csAt(j) * h;
  }
  const phEnd = addPhase({ gather: 1, kick: 0.5 * csAt(m - 1) * h, tOff: tAcc - Delta });
  const phInit = addPhase({ gather: 1, tOff: -Delta });
  const phDamp: number[][] = [];
  for (let p = 0; p < nSplit; p++) {
    phDamp.push(
      Array.from({ length: K }, (_, k) =>
        addPhase({ c: -Delta / nSplit / (k + 1), saveAv: p === 0 && k === 0 ? 1 : 0, last: p === nSplit - 1 && k === K - 1 ? 1 : 0 })
      )
    );
  }
  const bufQ = device.createBuffer({ size: phases.length * align, usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST });
  {
    const q = new ArrayBuffer(phases.length * align);
    phases.forEach((ph, r) => {
      const f = new Float32Array(q, r * align, 8);
      const u = new Uint32Array(q, r * align, 8);
      f[0] = ph.kick;
      f[1] = ph.drift;
      f[2] = ph.tOff;
      u[3] = ph.gather;
      f[4] = ph.c;
      u[5] = ph.saveAv;
      u[6] = ph.last;
    });
    device.queue.writeBuffer(bufQ, 0, q);
  }

  const module = device.createShaderModule({ code: shader });
  const info = await module.getCompilationInfo();
  const errs = info.messages.filter(msg => msg.type === "error");
  if (errs.length) throw new Error("WebGPU shader: " + errs.map(msg => `${msg.lineNum}:${msg.linePos} ${msg.message}`).join("; "));
  const storage = (binding: number, type: GPUBufferBindingType) => ({ binding, visibility: GPUShaderStage.COMPUTE, buffer: { type } });
  const layout0 = device.createBindGroupLayout({
    entries: [
      { binding: 0, visibility: GPUShaderStage.COMPUTE, buffer: { type: "uniform" as const } },
      storage(1, "read-only-storage"),
      storage(2, "read-only-storage"),
      storage(3, "storage"),
      storage(4, "storage"),
      storage(5, "storage"),
      storage(6, "storage"),
    ],
  });
  const layout1 = device.createBindGroupLayout({
    entries: [{ binding: 0, visibility: GPUShaderStage.COMPUTE, buffer: { type: "uniform", hasDynamicOffset: true, minBindingSize: 32 } }],
  });
  const plLayout = device.createPipelineLayout({ bindGroupLayouts: [layout0, layout1] });
  const pipe = (entryPoint: string) => device.createComputePipeline({ layout: plLayout, compute: { module, entryPoint } });
  const K_ = {
    elements: pipe("elements"),
    beads: pipe("beads"),
    pre: pipe("pre"),
    dampElements: pipe("dampElements"),
    dampBeads: pipe("dampBeads"),
    record: pipe("record"),
    impulse: pipe("impulse"),
    initHinges: pipe("initHinges"),
    resetT: pipe("resetT"),
  };
  const bind0 = device.createBindGroup({
    layout: layout0,
    entries: [bufH, bufP, bufI, bufS, bufE, bufO, bufC].map((buffer, binding) => ({ binding, resource: { buffer } })),
  });
  const bind1 = device.createBindGroup({ layout: layout1, entries: [{ binding: 0, resource: { buffer: bufQ, size: 32 } }] });

  const hdr = new ArrayBuffer(4 * HDR_WORDS);
  const hu = new Uint32Array(hdr);
  const hf = new Float32Array(hdr);
  const quantity = { displacement: 0, velocity: 1, acceleration: 2 }[sys.listenQuantity];
  const setHdr = (sample0: number, nImp: number, impOff: number) => {
    hu.set([n, ns, nb, nh, nd, sys.listenBeads.length, quantity, pulses.length]);
    hu.set([sample0, frameStride, frameSize, nImp, impOff, oFrames, 0, 0], 8);
    hu.set([oInvm, oMass, oFric, oS, oB, oH, oD, oPulse], 16);
    hu.set([iS, iB, iH, iD, iFr, iFe, iDr, iDe, iPr, iPe, iL, iImp], 24);
    hu.set([eFb, eFh, eDb, eDh], 36);
    hf.set([Delta, ...sys.listenDir], 40);
    device.queue.writeBuffer(bufH, 0, hdr);
  };
  const groups = (count: number) => Math.max(1, Math.ceil(count / 64));
  const gBeads = groups(n);
  const gElems = groups(ns + nb + nh);
  const gDamp = groups(nd + nb + nh);
  // only the kernels that read the phase table need it rebound
  let boundPhase = -1;
  const go = (pass: GPUComputePassEncoder, p: GPUComputePipeline, nGroups: number, phase = -1) => {
    pass.setPipeline(p);
    if (phase >= 0 && phase !== boundPhase) {
      pass.setBindGroup(1, bind1, [phase]);
      boundPhase = phase;
    }
    pass.dispatchWorkgroups(nGroups);
  };
  const encodeSample = (pass: GPUComputePassEncoder) => {
    go(pass, K_.pre, gBeads);
    if (quantity !== 2 || !hasDamp) go(pass, K_.record, 1);
    if (hasDamp) {
      phDamp.forEach((terms, p) => {
        if (p > 0) go(pass, K_.resetT, gBeads); // T = v again (pieces of a long damping step)
        terms.forEach((ph, k) => {
          go(pass, K_.dampElements, gDamp);
          go(pass, K_.dampBeads, gBeads, ph);
          if (quantity === 2 && p === 0 && k === 0) go(pass, K_.record, 1);
        });
      });
    }
    go(pass, K_.beads, gBeads, phStart);
    for (const ph of phMid) {
      go(pass, K_.elements, gElems);
      go(pass, K_.beads, gBeads, ph);
    }
    go(pass, K_.elements, gElems);
    go(pass, K_.beads, gBeads, phEnd);
  };

  // initial state: rest angles of the hinges, then the forces at t = 0
  setHdr(0, 0, 0);
  device.queue.writeBuffer(bufC, 0, Uint32Array.of(1));
  {
    const enc = device.createCommandEncoder();
    const pass = enc.beginComputePass();
    pass.setBindGroup(0, bind0);
    pass.setBindGroup(1, bind1, [0]);
    boundPhase = 0;
    go(pass, K_.initHinges, groups(nh));
    go(pass, K_.elements, gElems);
    go(pass, K_.beads, gBeads, phInit);
    pass.end();
    device.queue.submit([enc.finish()]);
  }

  // ---- outputs --------------------------------------------------------------------
  const audio = new Float32Array(nSamples);
  const frameTimes = new Float64Array(nFrames);
  const energy = new Float64Array(nFrames);
  const dissipated = new Float64Array(nFrames);
  const maxAnimBytes = opts.maxAnimBytes ?? 150e6;
  const animStride = frameStride * Math.max(1, Math.ceil((nFrames * 3 * n * 4) / Math.max(maxAnimBytes, 1)));
  const animTimes: number[] = [];
  const animFrames: Float32Array[] = [];
  const modes = opts.modes ?? null;
  const nModes = modes ? modes.freqs.length : 0;
  const modeEnergy = modes ? new Float32Array(nFrames * nModes) : null;
  const tracker = modes ? new ModeTracker(modes, mass, pos0) : null;
  const xf = new Float64Array(3 * n);
  const vf = new Float64Array(3 * n);
  const Fs = new Float64Array(3 * n);
  let frame = 0;
  let aborted: string | null = null;

  let sample = 0;
  let chunk = 16;
  let lastProgress = performance.now();
  while (sample < nSamples) {
    let len = Math.min(chunk, maxChunk, nSamples - sample);
    const nextImp = impSamples.find(s => s > sample);
    if (nextImp !== undefined) len = Math.min(len, nextImp - sample);
    const imps = impulses.get(sample);
    const tStart = performance.now();
    setHdr(sample, imps ? imps.length : 0, imps ? impStart.get(sample)! : 0);
    device.queue.writeBuffer(bufC, 0, Uint32Array.of(sample));
    const enc = device.createCommandEncoder();
    const pass = enc.beginComputePass();
    pass.setBindGroup(0, bind0);
    pass.setBindGroup(1, bind1, [0]);
    boundPhase = 0;
    if (imps) go(pass, K_.impulse, groups(imps.length));
    for (let q = 0; q < len; q++) encodeSample(pass);
    pass.end();
    enc.copyBufferToBuffer(bufO, 0, staging, 0, 4 * outWords);
    device.queue.submit([enc.finish()]);
    await staging.mapAsync(GPUMapMode.READ);
    const out = new Float32Array(staging.getMappedRange().slice(0));
    staging.unmap();
    const elapsed = performance.now() - tStart;

    audio.set(out.subarray(0, len), sample);
    for (let q = 0; q < len; q++) {
      if (!Number.isFinite(out[q])) {
        aborted = `The solution blew up near t = ${((sample + q) / fs).toFixed(4)} s. Reduce dt (or safety) or the strike strength.`;
        break;
      }
    }
    // frames in this chunk
    const first = Math.ceil(sample / frameStride) * frameStride;
    for (let fsmp = first; fsmp < sample + len && frame < nFrames && !aborted; fsmp += frameStride) {
      const o = oFrames + ((fsmp - first) / frameStride) * frameSize;
      let ke = 0;
      for (let q = 0; q < 3 * n; q++) {
        xf[q] = pos0[q] + out[o + q];
        vf[q] = out[o + 3 * n + q];
      }
      for (let i = 0; i < n; i++) ke += 0.5 * (pinned[i] ? 0 : mass[i]) * (vf[3 * i] ** 2 + vf[3 * i + 1] ** 2 + vf[3 * i + 2] ** 2);
      frameTimes[frame] = fsmp / fs;
      energy[frame] = ke + staticForces(sys, xf, Fs) - peRest;
      let diss = 0;
      for (let i = 0; i < n; i++) diss += out[o + 6 * n + i];
      dissipated[frame] = diss;
      if (tracker && modeEnergy) tracker.record(out.subarray(o, o + 3 * n), out.subarray(o + 3 * n, o + 6 * n), modeEnergy, frame * nModes, true);
      if (fsmp % animStride === 0) {
        animTimes.push(fsmp / fs);
        animFrames.push(out.slice(o, o + 3 * n));
      }
      opts.onFrame?.(frame, fsmp / fs, xf, vf);
      frame++;
    }
    if (aborted) break;
    sample += len;
    // adapt the chunk so each dispatch takes about TARGET_MS
    chunk = Math.max(16, Math.min(maxChunk, Math.round((len * TARGET_MS) / Math.max(elapsed, 1))));
    const now = performance.now();
    if (opts.onProgress && now - lastProgress > 200) {
      opts.onProgress(sample / fs);
      lastProgress = now;
    }
  }

  // energy drift after the last excitation
  const tExcitationEnd = Math.max(0, ...sys.hits.map(hh => hh.t + hh.duration));
  let energyDrift = 0;
  let ref = -1;
  for (let f = 0; f < frame; f++) {
    if (frameTimes[f] < tExcitationEnd) continue;
    const total = energy[f] + dissipated[f];
    if (ref < 0) ref = total;
    else if (ref > 0) energyDrift = Math.max(energyDrift, Math.abs(total - ref) / ref);
  }

  return {
    audio,
    sampleRate: fs,
    dt: h,
    substeps,
    omegaMax: prep.omegaMax,
    frameTimes: frameTimes.slice(0, frame),
    energy: energy.slice(0, frame),
    dissipated: dissipated.slice(0, frame),
    animTimes: Float64Array.from(animTimes),
    animFrames,
    modeEnergy: modeEnergy ? modeEnergy.slice(0, frame * nModes) : null,
    residualForce: prep.residualForce,
    maxTension: prep.maxTension,
    energyDrift,
    freqErrors: [1000, 2000, 5000, 10000, 20000]
      .filter(hz => hz < fs / 2)
      .map(hz => ({ hz, cents: frequencyErrorCents(sys.method, 2 * Math.PI * hz, h) })),
    elapsedMs: performance.now() - t0,
    aborted,
    warnings,
  };
}