8354482cymmodes: linear modes of a cymsim cymbal driven in the full nonlinear model, with the energy of every mode over timeJeremy Magland 1// Turning a selected linear mode into something the nonlinear simulator can run.
2//
3// Free decay: the shell is released from the mode shape at a chosen amplitude (given as the peak
4// transverse displacement at the antinode, in units of the local thickness) and rings down. The
5// in-plane displacements and rotations of the initial state come from the static condensation
6// (the worker's "complement" request), so the released state is the complete linearised mode and
7// not a pure-w field that would also launch fast membrane and rotation transients.
8//
9// Forcing: a distributed force F_n = Gamma * Mw_n * phi_n * env(t) * sin(omega_d t). Weighting by
10// the transverse lumped mass makes the force M-orthogonal to every other mode, so in linear
11// theory only the driven mode responds; whatever else appears in the nonlinear run got there
12// through the coupling. Gamma is chosen so that the *linear* steady-state amplitude at resonance
13// equals the requested one: Gamma = q_target * 2 sigma omega. With the damping scaled down this
14// target is exceeded (at zero damping the linear response grows without bound), which is stated
15// in the UI rather than hidden.
16import type { ModesData } from "./modes.ts";
17import type { DriveParams, ProjectionSpec } from "./model/simulator.ts";
19export type Mechanism = "decay" | "forced";
21export interface DriveSpec {
22 mode: number; // index into the modes file
23 mechanism: Mechanism;
24 amplitude: number; // peak w at the antinode, in units of the local thickness
25 detuneCents: number; // forced: drive frequency offset from the mode frequency
26 rampTime: number; // forced: seconds of linear onset
27 driveTime: number; // forced: drive stops after this many seconds (free decay afterwards)
28}
30/** Mode shape m on the nodes, dequantised, times ``factor`` (the density rescaling phiFactor). */
31export function modeShape(modes: ModesData, m: number, factor = 1): Float64Array {
32 const N = modes.nNodes, Nm = modes.nModes, out = new Float64Array(N);
33 const s = modes.phiScale[m] * factor;
34 for (let n = 0; n < N; n++) out[n] = modes.phiQ[n * Nm + m] * s;
35 return out;
36}
38/** Node of largest |phi|. */
39export function antinode(shape: Float64Array): { node: number; absPhi: number } {
40 let node = 0, absPhi = 0;
41 for (let n = 0; n < shape.length; n++) if (Math.abs(shape[n]) > absPhi) { absPhi = Math.abs(shape[n]); node = n; }
42 return { node, absPhi };
43}
45/** Modal-coordinate amplitude that puts the antinode at ``woverh`` thicknesses. */
46export function modalAmplitude(modes: ModesData, shape: Float64Array, woverh: number): number {
47 const { node, absPhi } = antinode(shape);
48 return (woverh * modes.h[node]) / absPhi;
49}
51/** Initial state (q, v) of the free decay: w = a phi, in-plane and rotations from the complement
52 * o (node-major u, v, thx, thy per node, computed for the *unit* shape), velocities zero. */
53export function modalInitialState(shape: Float64Array, o: Float64Array, a: number): { q: Float64Array; v: Float64Array } {
54 const N = shape.length, q = new Float64Array(5 * N), v = new Float64Array(5 * N);
55 for (let n = 0; n < N; n++) {
56 q[5 * n] = a * o[4 * n];
57 q[5 * n + 1] = a * o[4 * n + 1];
58 q[5 * n + 2] = a * shape[n];
59 q[5 * n + 3] = a * o[4 * n + 2];
60 q[5 * n + 4] = a * o[4 * n + 3];
61 }
62 return { q, v };
63}
65/** The distributed drive of the forced mechanism. ``Mw`` is the transverse lumped mass of the
66 * simulator's shell (current material), ``sigma`` the decay rate used for the amplitude target,
67 * ``omegaMode`` the (current-material) angular frequency of the mode. */
68export function buildDrive(spec: DriveSpec, modes: ModesData, shape: Float64Array, Mw: Float64Array,
69 sigma: number, omegaMode: number, dt: number): DriveParams {
70 const qTarget = modalAmplitude(modes, shape, spec.amplitude);
71 const gamma = qTarget * 2 * Math.max(sigma, 1e-6) * omegaMode;
72 const N = shape.length, fNode = new Float64Array(N);
73 for (let n = 0; n < N; n++) fNode[n] = gamma * Mw[n] * shape[n];
74 const omega = omegaMode * Math.pow(2, spec.detuneCents / 1200);
75 return {
76 fNode, omega,
77 rampSteps: Math.max(1, Math.round(spec.rampTime / dt)),
78 stopStep: Math.max(1, Math.round(spec.driveTime / dt)),
79 };
80}
82/** The projection matrix (mode-major, mass-weighted): with mass-normalised shapes and the same
83 * shell mass, q_m = sum_n phiM[m N + n] w_n is the modal coordinate. ``Mw`` must be the
84 * transverse lumped mass of the shell actually being simulated. */
85export function buildProjection(modes: ModesData, Mw: Float64Array, factor: number, every: number): ProjectionSpec {
86 const N = modes.nNodes, Nm = modes.nModes;
87 const phiM = new Float32Array(Nm * N);
88 for (let m = 0; m < Nm; m++) {
89 const s = modes.phiScale[m] * factor, row = m * N;
90 for (let n = 0; n < N; n++) phiM[row + n] = modes.phiQ[n * Nm + m] * s * Mw[n];
91 }
92 return { phiM, nModes: Nm, every };
93}