// Linear transverse modes of a flat bead system.
//
// When every bead rests in a plane z = const and every spring lies in that
// plane, small transverse (z) motion decouples from in-plane motion, and its
// stiffness matrix is a graph Laplacian weighted by spring tension over length:
// K = sum_s (f_s / L_s) (e_i - e_j)(e_i - e_j)^T. We return the eigenpairs of
// M^{-1/2} K M^{-1/2} over the free beads.

import type { BeadSystem, Modes } from "./types";

export type SymEig = (A: Float64Array, n: number) => { values: Float64Array; vectors: Float64Array };

export function isPlanar(sys: BeadSystem): boolean {
  const { n, pos0 } = sys;
  let zmin = Infinity;
  let zmax = -Infinity;
  let span = 0;
  for (let i = 0; i < n; i++) {
    zmin = Math.min(zmin, pos0[3 * i + 2]);
    zmax = Math.max(zmax, pos0[3 * i + 2]);
    span = Math.max(span, Math.abs(pos0[3 * i]), Math.abs(pos0[3 * i + 1]));
  }
  return zmax - zmin <= 1e-12 * Math.max(span, 1e-30);
}

export function transverseStiffness(sys: BeadSystem): { K: Float64Array; free: Int32Array } {
  const { n, pinned, si, sj, sk, sL0, sk3, pos0, mass } = sys;
  const freeIdx = new Int32Array(n).fill(-1);
  const free: number[] = [];
  for (let i = 0; i < n; i++) if (!pinned[i]) freeIdx[i] = free.push(i) - 1;
  const nf = free.length;
  const K = new Float64Array(nf * nf);
  for (let s = 0; s < si.length; s++) {
    const i3 = 3 * si[s];
    const j3 = 3 * sj[s];
    const L = Math.hypot(pos0[j3] - pos0[i3], pos0[j3 + 1] - pos0[i3 + 1], pos0[j3 + 2] - pos0[i3 + 2]);
    const e = L - sL0[s];
    const g = (sk[s] * e + sk3[s] * e * e * e) / L;
    const a = freeIdx[si[s]];
    const b = freeIdx[sj[s]];
    if (a >= 0) K[a * nf + a] += g;
    if (b >= 0) K[b * nf + b] += g;
    if (a >= 0 && b >= 0) {
      K[a * nf + b] -= g;
      K[b * nf + a] -= g;
    }
  }
  for (let a = 0; a < nf; a++) {
    for (let b = 0; b < nf; b++) K[a * nf + b] /= Math.sqrt(mass[free[a]] * mass[free[b]]);
  }
  return { K, free: Int32Array.from(free) };
}

export function computeModes(sys: BeadSystem, eig: SymEig): Modes {
  const { K, free } = transverseStiffness(sys);
  const nf = free.length;
  const { values, vectors } = eig(K, nf);
  const order = Array.from(values.keys()).sort((p, q) => values[p] - values[q]);
  const freqs = new Float64Array(nf);
  const psi = new Float64Array(nf * nf);
  order.forEach((src, k) => {
    freqs[k] = Math.sqrt(Math.max(values[src], 0)) / (2 * Math.PI);
    psi.set(vectors.subarray(src * nf, (src + 1) * nf), k * nf);
  });
  return { freqs, psi, free };
}