1// Linear transverse modes of a flat bead system.
2//
3// When every bead rests in a plane z = const and every spring lies in that
4// plane, small transverse (z) motion decouples from in-plane motion, and its
5// stiffness matrix is a graph Laplacian weighted by spring tension over length:
6// K = sum_s (f_s / L_s) (e_i - e_j)(e_i - e_j)^T. We return the eigenpairs of
7// M^{-1/2} K M^{-1/2} over the free beads.
9import type { BeadSystem, Modes } from "./types";
11export type SymEig = (A: Float64Array, n: number) => { values: Float64Array; vectors: Float64Array };
13export function isPlanar(sys: BeadSystem): boolean {
14 const { n, pos0 } = sys;
15 let zmin = Infinity;
16 let zmax = -Infinity;
17 let span = 0;
18 for (let i = 0; i < n; i++) {
19 zmin = Math.min(zmin, pos0[3 * i + 2]);
20 zmax = Math.max(zmax, pos0[3 * i + 2]);
21 span = Math.max(span, Math.abs(pos0[3 * i]), Math.abs(pos0[3 * i + 1]));
22 }
23 return zmax - zmin <= 1e-12 * Math.max(span, 1e-30);
24}
26export function transverseStiffness(sys: BeadSystem): { K: Float64Array; free: Int32Array } {
27 const { n, pinned, si, sj, sk, sL0, sk3, pos0, mass } = sys;
28 const freeIdx = new Int32Array(n).fill(-1);
29 const free: number[] = [];
30 for (let i = 0; i < n; i++) if (!pinned[i]) freeIdx[i] = free.push(i) - 1;
31 const nf = free.length;
32 const K = new Float64Array(nf * nf);
33 for (let s = 0; s < si.length; s++) {
34 const i3 = 3 * si[s];
35 const j3 = 3 * sj[s];
36 const L = Math.hypot(pos0[j3] - pos0[i3], pos0[j3 + 1] - pos0[i3 + 1], pos0[j3 + 2] - pos0[i3 + 2]);
37 const e = L - sL0[s];
38 const g = (sk[s] * e + sk3[s] * e * e * e) / L;
39 const a = freeIdx[si[s]];
40 const b = freeIdx[sj[s]];
41 if (a >= 0) K[a * nf + a] += g;
42 if (b >= 0) K[b * nf + b] += g;
43 if (a >= 0 && b >= 0) {
44 K[a * nf + b] -= g;
45 K[b * nf + a] -= g;
46 }
47 }
48 for (let a = 0; a < nf; a++) {
49 for (let b = 0; b < nf; b++) K[a * nf + b] /= Math.sqrt(mass[free[a]] * mass[free[b]]);
50 }
51 return { K, free: Int32Array.from(free) };
52}
54export function computeModes(sys: BeadSystem, eig: SymEig): Modes {
55 const { K, free } = transverseStiffness(sys);
56 const nf = free.length;
57 const { values, vectors } = eig(K, nf);
58 const order = Array.from(values.keys()).sort((p, q) => values[p] - values[q]);
59 const freqs = new Float64Array(nf);
60 const psi = new Float64Array(nf * nf);
61 order.forEach((src, k) => {
62 freqs[k] = Math.sqrt(Math.max(values[src], 0)) / (2 * Math.PI);
63 psi.set(vectors.subarray(src * nf, (src + 1) * nf), k * nf);
64 });
65 return { freqs, psi, free };
66}