1e29264cymsimweb: WebGPU companion of cymsim with YAML/form editor, live spectrogram and playbackJeremy Magland 1// Exciters and contact patches, as in cymsim/exciter.py: a rigid mass on P power-law springs
2// F = k delta^alpha along a line (stick) or on a small disk (mallet), released at its hit time
3// just above the moving surface. All hits are stacked into arrays padded to the largest patch.
4import type { Hit } from "../config/schema.ts";
5import type { CymbalProfile } from "./geometry.ts";
6import { planeResiduals } from "./linalg.ts";
7import { locate, type Mesh } from "./mesh.ts";
9export const HOLD_HEIGHT = 1e3;
11export interface ContactPatch {
12 hit: Hit; x: Float64Array; y: Float64Array; nodes: Uint32Array; weights: Float64Array; offsets: Float64Array;
13 kTotal: number; alpha: number; mass: number; gap: number; velocity: number; time: number;
14}
16export function patchEnergy(p: ContactPatch): number { return 0.5 * p.mass * p.velocity * p.velocity; }
17export function maxCompression(p: ContactPatch): number { return Math.pow(((p.alpha + 1) * patchEnergy(p)) / p.kTotal, 1 / (p.alpha + 1)); }
19/** Contact time against a rigid wall, 2 c(alpha) d_max / v with c = int_0^1 dx / sqrt(1 - x^(alpha+1)). */
20export function freeContactTime(p: ContactPatch): number {
21 // substitute x = 1 - s^2 to remove the endpoint singularity, then Simpson on s in [0, 1]
22 const a1 = p.alpha + 1, n = 2000;
23 const f = (s: number) => (s === 0 ? 2 / Math.sqrt(a1) : (2 * s) / Math.sqrt(1 - Math.pow(1 - s * s, a1)));
24 let sum = f(0) + f(1);
25 for (let i = 1; i < n; i++) sum += (i % 2 ? 4 : 2) * f(i / n);
26 const c = sum / (3 * n);
27 return (2 * c * maxCompression(p)) / p.velocity;
28}
30export function buildPatch(hit: Hit, mesh: Mesh, profile: CymbalProfile): ContactPatch {
31 const ex = hit.exciter;
32 const xs = hit.r * Math.cos(hit.phi), ys = hit.r * Math.sin(hit.phi);
33 let px: number[], py: number[], kTotal: number, alpha: number;
34 if (ex.kind === "stick") {
35 const ang = hit.phi + ex.patch_angle, P = Math.max(1, Math.floor(ex.patch_points));
36 px = []; py = [];
37 for (let i = 0; i < P; i++) {
38 const xi = (P === 1 ? 0 : -0.5 + i / (P - 1)) * ex.patch_length;
39 px.push(xs + xi * Math.cos(ang)); py.push(ys + xi * Math.sin(ang));
40 }
41 kTotal = ex.k; alpha = ex.alpha;
42 } else {
43 alpha = ex.alpha;
44 kTotal = ((alpha + 1) * 0.5 * ex.mass * hit.velocity ** 2) / ex.compression ** (alpha + 1);
45 px = [xs]; py = [ys];
46 for (let k = 0; k < 6; k++) { const th = (k * Math.PI) / 3; px.push(xs + ex.radius * Math.cos(th)); py.push(ys + ex.radius * Math.sin(th)); }
47 }
48 const P = px.length, nodes = new Uint32Array(3 * P), weights = new Float64Array(3 * P), z0p = new Float64Array(P);
49 for (let p = 0; p < P; p++) {
50 const loc = locate(mesh, px[p], py[p]);
51 nodes.set(loc.nodes, 3 * p); weights.set(loc.weights, 3 * p);
52 z0p[p] = profile.z0(px[p], py[p]);
53 }
54 const offsets = planeResiduals(px.map((x) => x - xs), py.map((y) => y - ys), z0p);
55 return { hit, x: Float64Array.from(px), y: Float64Array.from(py), nodes, weights, offsets, kTotal, alpha, mass: ex.mass,
56 gap: ex.gap, velocity: hit.velocity, time: hit.time };
57}
59export interface HitSchedule {
60 patches: ContactPatch[]; K: number; P: number;
61 nodes: Uint32Array; weights: Float64Array; offsets: Float64Array; // (K, P, 3), (K, P, 3), (K, P)
62 kPoint: Float64Array; alpha: Float64Array; mass: Float64Array; gap: Float64Array; velocity: Float64Array;
63 releaseStep: Uint32Array;
64}
66export function buildSchedule(patches: ContactPatch[], dt: number): HitSchedule {
67 const K = patches.length, P = Math.max(...patches.map((p) => p.x.length));
68 const nodes = new Uint32Array(K * P * 3), weights = new Float64Array(K * P * 3), offsets = new Float64Array(K * P).fill(-1e6);
69 patches.forEach((p, k) => {
70 const n = p.x.length;
71 nodes.set(p.nodes, (k * P) * 3); weights.set(p.weights, (k * P) * 3); offsets.set(p.offsets, k * P);
72 void n;
73 });
74 return {
75 patches, K, P, nodes, weights, offsets,
76 kPoint: Float64Array.from(patches, (p) => p.kTotal / p.x.length),
77 alpha: Float64Array.from(patches, (p) => p.alpha),
78 mass: Float64Array.from(patches, (p) => p.mass),
79 gap: Float64Array.from(patches, (p) => p.gap),
80 velocity: Float64Array.from(patches, (p) => p.velocity),
81 releaseStep: Uint32Array.from(patches, (p) => Math.round(p.time / dt)),
82 };
83}