5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 1// Modal synthesis of a strike. Each mode m of the linearized plate is a
2a1a6daDerive the losses from physical constants instead of a prescribed decay curveJeremy Magland 2// damped oscillator q'' + 2 sigma q' + omega^2 q = f, with sigma the sum of
3// the metal's internal friction, the thermoelastic loss, the radiation of
4// the mode (damping.ts) and the pads' loss factor on a stand. The stick is
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 5// the Hertz contact of the full simulation, stepped against the modal
6// system with an exact update of every oscillator under a piecewise
7// constant force. Once the stick has left, every mode is a damped sinusoid
8// with a known complex amplitude, and the delayed Rayleigh sum over the
9// plate collapses to one complex gain per mode. During the attack, while the
10// modes are being driven, the delayed sum is evaluated directly from a
11// table of the nodal accelerations. The parts that are sums over the modes
12// (gains, the acceleration table, the audio tail, the movie) accept a model
13// restricted to a subset of the modes, so that a pool of workers can split
14// the modes between them (synthPool.ts).
16import { type LinearSystem } from "../model/assemble.ts";
17import { type ModeSet } from "./modes.ts";
2a1a6daDerive the losses from physical constants instead of a prescribed decay curveJeremy Magland 18import { type Plate, type SimParams, STICK_RELEASE, listenerTerms, standRates, strikePatch, nodeMass } from "../params.ts";
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 19
20/** What the synthesis needs of a plate and its modes; plain arrays so that it can be posted to a worker. */
21export interface ModalModel {
22 nGrid: number;
23 /** Number of plate nodes (rows of wShapes). */
24 nRows: number;
25 count: number;
26 freqs: Float64Array;
27 /** Transverse component of the mass-normalized mode shapes, node-major: wShapes[row * count + m]. */
28 wShapes: Float32Array;
29 /** Row of each grid node in wShapes, or -1 outside the plate. */
30 wRow: Int32Array;
31 /** Grid node of each row. */
32 rowNode: Int32Array;
2a1a6daDerive the losses from physical constants instead of a prescribed decay curveJeremy Magland 33 /** Decay rate of each mode from the metal's internal friction [1/s]. */
34 sigmaMat: Float64Array;
35 /** Decay rate of each mode from thermoelastic damping [1/s]; filled by damping.ts, zero until then. */
36 sigmaTe: Float64Array;
37 /** Decay rate of each mode from its radiation into the air [1/s]; filled by damping.ts, zero until then. */
38 sigmaRad: Float64Array;
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 39 /** Decay rate added by the pads' loss factor to each mode [1/s]; zero on a free plate. */
40 padSigma: Float64Array;
41 /** Node masses on the grid [kg]. */
42 masses: Float64Array;
43}
45export function modalModel(pl: Plate, sys: LinearSystem, modes: ModeSet): ModalModel {
46 const nGrid = pl.grid.n;
47 const wRow = new Int32Array(nGrid).fill(-1);
48 const rows: number[] = [];
49 for (let k = 0; k < nGrid; k++) if (sys.nodeDof[k] >= 0) { wRow[k] = rows.length; rows.push(k); }
50 const nRows = rows.length, M = modes.count;
51 const wShapes = new Float32Array(nRows * M);
52 for (let r = 0; r < nRows; r++) {
53 const d = sys.nodeDof[rows[r]] + (sys.fields - 1);
54 for (let m = 0; m < M; m++) wShapes[r * M + m] = modes.shapes[m * sys.n + d];
55 }
56 const masses = new Float64Array(nGrid);
57 for (let j = 0; j < pl.grid.ny; j++) for (let i = 0; i < pl.grid.nx; i++) masses[j * pl.grid.nx + i] = nodeMass(pl, i, j);
2a1a6daDerive the losses from physical constants instead of a prescribed decay curveJeremy Magland 58 // Hysteretic damping: a loss factor eta on a store of energy U gives a mode
59 // the rate sigma = eta U / omega. The pads hold the energy phi' K_pad phi / 2
60 // and the plate the rest of the mode's strain energy, omega^2 / 2 for a
61 // mass-normalized mode; each has its own loss factor. The thermoelastic and
62 // the radiation rates need the shapes on the grid and are filled in by
63 // damping.ts.
64 const padSigma = new Float64Array(M), sigmaMat = new Float64Array(M);
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 65 const sr = standRates(pl);
2a1a6daDerive the losses from physical constants instead of a prescribed decay curveJeremy Magland 66 for (let m = 0; m < M; m++) {
67 let e = 0;
68 if (sr) {
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 69 for (const k of sr.nodes) {
70 const d = sys.nodeDof[k];
71 if (d < 0) continue;
72 const w = modes.shapes[m * sys.n + d + (sys.fields - 1)];
73 e += sr.kw[k] * w * w;
74 if (sys.fields === 3) {
75 const u = modes.shapes[m * sys.n + d], v = modes.shapes[m * sys.n + d + 1];
76 e += sr.ku[k] * (u * u + v * v);
77 }
78 }
79 }
2a1a6daDerive the losses from physical constants instead of a prescribed decay curveJeremy Magland 80 const omega = 2 * Math.PI * modes.freqs[m];
81 if (!(omega > 0)) continue;
82 padSigma[m] = sr ? (sr.eta * e) / (2 * omega) : 0;
83 sigmaMat[m] = Math.max(0, pl.mat.eta * (omega / 2 - e / (2 * omega)));
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 84 }
2a1a6daDerive the losses from physical constants instead of a prescribed decay curveJeremy Magland 85 return { nGrid, nRows, count: M, freqs: modes.freqs, wShapes, wRow, rowNode: Int32Array.from(rows), sigmaMat, sigmaTe: new Float64Array(M), sigmaRad: new Float64Array(M), padSigma, masses };
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 86}
88/** The model restricted to the modes [from, to). */
89export function sliceModel(model: ModalModel, from: number, to: number): ModalModel {
90 const M = to - from;
91 const wShapes = new Float32Array(model.nRows * M);
92 for (let r = 0; r < model.nRows; r++) wShapes.set(model.wShapes.subarray(r * model.count + from, r * model.count + to), r * M);
2a1a6daDerive the losses from physical constants instead of a prescribed decay curveJeremy Magland 93 return { ...model, count: M, freqs: model.freqs.slice(from, to), wShapes, sigmaMat: model.sigmaMat.slice(from, to), sigmaTe: model.sigmaTe.slice(from, to), sigmaRad: model.sigmaRad.slice(from, to), padSigma: model.padSigma.slice(from, to) };
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 94}
96/** The transverse shape of mode m on the grid (zero outside the plate). */
97export function modeShapeOnGrid(model: ModalModel, m: number): Float32Array {
98 const w = new Float32Array(model.nGrid);
99 for (let r = 0; r < model.nRows; r++) w[model.rowNode[r]] = model.wShapes[r * model.count + m];
100 return w;
101}
2a1a6daDerive the losses from physical constants instead of a prescribed decay curveJeremy Magland 103/** Per-mode decay rates: internal friction, thermoelastic, radiation and the pads, kept below the mode's frequency (underdamped). */
104export function modalRates(model: ModalModel): Float64Array {
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 105 const sigma = new Float64Array(model.count);
106 for (let m = 0; m < model.count; m++) {
107 const w = 2 * Math.PI * model.freqs[m];
2a1a6daDerive the losses from physical constants instead of a prescribed decay curveJeremy Magland 108 sigma[m] = Math.min(model.sigmaMat[m] + model.sigmaTe[m] + model.sigmaRad[m] + model.padSigma[m], 0.99 * w);
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 109 }
110 return sigma;
111}
113/** Damped angular frequency of a mode. */
114function dampedOmega(f: number, sigma: number): number {
115 const w = 2 * Math.PI * f;
116 return Math.sqrt(w * w - sigma * sigma);
117}
119/** The contact: modal gains g_m = sum over the patch of (m_k / M_patch) phi_m(k), so that the modal force is -F g_m and the contact coordinate is sum_m g_m q_m. */
120export function strikeGains(model: ModalModel, pl: Plate, p: SimParams): Float64Array {
121 const patch = strikePatch(pl, p);
122 const g = new Float64Array(model.count);
123 for (const k of patch.nodes) {
124 const r = model.wRow[k];
125 if (r < 0) continue;
126 const f = model.masses[k] / patch.mass;
127 for (let m = 0; m < model.count; m++) g[m] += f * model.wShapes[r * model.count + m];
128 }
129 return g;
130}
132/** The state of every mode at the end of the contact phase and the tables recorded during it. */
133export interface ContactResult {
134 /** Time at which the stick is withdrawn, STICK_RELEASE after the contact last opened [s]. */
135 tEnd: number;
136 /** Time at which the contact last opened [s]; every mode decays freely from here on. */
137 tContact: number;
138 /** Complex modal amplitudes A at tEnd: q(t) = Re[A exp(s (t - tEnd))], s = -sigma + i omega_d. */
139 aRe: Float64Array;
140 aIm: Float64Array;
141 /** Modal accelerations sampled at `tableRate` from t = 0 to tEnd, mode-major: qdd[m * tableLen + i]. */
142 qdd: Float32Array;
143 tableLen: number;
144 tableRate: number;
145 /** Modal displacements at the movie frame times before tEnd, mode-major: qFrames[m * nFrames + f]. */
146 qFrames: Float32Array;
147 nFrames: number;
148 /** Stick force history at the fine step [N] and the step [s], for display. */
149 force: Float32Array;
150 dt: number;
151 peakForce: number;
152 contactTime: number;
153}
155/** The contact result restricted to the modes [from, to). */
156export function sliceContact(c: ContactResult, from: number, to: number): ContactResult {
157 return {
158 ...c,
159 aRe: c.aRe.slice(from, to),
160 aIm: c.aIm.slice(from, to),
161 qdd: c.qdd.slice(from * c.tableLen, to * c.tableLen),
162 qFrames: c.qFrames.slice(from * c.nFrames, to * c.nFrames),
163 };
164}
166/**
167 * Steps the stick against the modal system until the contact has been open
168 * for STICK_RELEASE, recording the tables the synthesis needs. Every mode
169 * is advanced by its exact state transition under a constant force over the
170 * step; the stick by symplectic Euler; the force is the Hertz law of the
171 * penetration at the start of the step.
172 */
173export function simulateContact(model: ModalModel, sigma: Float64Array, g: Float64Array, p: SimParams, fps: number, dt = 2e-6, tableRate = 2 * p.fs): ContactResult {
174 const M = model.count;
175 const omega = new Float64Array(M), omegaD = new Float64Array(M);
176 const p11 = new Float64Array(M), p12 = new Float64Array(M), p21 = new Float64Array(M), p22 = new Float64Array(M);
177 for (let m = 0; m < M; m++) {
178 const w = 2 * Math.PI * model.freqs[m];
179 const s = sigma[m];
180 const wd = dampedOmega(model.freqs[m], s);
181 omega[m] = w; omegaD[m] = wd;
182 const e = Math.exp(-s * dt), c = Math.cos(wd * dt), sn = Math.sin(wd * dt);
183 p11[m] = e * (c + (s * sn) / wd);
184 p12[m] = (e * sn) / wd;
185 p21[m] = (-e * w * w * sn) / wd;
186 p22[m] = e * (c - (s * sn) / wd);
187 }
188 const q = new Float64Array(M), qd = new Float64Array(M);
189 const tMax = 0.2;
190 const maxSteps = Math.ceil(tMax / dt);
191 const tableCap = Math.ceil(tMax * tableRate) + 2;
192 const qddTable = new Float32Array(M * tableCap);
193 const frameCap = Math.ceil(tMax * fps) + 2;
194 const qFrames = new Float32Array(M * frameCap);
195 const force = new Float32Array(maxSteps + 1);
196 let tableLen = 0, nFrames = 0, peakForce = 0, contactTime = 0;
197 let z = 0, zd = -p.strikeVelocity; // stick tip, positive away from the plate; it arrives moving toward it
198 let lastOpen = 0, everContact = false, F = 0;
199 let step = 0;
200 const record = (t: number): void => {
201 while (tableLen < tableCap && tableLen / tableRate <= t + 1e-12) {
202 for (let m = 0; m < M; m++) qddTable[m * tableCap + tableLen] = -2 * sigma[m] * qd[m] - omega[m] * omega[m] * q[m] - F * g[m];
203 tableLen++;
204 }
205 while (nFrames < frameCap && nFrames / fps <= t + 1e-12) {
206 for (let m = 0; m < M; m++) qFrames[m * frameCap + nFrames] = q[m];
207 nFrames++;
208 }
209 };
210 for (step = 0; step < maxSteps; step++) {
211 const t = step * dt;
212 let wc = 0;
213 for (let m = 0; m < M; m++) wc += g[m] * q[m];
214 const delta = wc - z; // penetration: the plate has come past the tip
215 F = delta > 0 ? p.contactStiffness * Math.pow(delta, p.contactExponent) : 0;
216 if (F > 0) { everContact = true; lastOpen = t + dt; contactTime += dt; }
217 if (F > peakForce) peakForce = F;
218 force[step] = F;
219 record(t);
220 if (everContact && t - lastOpen >= STICK_RELEASE) break;
221 for (let m = 0; m < M; m++) {
222 const qp = (-F * g[m]) / (omega[m] * omega[m]);
223 const x = q[m] - qp, xd = qd[m];
224 q[m] = p11[m] * x + p12[m] * xd + qp;
225 qd[m] = p21[m] * x + p22[m] * xd;
226 }
227 zd += (F / p.stickMass) * dt;
228 z += zd * dt;
229 }
230 const tEnd = step * dt;
231 const aRe = new Float64Array(M), aIm = new Float64Array(M);
232 for (let m = 0; m < M; m++) {
233 aRe[m] = q[m];
234 aIm[m] = -(qd[m] + sigma[m] * q[m]) / omegaD[m];
235 }
236 const qdd = new Float32Array(M * tableLen);
237 for (let m = 0; m < M; m++) qdd.set(qddTable.subarray(m * tableCap, m * tableCap + tableLen), m * tableLen);
238 const frames = new Float32Array(M * nFrames);
239 for (let m = 0; m < M; m++) frames.set(qFrames.subarray(m * frameCap, m * frameCap + nFrames), m * nFrames);
240 return { tEnd, tContact: Math.min(lastOpen, tEnd), aRe, aIm, qdd, tableLen, tableRate, qFrames: frames, nFrames, force: force.slice(0, step + 1), dt, peakForce, contactTime };
241}
243// ---------- Listener ----------
245/** Listener terms with the common delay removed, so that only the spread of the delays across the plate remains. */
246export function listenerWeights(pl: Plate, p: SimParams): { weight: Float64Array; delay: Float64Array; maxDelay: number } {
247 const { weight, delay } = listenerTerms(pl, p);
248 let dmin = Infinity;
249 for (let k = 0; k < weight.length; k++) if (weight[k] > 0) dmin = Math.min(dmin, delay[k]);
250 let maxDelay = 0;
251 for (let k = 0; k < weight.length; k++) if (weight[k] > 0) { delay[k] -= dmin; maxDelay = Math.max(maxDelay, delay[k]); } else delay[k] = 0;
252 return { weight, delay, maxDelay };
253}
255/**
256 * The time from which the tail formula holds at the listener: every node's
257 * retarded time is past the last contact. The table path covers the output
258 * samples before it.
259 */
260export function tailStart(c: ContactResult, maxDelay: number): number {
261 return Math.min(c.tEnd, c.tContact + maxDelay + 2 / c.tableRate);
262}
264// ---------- Sums over the modes (accept a sliced model) ----------
266/**
267 * The complex radiation gain of each mode for the listener:
268 * G_m = sum_k weight_k phi_m(k) exp(-s_m delay_k), with s_m = -sigma + i omega_d.
269 */
270export function radiationGains(model: ModalModel, sigma: Float64Array, weight: Float64Array, delay: Float64Array): { re: Float64Array; im: Float64Array } {
271 const M = model.count;
272 const re = new Float64Array(M), im = new Float64Array(M);
273 const wd = new Float64Array(M);
274 for (let m = 0; m < M; m++) wd[m] = dampedOmega(model.freqs[m], sigma[m]);
275 for (let r = 0; r < model.nRows; r++) {
276 const k = model.rowNode[r];
277 const wk = weight[k];
278 if (wk === 0) continue;
279 const tau = delay[k];
280 const base = r * M;
281 for (let m = 0; m < M; m++) {
282 const a = wk * model.wShapes[base + m] * Math.exp(sigma[m] * tau);
283 re[m] += a * Math.cos(wd[m] * tau);
284 im[m] -= a * Math.sin(wd[m] * tau);
285 }
286 }
287 return { re, im };
288}
290/** The nodal accelerations of the model's modes over the first `len` samples of the contact table, row-major over the plate nodes: out[r * len + i] += sum_m phi_m(r) qdd_m(i). */
291export function nodalAccelerationTable(model: ModalModel, c: ContactResult, len: number, out?: Float32Array): Float32Array {
292 const L = c.tableLen, M = model.count;
293 const table = out ?? new Float32Array(model.nRows * len);
294 for (let r = 0; r < model.nRows; r++) {
295 const row = r * len, base = r * M;
296 for (let m = 0; m < M; m++) {
297 const phi = model.wShapes[base + m];
298 if (phi === 0) continue;
299 const src = m * L;
300 for (let i = 0; i < len; i++) table[row + i] += phi * c.qdd[src + i];
301 }
302 }
303 return table;
304}
306/**
307 * The listener pressure from `tStart` on: p(t) = sum_m Re[G_m s_m^2 A_m
308 * exp(s_m (t - tEnd))], generated by the two-term recurrence of a damped
309 * sinusoid, in blocks that stay in cache. Output samples before tStart are
310 * left alone.
311 */
312export function modalTail(model: ModalModel, sigma: Float64Array, c: ContactResult, gains: { re: Float64Array; im: Float64Array }, fs: number, tStart: number, nSamples: number, out?: Float32Array): Float32Array {
313 const y = out ?? new Float32Array(nSamples);
314 const dt = 1 / fs;
315 const block = 4096;
316 const M = model.count;
317 const n0 = Math.ceil(tStart * fs - 1e-9);
318 const a = new Float64Array(M), b = new Float64Array(M), p1 = new Float64Array(M), p2 = new Float64Array(M), limit = new Float64Array(M);
319 for (let m = 0; m < M; m++) {
320 const s = sigma[m], wd = dampedOmega(model.freqs[m], s);
321 // c = G s^2 A, then rotated to the first output sample: c exp(s (n0/fs - tEnd)).
322 const s2Re = s * s - wd * wd, s2Im = -2 * s * wd;
323 const saRe = s2Re * c.aRe[m] - s2Im * c.aIm[m], saIm = s2Re * c.aIm[m] + s2Im * c.aRe[m];
324 let cRe = gains.re[m] * saRe - gains.im[m] * saIm, cIm = gains.re[m] * saIm + gains.im[m] * saRe;
325 const t0 = n0 / fs - c.tEnd;
326 const e0 = Math.exp(-s * t0), c0 = Math.cos(wd * t0), s0 = Math.sin(wd * t0);
327 const r = e0 * (cRe * c0 - cIm * s0), i = e0 * (cRe * s0 + cIm * c0);
328 cRe = r; cIm = i;
329 const e = Math.exp(-s * dt), zRe = e * Math.cos(wd * dt), zIm = e * Math.sin(wd * dt);
330 a[m] = 2 * zRe; b[m] = -(e * e);
331 p2[m] = cRe; // p(n0)
332 p1[m] = cRe * zRe - cIm * zIm; // p(n0 + 1)
333 limit[m] = s > 0 ? Math.ceil((Math.log(1e7) / s) * fs) : Infinity; // 140 dB down: dropped
334 }
335 for (let start = n0; start < nSamples; start += block) {
336 const end = Math.min(nSamples, start + block);
337 for (let m = 0; m < M; m++) {
338 if (start - n0 > limit[m]) continue;
339 const am = a[m], bm = b[m];
340 let x1 = p1[m], x2 = p2[m];
341 let n = start;
342 if (start === n0) {
343 y[n] += x2; n++;
344 if (n < end) { y[n] += x1; n++; }
345 }
346 for (; n < end; n++) {
347 const x = am * x1 + bm * x2;
348 y[n] += x;
349 x2 = x1; x1 = x;
350 }
351 p1[m] = x1; p2[m] = x2;
352 }
353 }
354 return y;
355}
357/**
358 * Movie frames of the transverse displacement: frames[f * nGrid + k] +=
359 * sum_m phi_m(k) q_m(t_f), with q from the contact tables before tEnd and
360 * from the free decay after it.
361 */
362export function movieFrames(model: ModalModel, sigma: Float64Array, c: ContactResult, fps: number, nFrames: number, out?: Float32Array): Float32Array {
363 const frames = out ?? new Float32Array(nFrames * model.nGrid);
364 const M = model.count;
365 const q = new Float64Array(M);
366 const wd = new Float64Array(M);
367 for (let m = 0; m < M; m++) wd[m] = dampedOmega(model.freqs[m], sigma[m]);
368 for (let f = 0; f < nFrames; f++) {
369 const tt = f / fps - c.tEnd;
370 for (let m = 0; m < M; m++) {
371 if (f < c.nFrames) { q[m] = c.qFrames[m * c.nFrames + f]; continue; }
372 const e = Math.exp(-sigma[m] * tt);
373 q[m] = e * (c.aRe[m] * Math.cos(wd[m] * tt) - c.aIm[m] * Math.sin(wd[m] * tt));
374 }
375 const base = f * model.nGrid;
376 for (let r = 0; r < model.nRows; r++) {
377 const rb = r * M;
378 let s = 0;
379 for (let m = 0; m < M; m++) s += model.wShapes[rb + m] * q[m];
380 frames[base + model.rowNode[r]] += s;
381 }
382 }
383 return frames;
384}
386// ---------- Assembly of the parts ----------
388/**
389 * The listener pressure over the attack from the nodal acceleration table:
390 * p(t) = sum_k weight_k wdd_k(t - delay_k), linearly interpolated, for the
391 * output samples before tStart.
392 */
393export function attackPressure(model: ModalModel, table: Float32Array, len: number, tableRate: number, weight: Float64Array, delay: Float64Array, fs: number, tStart: number, out: Float32Array): void {
394 const nOut = Math.min(out.length, Math.ceil(tStart * fs - 1e-9));
395 for (let r = 0; r < model.nRows; r++) {
396 const k = model.rowNode[r];
397 const wk = weight[k];
398 if (wk === 0) continue;
399 const row = r * len;
400 const dk = delay[k];
401 for (let n = 0; n < nOut; n++) {
402 const tr = n / fs - dk;
403 if (tr < 0) continue;
404 const x = tr * tableRate;
405 const i = Math.floor(x);
406 if (i + 1 >= len) continue;
407 const fr = x - i;
408 out[n] += wk * ((1 - fr) * table[row + i] + fr * table[row + i + 1]);
409 }
410 }
411}
413export interface Synthesis {
414 audio: Float32Array;
415 frames: Float32Array;
416 nFrames: number;
417 contact: ContactResult;
418}
420/** Runs the whole synthesis on the calling thread (the tests, and browsers without workers). */
2a1a6daDerive the losses from physical constants instead of a prescribed decay curveJeremy Magland 421export function synthesize(model: ModalModel, pl: Plate, p: SimParams, withMovie = true): Synthesis {
422 const sigma = modalRates(model);
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 423 const g = strikeGains(model, pl, p);
424 const contact = simulateContact(model, sigma, g, p, p.fps);
425 const { weight, delay, maxDelay } = listenerWeights(pl, p);
426 const tStart = tailStart(contact, maxDelay);
427 const len = Math.min(contact.tableLen, Math.ceil(tStart * contact.tableRate) + 2);
428 const nSamples = Math.round(p.duration * p.fs);
429 const audio = new Float32Array(nSamples);
430 const table = nodalAccelerationTable(model, contact, len);
431 attackPressure(model, table, len, contact.tableRate, weight, delay, p.fs, tStart, audio);
432 const gains = radiationGains(model, sigma, weight, delay);
433 modalTail(model, sigma, contact, gains, p.fs, tStart, nSamples, audio);
434 const nFrames = withMovie ? Math.round(p.duration * p.fps) : 0;
435 const frames = withMovie ? movieFrames(model, sigma, contact, p.fps, nFrames) : new Float32Array(0);
436 return { audio, frames, nFrames, contact };
437}