1// Interactive tutorial: how the cymbal2 solver and renderer work, from the
2// physics of Samejima (2021) down to the step loop the app executes. The
3// prose cites the paper's equation numbers so it can be read side by side
4// with the paper itself.
6import { M, Eq } from './math';
7import GeometryFig from './figures/GeometryFig';
8import PitchGlideFig from './figures/PitchGlideFig';
9import CascadeFig from './figures/CascadeFig';
10import GridFig from './figures/GridFig';
11import WavenumberFig from './figures/WavenumberFig';
12import StabilityFig from './figures/StabilityFig';
13import HertzFig from './figures/HertzFig';
14import WasherFig from './figures/WasherFig';
16const PAPER_DOI = 'https://doi.org/10.1250/ast.42.314';
18const SECTIONS = [
19 { id: 'tut-overview', n: 1, title: 'What a render computes' },
20 { id: 'tut-nonlinear', n: 2, title: 'Why the model must be nonlinear' },
21 { id: 'tut-shell', n: 3, title: 'The shell and the von Kármán equations' },
22 { id: 'tut-boundaries', n: 4, title: 'The washer and the free edge' },
23 { id: 'tut-mallet', n: 5, title: 'The mallet: head, stick, grip, hand' },
24 { id: 'tut-grid', n: 6, title: 'Discretizing space: the grid' },
25 { id: 'tut-wavenumber', n: 7, title: 'One wavenumber at a time' },
26 { id: 'tut-time', n: 8, title: 'Marching in time' },
27 { id: 'tut-sound', n: 9, title: 'From motion to sound' },
28 { id: 'tut-departures', n: 10, title: 'Where the implementation departs from the paper' },
29 { id: 'tut-refs', n: 11, title: 'Credits and references' },
30];
32function PipelineDiagram() {
33 const box = (
34 x: number, y: number, w: number, h: number, name: string, sub: string[],
35 ) => (
36 <g key={name}>
37 <rect x={x} y={y} width={w} height={h} rx={6} fill="#23272f" stroke="#31363f" />
38 <text x={x + w / 2} y={y + 17} textAnchor="middle" fontSize="11.5"
39 fontFamily="ui-monospace, Menlo, monospace" fill="#e8a33d">{name}</text>
40 {sub.map((s, i) => (
41 <text key={i} x={x + w / 2} y={y + 32 + i * 13} textAnchor="middle" fontSize="10.5"
42 fill="#8b93a1">{s}</text>
43 ))}
44 </g>
45 );
46 const arrow = (x1: number, y1: number, x2: number, y2: number) => (
47 <g key={`${x1}-${y1}-${x2}-${y2}`}>
48 <line x1={x1} y1={y1} x2={x2} y2={y2} stroke="#8b93a1" strokeWidth={1.3}
49 markerEnd="url(#tutArrow)" />
50 </g>
51 );
52 return (
53 <svg viewBox="0 0 880 210" className="tut-pipe" role="img"
54 aria-label="Render pipeline: grid, per-wavenumber operators, assembled system and mallet feed the time loop, whose pickup signal goes through the output stage">
55 <defs>
56 <marker id="tutArrow" markerWidth="8" markerHeight="8" refX="7" refY="3"
57 orient="auto">
58 <path d="M0,0 L7,3 L0,6 z" fill="#8b93a1" />
59 </marker>
60 </defs>
61 {box(6, 24, 96, 42, 'parameters', ['shell, mallet,', 'stroke'])}
62 {box(134, 18, 142, 54, 'cymbal2_grid.m', ['polar grid,', 'azimuthal transform'])}
63 {box(308, 18, 150, 54, 'cymbal2_perm.m', ['radial operators', 'per wavenumber'])}
64 {box(490, 12, 168, 66, 'cymbal2_system.m', ['implicit update +', 'boundary rows,', 'inverted'])}
65 {box(134, 130, 168, 54, 'cymbal2_mallet.m', ['stick beam FEM,', 'grip, head mass'])}
66 {box(490, 120, 168, 70, 'cymbal2_sim.m', ['time loop: hand,', 'contact, von Kármán', 'shell, step by step'])}
67 {box(700, 126, 172, 58, 'cymbal2_output.m', ['DC block, decimate,', 'normalize → audio'])}
68 {arrow(102, 45, 130, 45)}
69 {arrow(276, 45, 304, 45)}
70 {arrow(458, 45, 486, 45)}
71 {arrow(574, 78, 574, 116)}
72 {arrow(302, 157, 486, 157)}
73 {arrow(658, 155, 696, 155)}
74 </svg>
75 );
76}
78function MalletDiagram() {
79 return (
80 <svg viewBox="0 0 640 210" className="tut-pipe" role="img"
81 aria-label="Mallet configuration: hand at a rotational spring drives an elastic beam whose tip carries the head, which meets the shell across a gap">
82 {/* hand and rotational spring */}
83 <circle cx={70} cy={70} r={16} fill="none" stroke="#e8a33d" strokeWidth={2}
84 strokeDasharray="20 6" />
85 <text x={70} y={34} textAnchor="middle" fontSize="11" fill="#e8a33d">rotational spring K_R</text>
86 <text x={70} y={122} textAnchor="middle" fontSize="11" fill="#8b93a1">hand: u_H(t) prescribed</text>
87 <line x1={70} y1={92} x2={70} y2={100} stroke="#8b93a1" markerEnd="none" />
88 {/* beam */}
89 <path d="M 86 70 C 220 62, 380 62, 500 78" fill="none" stroke="#d6dae2" strokeWidth={4} />
90 <text x={300} y={48} textAnchor="middle" fontSize="11" fill="#d6dae2">
91 stick: Euler–Bernoulli beam ζ(x, t)
92 </text>
93 {/* head */}
94 <circle cx={512} cy={86} r={13} fill="#e8a33d" />
95 <text x={512} y={62} textAnchor="middle" fontSize="11" fill="#e8a33d">head M_M, u_M</text>
96 {/* gap and force */}
97 <line x1={512} y1={101} x2={512} y2={124} stroke="#4dd0e1" strokeDasharray="3 3" />
98 <text x={528} y={116} fontSize="11" fill="#4dd0e1">ξ (compression)</text>
99 <text x={430} y={116} fontSize="11" fill="#4dd0e1">f = K_M [ξ]₊^α</text>
100 {/* shell */}
101 <path d="M 260 168 Q 420 118 630 158" fill="none" stroke="#c9a227" strokeWidth={4} />
102 <text x={330} y={190} textAnchor="middle" fontSize="11" fill="#c9a227">shell u(r, θ, t)</text>
103 </svg>
104 );
105}
107export default function Tutorial({ onClose }: { onClose: () => void }) {
108 return (
109 <div className="tutorial">
110 <nav className="tut-toc">
111 <div className="tut-toc-head">How it works</div>
112 {SECTIONS.map((s) => (
113 <a key={s.id} href={`#${s.id}`}>
114 {s.n}. {s.title}
115 </a>
116 ))}
117 <button onClick={onClose}>← Back to the simulator</button>
118 </nav>
120 <div className="tut-body">
121 <header className="tut-header">
122 <h1>How the cymbal is computed</h1>
123 <p className="tut-lede">
124 An interactive walk through the theory and the implementation behind this
125 simulator. The physical model is due to <strong>Toshiya Samejima</strong> (Faculty
126 of Design, Kyushu University):{' '}
127 <a href={PAPER_DOI} target="_blank" rel="noreferrer">
128 “Nonlinear physical modeling sound synthesis of cymbals involving dynamics
129 of washers and sticks/mallets,” <i>Acoustical Science and Technology</i>{' '}
130 <b>42</b>(6), 314–325 (2021)
131 </a>
132 , which itself builds on Bilbao's nonlinear shell synthesis (2010). Equation
133 numbers below, written (1), (20), and so on, are the paper's, so this page can be
134 read beside it. Where the implementation departs from the paper in order to run
135 stably or quickly, Section 10 says exactly how.
136 </p>
137 <button className="primary" onClick={onClose}>← Back to the simulator</button>
138 </header>
140 {/* ---------------- 1 ---------------- */}
141 <section id="tut-overview">
142 <h2>1. What a render computes</h2>
143 <p>
144 Nothing in this app is a recording, an oscillator bank, or a filter tuned to sound
145 like a cymbal. Pressing <i>Render</i> integrates the equations of motion of three
146 coupled mechanical systems, sample by sample, several hundred thousand times: a
147 thin bronze shell (the cymbal itself), the felt washer it is mounted on, and the
148 mallet striking it, the mallet being a flexible stick held by a moving hand with
149 an elastic head on its tip. The sound you hear is the vertical acceleration of one
150 point on the shell's surface, read off the simulation like a contact pickup.
151 </p>
152 <p>
153 The whole computation is written in MATLAB syntax, in the <code>.m</code> files
154 you can open (and edit) with the app's <i>Code</i> button, and runs in your
155 browser. Each stage lives in its own file:
156 </p>
157 <PipelineDiagram />
158 <p>
159 The first three boxes are precomputation: they build the grid and assemble, then
160 invert, the small matrices the time loop applies at every step. The loop itself,
161 in <code>cymbal2_sim.m</code>, is where the minutes go (or the seconds, when the
162 browser exposes WebGPU and the loop runs on the GPU). The sections below follow
163 this pipeline: first the physics being solved (Sections 2 to 5), then how it is
164 discretized and integrated (Sections 6 to 8), then how the motion becomes audio
165 (Section 9).
166 </p>
167 </section>
169 {/* ---------------- 2 ---------------- */}
170 <section id="tut-nonlinear">
171 <h2>2. Why the model must be nonlinear</h2>
172 <p>
173 A linear model of a struck plate is a fixed set of modes: the strike decides how
174 much energy each mode receives, and from then on each rings at its own constant
175 frequency and decays. Whatever spectrum exists at the moment of impact is the
176 spectrum forever after, only quieter. A real cymbal does something a linear model
177 cannot: struck hard, its sound <i>changes shape</i> after the strike. Pitches
178 glide, and the bright, noise-like wash of a crash builds up over several hundred
179 milliseconds, well after the stick has left. Chaigne, Touzé and Thomas
180 called cymbals “nonlinear percussion instruments” for this reason, and
181 Samejima's paper opens with the same observation.
182 </p>
183 <p>
184 The mechanism is geometric. A flat plate can bend without stretching, but a plate
185 deflected into a curve, or a shell deflected at all, must stretch its own
186 mid-surface, and the membrane tension this creates acts back on the bending
187 motion. The effect grows as the square of the deflection measured in plate
188 thicknesses: negligible for a whisper-quiet tap, dominant for a crash, where the
189 shell moves by tens of thicknesses. Two small experiments show what this coupling
190 does, and both are audible: each play button integrates the system's equations
191 right here in the page and renders the result. The first is the simplest system
192 that has it, a single oscillator whose spring stiffens with displacement:
193 </p>
194 <PitchGlideFig />
195 <p>
196 The second experiment shows the other, more consequential effect: nonlinearity
197 couples modes to each other, so energy placed low in the spectrum migrates upward.
198 Here is a chain of masses with a quadratic geometric term in its springs (the
199 Fermi–Pasta–Ulam system), started with all of its energy in its lowest
200 mode:
201 </p>
202 <CascadeFig />
203 <p>
204 Both effects, the amplitude-dependent detuning and the upward cascade, emerge in
205 the shell model of the next section from a single nonlinear term. Neither is put
206 in by hand. This is the sense in which the crash is computed rather than modeled:
207 hit the simulated shell gently and it rings like a bell; hit it hard and the
208 cascade fills the spectrum, because the equations say so.
209 </p>
210 </section>
212 {/* ---------------- 3 ---------------- */}
213 <section id="tut-shell">
214 <h2>3. The shell and the von Kármán equations</h2>
215 <p>
216 A real cymbal has a raised cup and a gently tapered bow. Following Bilbao (2010),
217 the paper replaces both by one <i>shallow spherical cap</i> of radius{' '}
218 <M>a</M>, uniform thickness <M>h</M>, and a small rise at the center. The rise
219 fixes the radius <M>R_s</M> of the sphere the cap sits on, and that single
220 curvature is all the geometry the model carries:
221 </p>
222 <GeometryFig />
223 <p>
224 The price of this simplification is real and audible: the linear spectrum is that
225 of an idealized cap, not of a hammered, lathed instrument, and no irregularity
226 appears anywhere in the model. What it buys is a formulation in which the strong
227 nonlinearity is actually computable. The transverse deflection{' '}
228 <M>u(r, \theta, t)</M> obeys the von Kármán shallow-shell equations
229 (the paper's Eqs. (1) and (2)):
230 </p>
231 <Eq tag="(1)">
232 {String.raw`\rho h\,\frac{\partial^2 u}{\partial t^2}
233 = -D\nabla^4 u \;+\; L(\phi, u)\;-\;\frac{1}{R_s}\nabla^2\phi
234 \;-\;2\sigma_0\rho h\,\frac{\partial u}{\partial t}
235 \;+\;2\sigma_1 a^2 \rho h\,\nabla^2\frac{\partial u}{\partial t}
236 \;+\;\epsilon f`}
237 </Eq>
238 <Eq tag="(2)">
239 {String.raw`\nabla^4\phi = -\frac{Eh}{2}\,L(u, u)\;+\;\frac{Eh}{R_s}\nabla^2 u`}
240 </Eq>
241 <p>
242 Reading Eq. (1) term by term: <M>\rho h</M> is the mass per unit area;{' '}
243 <M>-D\nabla^4 u</M> is ordinary linear plate bending, with{' '}
244 <M>D = Eh^3/12(1-\nu^2)</M> the bending stiffness and <M>\nabla^4</M> the
245 biharmonic operator (the Laplacian applied twice); the two damping terms carry the
246 fitted loss constants <M>\sigma_0</M> (frequency-independent) and{' '}
247 <M>\sigma_1</M> (growing with the squared wavenumber, so high frequencies die
248 faster, as they do on real cymbals); and <M>\epsilon f</M> is the mallet's contact
249 force, applied at the strike point.
250 </p>
251 <p>
252 The remaining two terms are the shell physics. The function <M>\phi</M> is the{' '}
253 <i>Airy stress function</i>: a single scalar field whose derivatives encode the
254 membrane (in-plane) stress everywhere in the shell. Eq. (2) says where that stress
255 comes from: deflecting the shell stretches it, through the curvature (the{' '}
256 <M>\nabla^2 u / R_s</M> term) and through the quadratic operator{' '}
257 <M>L(u,u)</M>. Eq. (1) says what the stress does: it feeds back on the bending
258 motion, again through the curvature and through <M>L(\phi, u)</M>. The bracket{' '}
259 <M>L</M> is a fixed bilinear form; in polar coordinates,
260 </p>
261 <Eq>
262 {String.raw`L(\phi, u) = \frac{\partial^2\phi}{\partial r^2}
263 \left(\frac{1}{r}\frac{\partial u}{\partial r} + \frac{1}{r^2}\frac{\partial^2 u}{\partial\theta^2}\right)
264 + \frac{\partial^2 u}{\partial r^2}
265 \left(\frac{1}{r}\frac{\partial\phi}{\partial r} + \frac{1}{r^2}\frac{\partial^2\phi}{\partial\theta^2}\right)
266 - 2\left(\frac{1}{r}\frac{\partial^2\phi}{\partial r\,\partial\theta} - \frac{1}{r^2}\frac{\partial\phi}{\partial\theta}\right)
267 \left(\frac{1}{r}\frac{\partial^2 u}{\partial r\,\partial\theta} - \frac{1}{r^2}\frac{\partial u}{\partial\theta}\right)`}
268 </Eq>
269 <p>
270 This closed loop, deflection creates stress (2), stress deflects (1), is the
271 engine of everything in Section 2. Because <M>L(u,u)</M> is quadratic in the
272 deflection, the feedback is negligible at small amplitude and dominant at large,
273 and because <M>L</M> mixes radial and azimuthal derivatives, it couples modes of
274 different shape, which is what drives the cascade. The paper works with a
275 dimensionless form of these equations (its Eqs. (3) to (5)); the app integrates
276 the dimensional form above, which is the same content in physical units.
277 </p>
278 </section>
280 {/* ---------------- 4 ---------------- */}
281 <section id="tut-boundaries">
282 <h2>4. The washer and the free edge</h2>
283 <p>
284 Fourth-order equations need two boundary conditions at each edge. At the outer rim{' '}
285 <M>r = a</M> the shell is free, so the bending moment and the (Kelvin–Kirchhoff)
286 shear force both vanish (the paper's Eqs. (6) and (7)):
287 </p>
288 <Eq tag="(6), (7)">
289 {String.raw`\left.\frac{\partial^2 u}{\partial r^2} + \nu\Bigl(\frac{1}{r}\frac{\partial u}{\partial r} + \frac{1}{r^2}\frac{\partial^2 u}{\partial\theta^2}\Bigr)\right|_{r=a} = 0,
290 \qquad
291 \left.\frac{\partial^3 u}{\partial r^3} + \frac{1}{r}\frac{\partial^2 u}{\partial r^2} - \frac{1}{r^2}\frac{\partial u}{\partial r}
292 + \frac{\nu - 3}{r^3}\frac{\partial^2 u}{\partial\theta^2}
293 + \frac{2-\nu}{r^2}\frac{\partial^3 u}{\partial r\,\partial\theta^2}\right|_{r=a} = 0`}
294 </Eq>
295 <p>
296 The inner boundary is the paper's first main contribution. Earlier work (Bilbao
297 2010) clamped the shell rigidly at its center, but a mounted cymbal is nothing
298 like clamped: it sits on a felt washer, and visibly bounces and rocks on it when
299 struck. Samejima models the washer as a mass–damper–spring per unit
300 length of the mounting circle <M>r = c</M>, and imposes it as a boundary
301 condition: the bending moment vanishes there (the shell may tilt freely), and the
302 shear force balances the washer's reaction (Eqs. (14) and (15)):
303 </p>
304 <Eq tag="(15)">
305 {String.raw`-D\,[\text{shear expression}]_{\,r=c}
306 \;=\; M_\mathrm{W}\frac{\partial^2 u}{\partial t^2}
307 \;+\; R_\mathrm{W}\frac{\partial u}{\partial t}
308 \;+\; K_\mathrm{W}\,u`}
309 </Eq>
310 <p>
311 The paper measures the stiffness of an actual washer,{' '}
312 <M>{String.raw`K_\mathrm{W} = 1.65\times 10^5\ \mathrm{N/m}`}</M> per meter, and sets the mass
313 and damping to zero. The consequences are on two scales. The whole shell acquires
314 rigid-body modes on its mount:
315 </p>
316 <WasherFig />
317 <p>
318 And the audible modes see a compliant inner edge instead of a rigid one, which
319 shifts the low-order resonances; the paper's Fig. 8 shows exactly this between its
320 clamped and impedance-supported calculations, and its Fig. 9 shows the mode shapes
321 that appear only when the support may move, such as the (0,0) mode in which the
322 cymbal moves as one piece. The app defaults to the washer support;{' '}
323 <code>support_clamped</code> restores the rigid clamp of the earlier literature
324 (Eqs. (12), (13)) for comparison.
325 </p>
326 </section>
328 {/* ---------------- 5 ---------------- */}
329 <section id="tut-mallet">
330 <h2>5. The mallet: head, stick, grip, hand</h2>
331 <p>
332 The paper's second main contribution is taking the striking mechanism seriously.
333 Most synthesis work injects a prescribed force pulse (a raised cosine, say) at the
334 strike point. Here the force is not prescribed at all: it is whatever comes out of
335 a mallet, itself a dynamical system, colliding with the moving shell.
336 </p>
337 <MalletDiagram />
338 <p>
339 The <b>head</b> (mass <M>{String.raw`M_\mathrm{M}`}</M>, 28 g for the paper's measured mallet)
340 is an elastic body whose compression <M>\xi</M> against the shell produces a
341 Hertz-type force (Eq. (32)):
342 </p>
343 <Eq tag="(32)">
344 {String.raw`f_{\mathrm{M}1} = K_\mathrm{M}\,\bigl[\,\xi\,\bigr]_+^{\,\alpha},
345 \qquad \xi = u - (u_\mathrm{M} - d_\mathrm{M})`}
346 </Eq>
347 <p>
348 where <M>[\cdot]_+</M> keeps only positive values: the head pushes when
349 compressed and lets go when it separates, and the same force acts back on the
350 shell (Eq. (33)). This one-sided power law is the standard model for
351 hammer–string and mallet–membrane contact (Chaigne & Askenfelt
352 1994; Rhaouti et al. 1999), which the paper adopts for the cymbal:
353 </p>
354 <HertzFig />
355 <p>
356 The <b>stick</b> is an Euler–Bernoulli beam (Eq. (35)), discretized by the
357 finite element method into cubic beam elements (Eq. (38)), with the head attached
358 at its tip. The <b>grip</b> is the boundary condition at the other end, and it is
359 the paper's third experiment: a tight grip is a clamp (displacement and rotation
360 both follow the hand), a loose grip lets the stick pivot against a rotational
361 spring <M>{String.raw`K_\mathrm{R}`}</M>, measured at 140 N·m/rad. Because the beam is
362 compliant, the tip whips: it arrives faster than the hand that drives it (by a
363 factor near 1.3 at the default stroke), flexural waves reflect along the stick
364 during contact, and the contact force develops the multiple peaks of the paper's
365 Fig. 12. All of this is audible in the attack, and the paper connects it to the
366 measured effect of grip tightness on cymbal tone (Takashima et al. 2018).
367 </p>
368 <p>
369 The <b>hand</b> closes the system: the grip end's displacement{' '}
370 <M>{String.raw`u_\mathrm{H}(t)`}</M> is prescribed. The paper drives it at constant acceleration
371 to a terminal speed; this app's hand is a small state machine that accelerates,
372 coasts, holds briefly on contact, and retracts to rest, so that repeated strokes
373 and rolls are possible (Section 10 has the details).
374 </p>
375 </section>
377 {/* ---------------- 6 ---------------- */}
378 <section id="tut-grid">
379 <h2>6. Discretizing space: the grid</h2>
380 <p>
381 To compute with the equations of Section 3, the shell's annulus (washer circle to
382 rim) is sampled on a polar grid: <M>N_r</M> intervals in radius,{' '}
383 <M>N_\theta</M> points around the circumference (the paper's Fig. 4, and{' '}
384 <code>cymbal2_grid.m</code> here). Spatial derivatives become central differences:
385 for example
386 </p>
387 <Eq>
388 {String.raw`\left.\frac{\partial u}{\partial r}\right|_{l}\approx\frac{u_{l+1}-u_{l-1}}{2h_r},
389 \qquad
390 \left.\frac{\partial^2 u}{\partial r^2}\right|_{l}\approx\frac{u_{l+1}-2u_l+u_{l-1}}{h_r^2}`}
391 </Eq>
392 <p>
393 The biharmonic <M>\nabla^4</M> reaches two neighbors deep, so evaluating it on the
394 first real line inside a boundary needs values on lines that lie <i>outside</i>{' '}
395 the shell. The grid therefore carries two <i>virtual</i> lines beyond each
396 boundary, and the four boundary conditions of Section 4 (two per edge) are exactly
397 the equations that determine the values there. Nothing is wasted and nothing is
398 free: the boundary conditions are enforced by construction, at the cost of
399 carrying the virtual unknowns in every solve.
400 </p>
401 <GridFig />
402 <p>
403 The grid also sets what the model can resolve. The paper estimates that with its
404 production grid (24 by 60) the computation is reliable up to roughly the (6,5)
405 mode near 2 kHz; the crash extends far above that, but its top octaves ride on
406 grid modes that are dispersion-warped (Section 7 shows why), which is the main
407 respect in which the simulated crash is duller than a real one. Raising both
408 counts sharpens it, at a cost that grows with their product.
409 </p>
410 </section>
412 {/* ---------------- 7 ---------------- */}
413 <section id="tut-wavenumber">
414 <h2>7. One wavenumber at a time</h2>
415 <p>
416 The implicit time step of the next section needs two elliptic solves per sample:
417 one to find the stress function <M>\phi</M> from Eq. (2), one to advance{' '}
418 <M>u</M>. Done directly, each is a linear system over the whole grid, about 1,500
419 unknowns at the default resolution, several hundred thousand times per render.
420 The paper assembles exactly these full-grid operators (its Eqs. (18), (19) and
421 Appendix, as Kronecker products). The implementation here exploits a symmetry
422 instead, and this is its one structural departure from the paper (Section 10):
423 every linear operator and boundary row in the model is axisymmetric, unchanged by
424 rotating the cymbal. Such operators cannot mix azimuthal wavenumbers, so
425 decomposing everything into <M>\cos m\theta</M> and <M>\sin m\theta</M> patterns
426 splits each big solve into <M>N_\theta</M> independent radial solves of size{' '}
427 <M>N_r + 2</M>, each of which is precomputed and inverted once before the loop
428 starts. The result of each step is identical (the decomposition is exact, and the
429 self-test verifies it); the cost falls by about two orders of magnitude.
430 </p>
431 <WavenumberFig />
432 <p>
433 One subtlety is worth naming. On the discrete ring, differentiating twice in{' '}
434 <M>\theta</M> multiplies the wavenumber-<M>m</M> pattern by{' '}
435 <M>{String.raw`\lambda_2(m) = (2\cos m h_\theta - 2)/h_\theta^2`}</M> rather than
436 by <M>-m^2</M>. At low order the two agree; toward the grid cutoff{' '}
437 <M>m = N_\theta/2</M> the discrete value falls short, so the highest
438 nodal-diameter modes sit flat of where physics would put them. The same happens
439 radially. This numerical dispersion is harmless for the low modes that carry the
440 pitch and wobble, and increasingly wrong for the top of the crash, which is why
441 the tutorial keeps saying the top octaves are approximate.
442 </p>
443 <p>
444 The nonlinear terms do mix wavenumbers, deliberately: they are the cascade. They
445 are therefore evaluated in physical space, on the grid, and transformed to the
446 wavenumber basis only for the solves. Every time step crosses between the two
447 representations twice.
448 </p>
449 </section>
451 {/* ---------------- 8 ---------------- */}
452 <section id="tut-time">
453 <h2>8. Marching in time</h2>
454 <p>
455 Time advances in steps of <M>{String.raw`h_t = 1/(f_s \cdot \text{oversample})`}</M>,
456 about 10 μs at the default 96 kHz simulation rate. The scheme, the paper's Eq.
457 (20), splits the physics by stiffness. The <i>linear</i> terms, bending, curvature
458 coupling, damping, are integrated <i>implicitly</i>: the update solves
459 </p>
460 <Eq tag="(20)">
461 {String.raw`\boldsymbol{A}\,\boldsymbol{u}_{k+1} = \boldsymbol{B}\,\boldsymbol{u}_k + \boldsymbol{C}\,\boldsymbol{u}_{k-1}
462 \;+\; h_t^2\,[\,\text{nonlinear terms} + \text{contact force}\,]_k`}
463 </Eq>
464 <p>
465 where <M>\boldsymbol A</M>, <M>\boldsymbol B</M>, <M>\boldsymbol C</M> contain the
466 biharmonic and Laplacian operators weighted by a free parameter <M>g</M> (the
467 paper uses <M>g = 0.2498</M>), and the boundary rows of Sections 4 replace the
468 equation-of-motion rows on the boundary and virtual lines. The washer's own
469 dynamics, Eq. (15), is folded into two of those rows, so the support costs
470 nothing extra. Why implicit? Because the alternative fails:
471 </p>
472 <StabilityFig />
473 <p>
474 The <i>nonlinear</i> terms and the contact force are evaluated explicitly, at the
475 current step, because <M>L(\phi, u)</M> changes every sample and cannot be
476 pre-inverted. This is the standard compromise (Bilbao 2009), and it has a
477 consequence you will meet if you push the hand speed slider: an explicit quadratic
478 term is only conditionally stable, and the condition tightens with amplitude.
479 Strokes whose tip arrives above roughly 10 m/s drive deflections past forty plate
480 thicknesses and the integration diverges, at any rate we tried. This is also
481 about where the shallow-shell expansion itself stops being valid, so the model
482 refuses more or less where its physics runs out; the renderer detects the
483 divergence and reports it rather than emitting garbage. Less intuitively, the
484 azimuthal resolution matters more than the time step: the cascade climbs the
485 nodal-diameter ladder, and if <code>n_theta</code> cuts the ladder short, energy
486 piles up at the cutoff and blows up where a full grid rides it out.
487 </p>
488 <p>With the operators prebuilt, one pass through the loop in{' '}
489 <code>cymbal2_sim.m</code> is:
490 </p>
491 <ol className="tut-steps">
492 <li>
493 <b>Contact.</b> Compression <M>\xi</M> from the shell and beam-tip positions;
494 force <M>{String.raw`f = K_\mathrm{M}[\xi]_+^\alpha`}</M> if touching, zero if
495 not.
496 </li>
497 <li>
498 <b>Mallet.</b> Advance the stick beam one step (its own implicit two-step
499 scheme), driven by the hand at one end and the contact force at the tip.
500 </li>
501 <li>
502 <b>Hand.</b> Advance the stroke state machine (accelerate, coast, hold,
503 retract).
504 </li>
505 <li>
506 <b>Stress function.</b> Evaluate <M>L(u,u)</M> and <M>\nabla^2 u</M> on the
507 grid, transform to wavenumbers, solve Eq. (2) for <M>\phi</M> (one small radial
508 solve per wavenumber), transform back.
509 </li>
510 <li>
511 <b>Nonlinear forcing.</b> Evaluate <M>L(\phi, u)</M> and the membrane term on
512 the grid; inject the contact force at the strike node.
513 </li>
514 <li>
515 <b>Shell.</b> Transform the total forcing, apply the prebuilt{' '}
516 <M>{String.raw`\boldsymbol A^{-1}\boldsymbol B`}</M>,{' '}
517 <M>{String.raw`\boldsymbol A^{-1}\boldsymbol C`}</M>,{' '}
518 <M>{String.raw`\boldsymbol A^{-1}`}</M> per wavenumber, transform back:{' '}
519 <M>{String.raw`\boldsymbol u_{k+1}`}</M>.
520 </li>
521 <li>
522 <b>Pickup.</b> Record the second time difference of <M>u</M> at the pickup
523 node: the acceleration sample.
524 </li>
525 </ol>
526 <p>
527 A full render at the defaults is about 580,000 such steps. In the browser the
528 same loop runs either through the numbl interpreter or, when WebGPU is available,
529 as a sequence of GPU compute passes per step, which is roughly twenty times
530 faster; the diagnostics panel reports which path a render took.
531 </p>
532 </section>
534 {/* ---------------- 9 ---------------- */}
535 <section id="tut-sound">
536 <h2>9. From motion to sound</h2>
537 <p>
538 The rendered signal is the shell's normal acceleration at the pickup point,
539 just as in the paper's numerical experiments (its output point is{' '}
540 <M>\bar r = 0.6</M> on the meridian opposite the strike, and that is the app's
541 default too). <code>cymbal2_output.m</code> conditions it: a gentle DC blocker
542 removes the quasi-static offset the washer relaxation leaves; a windowed-sinc
543 filter anti-aliases and decimates from the simulation rate down to the output
544 rate; a short fade closes the tail; and the result is normalized to the requested
545 peak level.
546 </p>
547 <p>
548 It is worth being precise about what this signal is and is not. It is the motion
549 of one point of the metal: a contact pickup, not a microphone. Radiating the
550 sound into air, with the interference and directivity of a real plate source,
551 would require the full surface velocity history and an acoustic model on top, and
552 is not part of the paper's method. The audible consequence is that low-order
553 modes which radiate poorly from a real cymbal are present here at full strength,
554 so renders carry somewhat more low end than a room recording of the same
555 instrument would. With that caveat, the paper's trial calculations, and this
556 implementation against the demonstration recording published by the paper's
557 authors, reproduce the features that define the instrument: the clear low-order
558 peaks, the slow buildup of high-frequency energy after a hard stroke, the
559 difference between a washer-mounted and a clamped cymbal, and the effect of the
560 grip on the attack. The spectrogram panel of the simulator is the place to watch
561 the first of these; render the <i>Paper demo</i> preset and look for the energy
562 arriving above 5 kHz over the first half second.
563 </p>
564 </section>
566 {/* ---------------- 10 ---------------- */}
567 <section id="tut-departures">
568 <h2>10. Where the implementation departs from the paper</h2>
569 <p>
570 The physics above is Samejima's throughout. Turning a published scheme into code
571 that must run stably for arbitrary knob settings forced a small number of
572 departures, and honesty requires naming them; the repository README states each
573 with its full justification.
574 </p>
575 <ul className="tut-list">
576 <li>
577 <b>The head is lumped onto the beam tip.</b> The paper integrates the head as
578 its own one-degree body (Eq. (34)) coupled to the beam through a rigid
579 connection (Eqs. (36), (37)). Coupling that explicit update to the implicit
580 beam turns out to be unstable as a discrete pair. Since the connection is
581 rigid, substituting it into Newton's law puts the head mass on the tip node
582 exactly, which is the same physics and inherits the beam scheme's stability.
583 </li>
584 <li>
585 <b>The implicit weighting of the washer and beam schemes.</b> The paper prints
586 the split with weight <M>(1-\phi)/2</M> on the old and new steps; that choice
587 is unstable whenever the mass term does not dominate (the massless washer
588 included). Moving the implicit weight <M>\phi/2</M> onto the average of the new
589 and old steps makes <M>\phi = 0.7</M> unconditionally stable, and that is what
590 runs here.
591 </li>
592 <li>
593 <b>The hand can retract.</b> The paper's hand accelerates to a terminal speed
594 and never reverses, which is fine for one stroke and cannot play a second. The
595 app's hand accelerates at an arm-like rate, coasts, and snaps back after
596 contact, so figures and rolls are possible.
597 </li>
598 <li>
599 <b>The solves are block-diagonalized over wavenumber</b> (Section 7). This
600 changes the cost, not the result.
601 </li>
602 <li>
603 <b>Guard rails.</b> The loop detects divergence and stops with a message rather
604 than emitting garbage, and the system builder truncates azimuthal columns whose
605 inner boundary layer the radial grid cannot resolve (none, at the default
606 grid).
607 </li>
608 </ul>
609 </section>
611 {/* ---------------- 11 ---------------- */}
612 <section id="tut-refs">
613 <h2>11. Credits and references</h2>
614 <p>
615 The physical model implemented by this app, the washer boundary condition, the
616 stick/mallet system, and the numerical scheme, are the work of{' '}
617 <b>Toshiya Samejima</b> (Faculty of Design, Kyushu University), published as the
618 paper below. The defaults are the paper's measured instrument: an 18-inch medium
619 crash on its measured washer, struck by its measured mallet. This tutorial and
620 the implementation are the app's own; errors in either should be laid at the
621 app's door, not the paper's.
622 </p>
623 <ol className="tut-refs">
624 <li>
625 T. Samejima, “Nonlinear physical modeling sound synthesis of cymbals
626 involving dynamics of washers and sticks/mallets,” <i>Acoustical Science
627 and Technology</i> <b>42</b>(6), 314–325 (2021).{' '}
628 <a href={PAPER_DOI} target="_blank" rel="noreferrer">doi:10.1250/ast.42.314</a>
629 </li>
630 <li>
631 S. Bilbao, “Percussion synthesis based on models of nonlinear shell
632 vibration,” <i>IEEE Transactions on Audio, Speech, and Language
633 Processing</i> <b>18</b>, 872–880 (2010).
634 </li>
635 <li>
636 S. Bilbao, <i>Numerical Sound Synthesis: Finite Difference Schemes and
637 Simulation in Musical Acoustics</i> (Wiley, 2009).
638 </li>
639 <li>
640 A. Chaigne, C. Touzé and O. Thomas, “Nonlinear vibrations and
641 chaos in gongs and cymbals,” <i>Acoustical Science and Technology</i>{' '}
642 <b>26</b>, 403–409 (2005).
643 </li>
644 <li>
645 A. Chaigne and A. Askenfelt, “Numerical simulations of piano strings. I.
646 A physical model for a struck string using finite difference methods,”{' '}
647 <i>Journal of the Acoustical Society of America</i> <b>95</b>, 1112–1118
648 (1994).
649 </li>
650 <li>
651 L. Rhaouti, A. Chaigne and P. Joly, “Time-domain modeling and numerical
652 simulation of a kettledrum,” <i>Journal of the Acoustical Society of
653 America</i> <b>105</b>, 3545–3562 (1999).
654 </li>
655 <li>
656 S. Takashima, N. Wakatsuki and K. Mizutani, “Measurement of tone change
657 of cymbal depending on grip tightness,” <i>Proceedings of Musical
658 Acoustics, Acoustical Society of Japan</i> <b>37</b>(6), 71–74 (2018, in
659 Japanese).
660 </li>
661 </ol>
662 <button className="primary" onClick={onClose}>← Back to the simulator</button>
663 </section>
664 </div>
665 </div>
666 );
667}