1import type { BeadSystem, Excitation, SimOptions, SimResult } from './types';
3// Upper bound on the highest angular frequency of the linearized system
4// (Gershgorin on M^-1 K at rest). A spring under tension T has stiffness k along
5// its axis and T/l across it, so each spring contributes 2 max(k, |T|/l) / m.
6export function omegaMaxBound(sys: BeadSystem): number {
7 const { n, si, sj, sk, sL0 } = sys;
8 const rowSum = new Float64Array(n);
9 const p = (node: number, c: number) =>
10 node < n ? sys.rest[3 * node + c] : sys.anchors[node - n].pos[c];
11 for (let s = 0; s < si.length; s++) {
12 const a = si[s];
13 const b = sj[s];
14 let l2 = 0;
15 for (let c = 0; c < 3; c++) l2 += (p(b, c) - p(a, c)) ** 2;
16 const l = Math.sqrt(l2);
17 const T = sk[s] * (l - sL0[s]);
18 const kmax = Math.max(sk[s], Math.abs(T) / l);
19 if (a < n) rowSum[a] += 2 * kmax;
20 if (b < n) rowSum[b] += 2 * kmax;
21 }
22 let w2 = 0;
23 for (let i = 0; i < n; i++) if (!sys.pinned[i]) w2 = Math.max(w2, rowSum[i] / sys.mass[i]);
24 return Math.sqrt(w2);
25}
27export type Progress = (fraction: number) => void;
29// Velocity Verlet (kick-drift-kick), explicit. Spring forces are exact
30// (geometric nonlinearity kept); damping uses the half-step velocity in the
31// kick, and the full-step velocity whenever an audio sample is recorded.
32export function simulate(
33 sys: BeadSystem,
34 exc: Excitation,
35 opt: SimOptions,
36 progress?: Progress,
37): SimResult {
38 const t0 = performance.now();
39 const { n, si, sj, sk, sL0, sc, mass, pinned, anchors } = sys;
40 const nn = n + anchors.length;
41 const ns = si.length;
43 const omegaMax = omegaMaxBound(sys);
44 const dtMax = (opt.cfl * 2) / omegaMax;
45 const stepsPerSample = Math.max(1, Math.ceil(1 / (opt.fs * dtMax)));
46 const dt = 1 / (opt.fs * stepsPerSample);
47 const nSamples = Math.round(opt.duration * opt.fs);
48 const samplesPerFrame = Math.max(1, Math.round(opt.fs / opt.frameRate));
49 const nFrames = Math.floor((nSamples - 1) / samplesPerFrame) + 1;
51 const X = new Float64Array(3 * nn);
52 const V = new Float64Array(3 * nn);
53 const F = new Float64Array(3 * n); // force on beads
54 const A = new Float64Array(3 * n);
55 for (let i = 0; i < n; i++) {
56 for (let c = 0; c < 3; c++) {
57 const k = 3 * i + c;
58 X[k] = sys.rest[k] + (pinned[i] ? 0 : exc.u0[k]);
59 V[k] = pinned[i] ? 0 : exc.v0[k];
60 }
61 }
62 const slide: number[] = [];
63 for (let a = 0; a < anchors.length; a++) {
64 for (let c = 0; c < 3; c++) X[3 * (n + a) + c] = anchors[a].pos[c];
65 if (anchors[a].slideBead >= 0) slide.push(a);
66 }
67 const syncAnchors = () => {
68 for (const a of slide) {
69 const b = anchors[a].slideBead;
70 const k = 3 * (n + a);
71 X[k + 1] = X[3 * b + 1];
72 X[k + 2] = X[3 * b + 2];
73 V[k + 1] = V[3 * b + 1];
74 V[k + 2] = V[3 * b + 2];
75 }
76 };
77 syncAnchors();
79 const force = exc.force;
80 const fShape = force?.shape;
81 const pulse = (t: number) =>
82 force && t >= 0 && t <= force.tau
83 ? force.F0 * 0.5 * (1 - Math.cos((2 * Math.PI * t) / force.tau))
84 : 0;
86 // Fills F with total force on beads; returns damping power (W, >= 0)
87 // and external power (W), both at velocity V.
88 let pDamp = 0;
89 let pExt = 0;
90 const computeForces = (t: number) => {
91 F.fill(0);
92 pDamp = 0;
93 for (let s = 0; s < ns; s++) {
94 const a = 3 * si[s];
95 const b = 3 * sj[s];
96 const dx = X[b] - X[a];
97 const dy = X[b + 1] - X[a + 1];
98 const dz = X[b + 2] - X[a + 2];
99 const l = Math.sqrt(dx * dx + dy * dy + dz * dz);
100 const g = (sk[s] * (l - sL0[s])) / l;
101 let fx = g * dx;
102 let fy = g * dy;
103 let fz = g * dz;
104 const c = sc[s];
105 if (c !== 0) {
106 const rx = V[b] - V[a];
107 const ry = V[b + 1] - V[a + 1];
108 const rz = V[b + 2] - V[a + 2];
109 fx += c * rx;
110 fy += c * ry;
111 fz += c * rz;
112 pDamp += c * (rx * rx + ry * ry + rz * rz);
113 }
114 if (si[s] < n) {
115 F[a] += fx;
116 F[a + 1] += fy;
117 F[a + 2] += fz;
118 }
119 if (sj[s] < n) {
120 F[b] -= fx;
121 F[b + 1] -= fy;
122 F[b + 2] -= fz;
123 }
124 }
125 const s0 = sys.sigma0;
126 const p = fShape ? pulse(t) : 0;
127 pExt = 0;
128 for (let i = 0; i < n; i++) {
129 const k = 3 * i;
130 if (pinned[i]) {
131 F[k] = F[k + 1] = F[k + 2] = 0;
132 continue;
133 }
134 if (s0 !== 0) {
135 const m2 = 2 * s0 * mass[i];
136 for (let c = 0; c < 3; c++) {
137 F[k + c] -= m2 * V[k + c];
138 pDamp += m2 * V[k + c] * V[k + c];
139 }
140 }
141 if (p !== 0) {
142 for (let c = 0; c < 3; c++) {
143 const fe = p * fShape![k + c];
144 F[k + c] += fe;
145 pExt += fe * V[k + c];
146 }
147 }
148 }
149 for (let i = 0; i < n; i++) {
150 const im = 1 / mass[i];
151 A[3 * i] = F[3 * i] * im;
152 A[3 * i + 1] = F[3 * i + 1] * im;
153 A[3 * i + 2] = F[3 * i + 2] * im;
154 }
155 };
157 const damped = sys.sigma0 !== 0 || sc.some((c) => c !== 0);
158 let maxStretch = 0;
159 const energyNow = () => {
160 let e = 0;
161 for (let i = 0; i < n; i++) {
162 const k = 3 * i;
163 e += 0.5 * mass[i] * (V[k] * V[k] + V[k + 1] * V[k + 1] + V[k + 2] * V[k + 2]);
164 }
165 for (let s = 0; s < ns; s++) {
166 const a = 3 * si[s];
167 const b = 3 * sj[s];
168 const l = Math.hypot(X[b] - X[a], X[b + 1] - X[a + 1], X[b + 2] - X[a + 2]);
169 const d = l - sL0[s];
170 e += 0.5 * sk[s] * d * d;
171 maxStretch = Math.max(maxStretch, Math.abs(d) / sL0[s]);
172 }
173 return e;
174 };
176 const audio = new Float32Array(nSamples);
177 const frames = new Float32Array(nFrames * 3 * n);
178 const frameTimes = new Float64Array(nFrames);
179 const energy = new Float64Array(nFrames);
180 const balance = new Float64Array(nFrames);
181 const w = opt.pickup;
183 // energy is reported relative to the equilibrium state (pre-tension energy subtracted)
184 let Eeq = 0;
185 for (let s = 0; s < ns; s++) {
186 const pa = (node: number, c: number) => (node < n ? sys.rest[3 * node + c] : anchors[node - n].pos[c]);
187 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));
188 Eeq += 0.5 * sk[s] * (l - sL0[s]) ** 2;
189 }
190 const E0 = energyNow();
191 let dissipated = 0;
192 let workIn = 0;
193 let t = 0;
194 computeForces(0);
195 let frame = 0;
196 const progressEvery = Math.max(1, Math.floor(nSamples / 100));
198 for (let smp = 0; smp < nSamples; smp++) {
199 // record at t = smp / fs; A holds forces at full-step velocity here
200 let acc = 0;
201 for (let k = 0; k < 3 * n; k++) acc += w[k] * A[k];
202 audio[smp] = acc;
203 if (smp % samplesPerFrame === 0 && frame < nFrames) {
204 const off = frame * 3 * n;
205 for (let k = 0; k < 3 * n; k++) frames[off + k] = X[k] - sys.rest[k];
206 frameTimes[frame] = t;
207 energy[frame] = energyNow() - Eeq;
208 balance[frame] = energy[frame] + Eeq + dissipated - workIn - E0;
209 frame++;
210 }
211 if (smp === nSamples - 1) break;
213 for (let st = 0; st < stepsPerSample; st++) {
214 // kick (half), drift
215 for (let k = 0; k < 3 * n; k++) {
216 V[k] += 0.5 * dt * A[k];
217 X[k] += dt * V[k];
218 }
219 syncAnchors();
220 t += dt;
221 // forces at new position with half-step velocity; powers use it too
222 computeForces(t);
223 dissipated += dt * pDamp;
224 workIn += dt * pExt;
225 for (let k = 0; k < 3 * n; k++) V[k] += 0.5 * dt * A[k];
226 syncAnchors();
227 }
228 // at a sample time: recompute with the full-step velocity
229 if (damped) computeForces(t);
230 if (progress && smp % progressEvery === 0) progress(smp / nSamples);
231 }
233 return {
234 audio,
235 fs: opt.fs,
236 frames,
237 frameTimes,
238 energy,
239 balance,
240 dt,
241 stepsPerSample,
242 omegaMax,
243 wallSeconds: (performance.now() - t0) / 1000,
244 maxStretch,
245 };
246}