import type { BeadSystem, Bending, Excitation, Mallet, SimOptions, SimResult } from './types';
import { hingeAngleGrad } from './hinge';

// Normal speed of a mallet into the surface, and its linearized contact
// stiffness at the largest compression a rigid surface would see.
function contactStiffness(m: Mallet): number {
  const vn = m.vel0[0] * m.normal[0] + m.vel0[1] * m.normal[1] + m.vel0[2] * m.normal[2];
  if (!(vn > 0)) return 0;
  const dmax = (((m.p + 1) * m.mass * vn * vn) / (2 * m.K)) ** (1 / (m.p + 1));
  return m.p * m.K * dmax ** (m.p - 1);
}

function nodePositions(sys: BeadSystem): Float64Array {
  const X = new Float64Array(3 * (sys.n + sys.anchors.length));
  X.set(sys.rest);
  sys.anchors.forEach((a, k) => X.set(a.pos, 3 * (sys.n + k)));
  return X;
}

// Linearized stiffness at rest, y = K x (x, y: 3n, beads only; anchors held).
// A spring under tension T has stiffness k along its axis and T/l across it;
// bending contributes G^T M G per triangle (G = hinge-angle gradients).
function restStiffness(sys: BeadSystem): (x: Float64Array, y: Float64Array) => void {
  const { n, si, sj, sk, sL0 } = sys;
  const X = nodePositions(sys);
  const ns = si.length;
  const sd = new Float64Array(3 * ns); // unit axis
  const ka = new Float64Array(ns); // axial stiffness minus transverse
  const kt = new Float64Array(ns); // transverse T / l
  for (let s = 0; s < ns; s++) {
    const a = 3 * si[s], b = 3 * sj[s];
    const dx = X[b] - X[a], dy = X[b + 1] - X[a + 1], dz = X[b + 2] - X[a + 2];
    const l = Math.sqrt(dx * dx + dy * dy + dz * dz);
    sd[3 * s] = dx / l;
    sd[3 * s + 1] = dy / l;
    sd[3 * s + 2] = dz / l;
    kt[s] = (sk[s] * (l - sL0[s])) / l;
    ka[s] = sk[s] - kt[s];
  }
  const b = sys.bending;
  const nh = b ? b.w0.length : 0;
  const G = new Float64Array(12 * nh);
  const nodes = new Int32Array(4 * nh);
  if (b) {
    const g = new Float64Array(13);
    for (let e = 0; e < nh; e++) {
      hingeAngleGrad(X, b.w0[e], b.w1[e], b.e0[e], b.e1[e], g);
      for (let q = 0; q < 12; q++) G[12 * e + q] = b.sign[e] * g[1 + q];
      nodes.set([b.w0[e], b.w1[e], b.e0[e], b.e1[e]], 4 * e);
    }
  }
  const phi = new Float64Array(nh);
  const mom = new Float64Array(nh);
  const get = (x: Float64Array, node: number, c: number) => (node < n ? x[3 * node + c] : 0);
  return (x, y) => {
    y.fill(0);
    for (let s = 0; s < ns; s++) {
      const a = si[s], c = sj[s];
      const rx = get(x, c, 0) - get(x, a, 0), ry = get(x, c, 1) - get(x, a, 1), rz = get(x, c, 2) - get(x, a, 2);
      const ex = sd[3 * s], ey = sd[3 * s + 1], ez = sd[3 * s + 2];
      const pa = ka[s] * (ex * rx + ey * ry + ez * rz);
      const fx = pa * ex + kt[s] * rx, fy = pa * ey + kt[s] * ry, fz = pa * ez + kt[s] * rz;
      if (a < n) {
        y[3 * a] -= fx;
        y[3 * a + 1] -= fy;
        y[3 * a + 2] -= fz;
      }
      if (c < n) {
        y[3 * c] += fx;
        y[3 * c + 1] += fy;
        y[3 * c + 2] += fz;
      }
    }
    if (!b) return;
    for (let e = 0; e < nh; e++) {
      let v = 0;
      for (let k = 0; k < 4; k++) {
        const nd = nodes[4 * e + k];
        if (nd < n) v += G[12 * e + 3 * k] * x[3 * nd] + G[12 * e + 3 * k + 1] * x[3 * nd + 1] + G[12 * e + 3 * k + 2] * x[3 * nd + 2];
      }
      phi[e] = v;
    }
    mom.fill(0);
    for (let t = 0; t < b.tri.length / 3; t++)
      for (let a = 0; a < 3; a++) {
        const ea = b.tri[3 * t + a];
        if (ea < 0) continue;
        for (let c = 0; c < 3; c++) {
          const ec = b.tri[3 * t + c];
          if (ec >= 0) mom[ea] += b.M[9 * t + 3 * a + c] * phi[ec];
        }
      }
    for (let e = 0; e < nh; e++)
      for (let k = 0; k < 4; k++) {
        const nd = nodes[4 * e + k];
        if (nd < n) for (let c = 0; c < 3; c++) y[3 * nd + c] += mom[e] * G[12 * e + 3 * k + c];
      }
  };
}

// Highest angular frequency of the linearized system at rest: power iteration
// on M^-1/2 K M^-1/2 (seeded, so the same system always gives the same value).
// The Rayleigh quotient approaches the top eigenvalue from below; the contact
// stiffness of a mallet is added on top (Weyl: lambda(A + B) <= lambda(A) + lambda(B)).
export function omegaMaxEstimate(sys: BeadSystem, mallet?: Mallet, iterations = 400): number {
  const { n, pinned } = sys;
  const K = restStiffness(sys);
  const rs = Float64Array.from(sys.mass, (m, i) => (pinned[i] ? 0 : 1 / Math.sqrt(m)));
  let seed = 12345;
  const rand = () => ((seed = (seed * 1103515245 + 12345) % 2147483648) / 2147483648) - 0.5;
  const v = new Float64Array(3 * n);
  const z = new Float64Array(3 * n);
  const y = new Float64Array(3 * n);
  for (let k = 0; k < 3 * n; k++) v[k] = rs[Math.floor(k / 3)] > 0 ? rand() : 0;
  let lam = 0;
  for (let it = 0; it < iterations; it++) {
    let nv = 0;
    for (let k = 0; k < 3 * n; k++) nv += v[k] * v[k];
    nv = Math.sqrt(nv);
    if (!(nv > 0)) break;
    for (let k = 0; k < 3 * n; k++) z[k] = (v[k] / nv) * rs[Math.floor(k / 3)];
    K(z, y);
    lam = 0;
    for (let k = 0; k < 3 * n; k++) {
      y[k] *= rs[Math.floor(k / 3)];
      lam += (y[k] * v[k]) / nv;
    }
    v.set(y);
  }
  let w2 = Math.max(lam, 0);
  if (mallet) {
    let inv = 1 / mallet.mass;
    for (let i = 0; i < n; i++) if (!pinned[i] && mallet.weights[i] > 0) inv += (mallet.weights[i] * mallet.weights[i]) / sys.mass[i];
    w2 += contactStiffness(mallet) * inv;
  }
  return Math.sqrt(w2);
}

export interface TimeStep {
  omegaMax: number;
  dt: number;
  stepsPerSample: number;
}

// dt = cfl * 2 / omega_max, rounded down to a whole number of steps per audio sample
export function timeStep(sys: BeadSystem, mallet: Mallet | undefined, fs: number, cfl: number): TimeStep {
  const omegaMax = omegaMaxEstimate(sys, mallet);
  const stepsPerSample = Math.max(1, Math.ceil(omegaMax / (fs * cfl * 2)));
  return { omegaMax, dt: 1 / (fs * stepsPerSample), stepsPerSample };
}

// bending energy at positions X, optionally only over triangles with keep(t)
export function bendingEnergy(b: Bending, X: Float64Array, keep?: (t: number) => boolean): number {
  const g = new Float64Array(13);
  const phi = Float64Array.from(b.w0, (_, e) => b.sign[e] * (hingeAngleGrad(X, b.w0[e], b.w1[e], b.e0[e], b.e1[e], g) - b.theta0[e]));
  let E = 0;
  for (let t = 0; t < b.tri.length / 3; t++) {
    if (keep && !keep(t)) continue;
    for (let a = 0; a < 3; a++)
      for (let c = 0; c < 3; c++) {
        const ea = b.tri[3 * t + a];
        const ec = b.tri[3 * t + c];
        if (ea >= 0 && ec >= 0) E += 0.5 * phi[ea] * b.M[9 * t + 3 * a + c] * phi[ec];
      }
  }
  return E;
}

export type Progress = (fraction: number) => void;
// receives the result so far (a prefix of every output) while a run is in progress
export type PartialResult = (r: SimResult) => void;

// wall time before the first partial result, and between later ones (ms)
export const PARTIAL_FIRST_MS = 1000;
export const PARTIAL_EVERY_MS = 2000;

// Velocity Verlet (kick-drift-kick), explicit. Spring and hinge forces are
// exact (geometric nonlinearity kept); damping and friction use the half-step
// velocity in the kick, and the full-step velocity whenever an audio sample is
// recorded.
export function simulate(sys: BeadSystem, exc: Excitation, opt: SimOptions, progress?: Progress, partial?: PartialResult): SimResult {
  const t0 = performance.now();
  const { n, si, sj, sk, sL0, sc, mass, pinned, anchors } = sys;
  const nn = n + anchors.length;
  const ns = si.length;
  const bend = sys.bending;
  const nh = bend ? bend.w0.length : 0;
  const nt = bend ? bend.tri.length / 3 : 0;
  const hG = new Float64Array(12 * nh); // gradients of phi
  const hPhi = new Float64Array(nh);
  const hRate = new Float64Array(nh);
  const hMom = new Float64Array(nh);
  const hMomEl = new Float64Array(nh); // elastic part of the moments

  const mallet = exc.mallet;
  const { omegaMax, dt, stepsPerSample } = timeStep(sys, mallet, opt.fs, opt.cfl);
  const nSamples = Math.round(opt.duration * opt.fs);
  const samplesPerFrame = Math.max(1, Math.round(opt.fs / opt.frameRate));
  const nFrames = Math.floor((nSamples - 1) / samplesPerFrame) + 1;

  const X = new Float64Array(3 * nn);
  const V = new Float64Array(3 * nn);
  const F = new Float64Array(3 * n); // force on beads
  const Fel = new Float64Array(3 * n); // elastic (position-dependent) part
  const A = new Float64Array(3 * n);
  for (let i = 0; i < n; i++) {
    for (let c = 0; c < 3; c++) {
      const k = 3 * i + c;
      X[k] = sys.rest[k] + (pinned[i] ? 0 : exc.u0[k]);
      V[k] = pinned[i] ? 0 : exc.v0[k];
    }
  }
  const slide: number[] = [];
  for (let a = 0; a < anchors.length; a++) {
    for (let c = 0; c < 3; c++) X[3 * (n + a) + c] = anchors[a].pos[c];
    if (anchors[a].slideBead >= 0) slide.push(a);
  }
  const syncAnchors = () => {
    for (const a of slide) {
      const b = anchors[a].slideBead;
      const k = 3 * (n + a);
      X[k + 1] = X[3 * b + 1];
      X[k + 2] = X[3 * b + 2];
      V[k + 1] = V[3 * b + 1];
      V[k + 2] = V[3 * b + 2];
    }
  };
  syncAnchors();

  const force = exc.force;
  const fShape = force?.shape;
  const pulse = (t: number) =>
    force && t >= 0 && t <= force.tau ? force.F0 * 0.5 * (1 - Math.cos((2 * Math.PI * t) / force.tau)) : 0;

  // mallet state: 3D position q, velocity qv, acceleration qa
  const q = new Float64Array(3);
  const qv = new Float64Array(mallet ? mallet.vel0 : [0, 0, 0]);
  const qa = new Float64Array(3);
  let fContact = 0;
  let eContact = 0;
  let pFriction = 0;
  const patch: number[] = [];
  let mEff = 1;
  if (mallet) {
    let inv = 0;
    for (let i = 0; i < n; i++)
      if (mallet.weights[i] > 0) {
        patch.push(i);
        inv += (mallet.weights[i] * mallet.weights[i]) / mass[i];
      }
    mEff = 1 / (1 / mallet.mass + inv);
  }
  const [nx, ny, nz] = mallet?.normal ?? [0, 0, 0];
  const contact = () => {
    fContact = 0;
    eContact = 0;
    pFriction = 0;
    qa.fill(0);
    if (!mallet) return;
    const w = mallet.weights;
    let ux = 0, uy = 0, uz = 0, vx = 0, vy = 0, vz = 0;
    for (const i of patch) {
      const k = 3 * i;
      ux += w[i] * (X[k] - sys.rest[k]);
      uy += w[i] * (X[k + 1] - sys.rest[k + 1]);
      uz += w[i] * (X[k + 2] - sys.rest[k + 2]);
      vx += w[i] * V[k];
      vy += w[i] * V[k + 1];
      vz += w[i] * V[k + 2];
    }
    const d = (q[0] - ux) * nx + (q[1] - uy) * ny + (q[2] - uz) * nz;
    if (!(d > 0)) return;
    if (mallet.reach !== undefined) {
      const ox = q[0] - ux - d * nx, oy = q[1] - uy - d * ny, oz = q[2] - uz - d * nz;
      if (ox * ox + oy * oy + oz * oz > mallet.reach * mallet.reach) return;
    }
    const fn = mallet.K * d ** mallet.p;
    fContact = fn;
    eContact = (fn * d) / (mallet.p + 1);
    // force on the mallet: -fn normal + friction against tangential slip
    let fx = -fn * nx;
    let fy = -fn * ny;
    let fz = -fn * nz;
    if (mallet.mu > 0) {
      let rx = qv[0] - vx, ry = qv[1] - vy, rz = qv[2] - vz;
      const rn = rx * nx + ry * ny + rz * nz;
      rx -= rn * nx;
      ry -= rn * ny;
      rz -= rn * nz;
      const slip = Math.hypot(rx, ry, rz);
      // regularized Coulomb; the slip scale keeps the explicit step stable
      const eps = Math.max(1e-3, (4 * mallet.mu * fn * dt) / mEff);
      const g = (mallet.mu * fn) / Math.sqrt(slip * slip + eps * eps);
      fx -= g * rx;
      fy -= g * ry;
      fz -= g * rz;
      pFriction = g * slip * slip;
    }
    qa[0] = fx / mallet.mass;
    qa[1] = fy / mallet.mass;
    qa[2] = fz / mallet.mass;
    for (const i of patch) {
      F[3 * i] -= w[i] * fx;
      F[3 * i + 1] -= w[i] * fy;
      F[3 * i + 2] -= w[i] * fz;
    }
  };

  const hOut = new Float64Array(13);
  let eHinge = 0;
  const anyDashpot = sc.some((c) => c !== 0);
  // Position-dependent forces: spring and bending elasticity into Fel; hinge
  // angles, gradients and elastic moments; rigid-body geometry.
  const elasticForces = () => {
    Fel.fill(0);
    for (let s = 0; s < ns; s++) {
      const a = 3 * si[s];
      const b = 3 * sj[s];
      const dx = X[b] - X[a];
      const dy = X[b + 1] - X[a + 1];
      const dz = X[b + 2] - X[a + 2];
      const l = Math.sqrt(dx * dx + dy * dy + dz * dz);
      const g = (sk[s] * (l - sL0[s])) / l;
      if (si[s] < n) {
        Fel[a] += g * dx;
        Fel[a + 1] += g * dy;
        Fel[a + 2] += g * dz;
      }
      if (sj[s] < n) {
        Fel[b] -= g * dx;
        Fel[b + 1] -= g * dy;
        Fel[b + 2] -= g * dz;
      }
    }
    eHinge = 0;
    if (bend) {
      for (let e = 0; e < nh; e++) {
        const th = hingeAngleGrad(X, bend.w0[e], bend.w1[e], bend.e0[e], bend.e1[e], hOut);
        const sg = bend.sign[e];
        hPhi[e] = sg * (th - bend.theta0[e]);
        const o = 12 * e;
        for (let q = 0; q < 12; q++) hG[o + q] = sg * hOut[1 + q];
      }
      hMomEl.fill(0);
      const tri = bend.tri;
      const Mt = bend.M;
      for (let t = 0; t < nt; t++) {
        for (let a = 0; a < 3; a++) {
          const ea = tri[3 * t + a];
          if (ea < 0) continue;
          let m = 0;
          for (let c = 0; c < 3; c++) {
            const ec = tri[3 * t + c];
            if (ec >= 0) m += Mt[9 * t + 3 * a + c] * hPhi[ec];
          }
          eHinge += 0.5 * hPhi[ea] * m;
          hMomEl[ea] += m;
        }
      }
    }
    if (sys.freeRigid) rigidGeometry();
  };

  // Fills F with total force on beads and A with accelerations, from the
  // elastic part (elasticForces must have run at this position) and the
  // velocity-dependent part at V; sets damping power (W, >= 0) and external
  // power (W), both at V.
  let pDamp = 0;
  let pExt = 0;
  const dynamicForces = (t: number) => {
    F.set(Fel);
    contact();
    pDamp = pFriction;
    if (anyDashpot)
      for (let s = 0; s < ns; s++) {
        const c = sc[s];
        if (c === 0) continue;
        const a = 3 * si[s];
        const b = 3 * sj[s];
        const rx = V[b] - V[a];
        const ry = V[b + 1] - V[a + 1];
        const rz = V[b + 2] - V[a + 2];
        pDamp += c * (rx * rx + ry * ry + rz * rz);
        if (si[s] < n) {
          F[a] += c * rx;
          F[a + 1] += c * ry;
          F[a + 2] += c * rz;
        }
        if (sj[s] < n) {
          F[b] -= c * rx;
          F[b + 1] -= c * ry;
          F[b + 2] -= c * rz;
        }
      }
    if (bend) {
      const s1 = bend.sigma1;
      hMom.set(hMomEl);
      if (s1 !== 0) {
        for (let e = 0; e < nh; e++) {
          const o = 12 * e;
          const w0 = 3 * bend.w0[e], w1 = 3 * bend.w1[e], e0 = 3 * bend.e0[e], e1 = 3 * bend.e1[e];
          hRate[e] =
            hG[o] * V[w0] + hG[o + 1] * V[w0 + 1] + hG[o + 2] * V[w0 + 2] +
            hG[o + 3] * V[w1] + hG[o + 4] * V[w1 + 1] + hG[o + 5] * V[w1 + 2] +
            hG[o + 6] * V[e0] + hG[o + 7] * V[e0 + 1] + hG[o + 8] * V[e0 + 2] +
            hG[o + 9] * V[e1] + hG[o + 10] * V[e1 + 1] + hG[o + 11] * V[e1 + 2];
        }
        const tri = bend.tri;
        const Mt = bend.M;
        for (let t = 0; t < nt; t++) {
          for (let a = 0; a < 3; a++) {
            const ea = tri[3 * t + a];
            if (ea < 0) continue;
            let r = 0;
            for (let c = 0; c < 3; c++) {
              const ec = tri[3 * t + c];
              if (ec >= 0) r += Mt[9 * t + 3 * a + c] * hRate[ec];
            }
            pDamp += 2 * s1 * hRate[ea] * r;
            hMom[ea] += 2 * s1 * r;
          }
        }
      }
      for (let e = 0; e < nh; e++) {
        const m = hMom[e];
        if (m === 0) continue;
        const o = 12 * e;
        const w0 = 3 * bend.w0[e], w1 = 3 * bend.w1[e], e0 = 3 * bend.e0[e], e1 = 3 * bend.e1[e];
        F[w0] -= m * hG[o];
        F[w0 + 1] -= m * hG[o + 1];
        F[w0 + 2] -= m * hG[o + 2];
        F[w1] -= m * hG[o + 3];
        F[w1 + 1] -= m * hG[o + 4];
        F[w1 + 2] -= m * hG[o + 5];
        F[e0] -= m * hG[o + 6];
        F[e0 + 1] -= m * hG[o + 7];
        F[e0 + 2] -= m * hG[o + 8];
        F[e1] -= m * hG[o + 9];
        F[e1 + 1] -= m * hG[o + 10];
        F[e1 + 2] -= m * hG[o + 11];
      }
    }
    const s0 = sys.sigma0;
    const p = fShape ? pulse(t) : 0;
    pExt = 0;
    for (let i = 0; i < n; i++) {
      const k = 3 * i;
      if (pinned[i]) {
        F[k] = F[k + 1] = F[k + 2] = 0;
        continue;
      }
      if (s0 !== 0) {
        const m2 = 2 * s0 * mass[i];
        for (let c = 0; c < 3; c++) {
          F[k + c] -= m2 * V[k + c];
          pDamp += m2 * V[k + c] * V[k + c];
        }
      }
      if (p !== 0) {
        for (let c = 0; c < 3; c++) {
          const fe = p * fShape![k + c];
          F[k + c] += fe;
          pExt += fe * V[k + c];
        }
      }
    }
    for (let i = 0; i < n; i++) {
      const im = 1 / mass[i];
      A[3 * i] = F[3 * i] * im;
      A[3 * i + 1] = F[3 * i + 1] * im;
      A[3 * i + 2] = F[3 * i + 2] * im;
    }
    if (sys.freeRigid) removeRigid(A);
  };
  const computeForces = (t: number) => {
    elasticForces();
    dynamicForces(t);
  };

  // centre of mass and inverse inertia tensor at the current positions
  let totalMass = 0;
  for (let i = 0; i < n; i++) totalMass += mass[i];
  const cm = new Float64Array(3);
  const Iinv = new Float64Array(6);
  const rigidGeometry = () => {
    let cx = 0, cy = 0, cz = 0;
    for (let i = 0; i < n; i++) {
      cx += mass[i] * X[3 * i];
      cy += mass[i] * X[3 * i + 1];
      cz += mass[i] * X[3 * i + 2];
    }
    cx /= totalMass;
    cy /= totalMass;
    cz /= totalMass;
    let ixx = 0, iyy = 0, izz = 0, ixy = 0, ixz = 0, iyz = 0;
    for (let i = 0; i < n; i++) {
      const m = mass[i];
      const rx = X[3 * i] - cx, ry = X[3 * i + 1] - cy, rz = X[3 * i + 2] - cz;
      ixx += m * (ry * ry + rz * rz);
      iyy += m * (rx * rx + rz * rz);
      izz += m * (rx * rx + ry * ry);
      ixy -= m * rx * ry;
      ixz -= m * rx * rz;
      iyz -= m * ry * rz;
    }
    cm[0] = cx;
    cm[1] = cy;
    cm[2] = cz;
    invert3(ixx, ixy, ixz, iyy, iyz, izz, Iinv);
  };

  // subtract the rigid-body part (mass-weighted) of a per-bead field u, about
  // the geometry from rigidGeometry; returns the kinetic energy of that part
  // when u is a velocity
  const removeRigid = (u: Float64Array): number => {
    const [cx, cy, cz] = cm;
    let px = 0, py = 0, pz = 0, lx = 0, ly = 0, lz = 0;
    for (let i = 0; i < n; i++) {
      const m = mass[i];
      const rx = X[3 * i] - cx, ry = X[3 * i + 1] - cy, rz = X[3 * i + 2] - cz;
      const ux = u[3 * i], uy = u[3 * i + 1], uz = u[3 * i + 2];
      px += m * ux;
      py += m * uy;
      pz += m * uz;
      lx += m * (ry * uz - rz * uy);
      ly += m * (rz * ux - rx * uz);
      lz += m * (rx * uy - ry * ux);
    }
    const [a, b, c, d, e, f] = Iinv;
    const wx = a * lx + b * ly + c * lz;
    const wy = b * lx + d * ly + e * lz;
    const wz = c * lx + e * ly + f * lz;
    const vx = px / totalMass, vy = py / totalMass, vz = pz / totalMass;
    for (let i = 0; i < n; i++) {
      const rx = X[3 * i] - cx, ry = X[3 * i + 1] - cy, rz = X[3 * i + 2] - cz;
      u[3 * i] -= vx + (wy * rz - wz * ry);
      u[3 * i + 1] -= vy + (wz * rx - wx * rz);
      u[3 * i + 2] -= vz + (wx * ry - wy * rx);
    }
    return 0.5 * (px * vx + py * vy + pz * vz) + 0.5 * (wx * lx + wy * ly + wz * lz);
  };
  if (sys.freeRigid) {
    rigidGeometry();
    removeRigid(V); // an initial strike carries no rigid motion either
  }

  const velocityForces = sys.sigma0 !== 0 || anyDashpot || (bend?.sigma1 ?? 0) !== 0 || (mallet?.mu ?? 0) > 0;
  let maxStretch = 0;
  const malletEnergy = () => (mallet ? 0.5 * mallet.mass * (qv[0] ** 2 + qv[1] ** 2 + qv[2] ** 2) + eContact : 0);
  // total mechanical energy at the current state; computeForces must have run
  // at this state (it sets the hinge and contact energies)
  const energyNow = () => {
    let e = 0;
    for (let i = 0; i < n; i++) {
      const k = 3 * i;
      e += 0.5 * mass[i] * (V[k] * V[k] + V[k + 1] * V[k + 1] + V[k + 2] * V[k + 2]);
    }
    for (let s = 0; s < ns; s++) {
      const a = 3 * si[s];
      const b = 3 * sj[s];
      const l = Math.hypot(X[b] - X[a], X[b + 1] - X[a + 1], X[b + 2] - X[a + 2]);
      const d = l - sL0[s];
      e += 0.5 * sk[s] * d * d;
      if (sj[s] < n && si[s] < n) maxStretch = Math.max(maxStretch, Math.abs(d) / sL0[s]);
    }
    return e + eHinge + malletEnergy();
  };

  const audio = new Float32Array(nSamples);
  const frames = new Float32Array(nFrames * 3 * n);
  const frameTimes = new Float64Array(nFrames);
  const energy = new Float64Array(nFrames);
  const balance = new Float64Array(nFrames);
  const contactForce = mallet ? new Float32Array(nSamples) : undefined;
  const malletPos = mallet ? new Float32Array(3 * nFrames) : undefined;
  const w = opt.pickup;
  const delay = opt.pickupDelay;

  // energy is reported relative to the equilibrium state (pre-tension energy subtracted)
  let Eeq = 0;
  for (let s = 0; s < ns; s++) {
    const pa = (node: number, c: number) => (node < n ? sys.rest[3 * node + c] : anchors[node - n].pos[c]);
    const l = Math.hypot(pa(sj[s], 0) - pa(si[s], 0), pa(sj[s], 1) - pa(si[s], 1), pa(sj[s], 2) - pa(si[s], 2));
    Eeq += 0.5 * sk[s] * (l - sL0[s]) ** 2;
  }
  computeForces(0);
  const E0 = energyNow();
  let dissipated = 0;
  let workIn = 0;
  let t = 0;
  let frame = 0;
  const progressEvery = Math.max(1, Math.floor(nSamples / 100));
  let nextPartial = t0 + PARTIAL_FIRST_MS;
  // The outputs up to sample count S and frame count F: copies, except for the
  // finished run (a partial result may be transferred away). With a delayed
  // pickup, audio[i] only gets terms from samples <= i, so the first S are complete.
  const result = (S: number, F: number): SimResult => {
    const all = S === nSamples;
    return {
      audio: all ? audio : audio.slice(0, S),
      fs: opt.fs,
      frames: all ? frames : frames.slice(0, F * 3 * n),
      frameTimes: all ? frameTimes : frameTimes.slice(0, F),
      energy: all ? energy : energy.slice(0, F),
      balance: all ? balance : balance.slice(0, F),
      dt,
      stepsPerSample,
      omegaMax,
      wallSeconds: (performance.now() - t0) / 1000,
      maxStretch,
      contactForce: all ? contactForce : contactForce?.slice(0, S),
      malletPos: all ? malletPos : malletPos?.slice(0, 3 * F),
    };
  };

  for (let smp = 0; smp < nSamples; smp++) {
    // record at t = smp / fs; A holds forces at full-step velocity here
    if (delay) {
      for (let i = 0; i < n; i++) {
        const k = 3 * i;
        const c = w[k] * A[k] + w[k + 1] * A[k + 1] + w[k + 2] * A[k + 2];
        if (c === 0) continue;
        const pos = smp + delay[i];
        const i0 = Math.floor(pos);
        const fr = pos - i0;
        if (i0 < nSamples) audio[i0] += (1 - fr) * c;
        if (i0 + 1 < nSamples) audio[i0 + 1] += fr * c;
      }
    } else {
      let acc = 0;
      for (let k = 0; k < 3 * n; k++) acc += w[k] * A[k];
      audio[smp] = acc;
    }
    if (contactForce) contactForce[smp] = fContact;
    if (smp % samplesPerFrame === 0 && frame < nFrames) {
      if (malletPos) malletPos.set(q, 3 * frame);
      const off = frame * 3 * n;
      for (let k = 0; k < 3 * n; k++) frames[off + k] = X[k] - sys.rest[k];
      frameTimes[frame] = t;
      const eTotal = energyNow();
      // the plot shows the instrument's own energy; the balance includes the mallet
      energy[frame] = eTotal - Eeq - malletEnergy();
      balance[frame] = eTotal + dissipated - workIn - E0;
      frame++;
    }
    if (smp === nSamples - 1) break;
    if (partial && smp % 32 === 0 && performance.now() >= nextPartial) {
      partial(result(smp + 1, frame));
      nextPartial = performance.now() + PARTIAL_EVERY_MS;
    }

    for (let st = 0; st < stepsPerSample; st++) {
      // kick (half), drift
      for (let k = 0; k < 3 * n; k++) {
        V[k] += 0.5 * dt * A[k];
        X[k] += dt * V[k];
      }
      if (mallet)
        for (let c = 0; c < 3; c++) {
          qv[c] += 0.5 * dt * qa[c];
          q[c] += dt * qv[c];
        }
      syncAnchors();
      t += dt;
      // forces at new position with half-step velocity; powers use it too
      computeForces(t);
      dissipated += dt * pDamp;
      workIn += dt * pExt;
      for (let k = 0; k < 3 * n; k++) V[k] += 0.5 * dt * A[k];
      if (mallet) for (let c = 0; c < 3; c++) qv[c] += 0.5 * dt * qa[c];
      // the rigid-free subspace turns with the body; drop what drifts out of it
      if (sys.freeRigid) dissipated += removeRigid(V);
      syncAnchors();
    }
    // at a sample time: recompute with the full-step velocity (positions unchanged)
    if (velocityForces) dynamicForces(t);
    if (progress && smp % progressEvery === 0) progress(smp / nSamples);
  }

  return result(nSamples, nFrames);
}

// inverse of the symmetric 3x3 matrix [a b c; b d e; c e f], as its six
// distinct entries (zero if singular)
function invert3(a: number, b: number, c: number, d: number, e: number, f: number, out: Float64Array) {
  const A00 = d * f - e * e, A01 = c * e - b * f, A02 = b * e - c * d;
  const A11 = a * f - c * c, A12 = b * c - a * e, A22 = a * d - b * b;
  const det = a * A00 + b * A01 + c * A02;
  if (!(Math.abs(det) > 0)) {
    out.fill(0);
    return;
  }
  out[0] = A00 / det;
  out[1] = A01 / det;
  out[2] = A02 / det;
  out[3] = A11 / det;
  out[4] = A12 / det;
  out[5] = A22 / det;
}