Time-domain nonlinear shell cymbal after Samejima (2021)
45 changed files+8265−0
.github/workflows/deploy.ymladded+70−0View file
@@ -0,0 +1,70 @@
1+name: Deploy to Pages
2+
3+# Builds webapp/ and publishes it as this repository's site.
4+#
5+# Two notes on the vault:
6+# - there is no cache service, so setup-node runs without `cache: npm`;
7+# - a mochi site is served at /<collection>/<repo>/site/ or at the root of an
8+# origin of its own. The Vite build uses a relative base (`base: './'`), so
9+# it works at either without being told the deploy path, and configure-pages
10+# is therefore only here for parity with the deploy step.
11+
12+on:
13+ push:
14+ branches: [main]
15+ workflow_dispatch:
16+
17+permissions:
18+ contents: read
19+ pages: write
20+ id-token: write
21+
22+concurrency:
23+ group: pages
24+ cancel-in-progress: false
25+
26+jobs:
27+ build:
28+ runs-on: ubuntu-latest
29+ steps:
30+ - uses: actions/checkout@v4
31+
32+ - uses: actions/setup-node@v4
33+ with:
34+ node-version: 22
35+
36+ - name: Install
37+ working-directory: webapp
38+ run: npm ci
39+
40+ # Runs matlab/cymbal2_validate.m through numbl: the flat-plate and rigid-body
41+ # closed forms, the paper's linear-FEM reference, the stability of every
42+ # azimuthal column of the update, the stick statics, reproducibility
43+ # of the render, and the pinned end-to-end diagnostics. Fails the build on any
44+ # mismatch.
45+ - name: Validate the model
46+ working-directory: webapp
47+ run: npm run selftest
48+
49+ - name: Build
50+ working-directory: webapp
51+ run: npm run build
52+
53+ - name: Setup Pages
54+ uses: actions/configure-pages@v5
55+
56+ - name: Upload artifact
57+ uses: actions/upload-pages-artifact@v3
58+ with:
59+ path: webapp/dist
60+
61+ deploy:
62+ needs: build
63+ runs-on: ubuntu-latest
64+ environment:
65+ name: github-pages
66+ url: ${{ steps.deployment.outputs.page_url }}
67+ steps:
68+ - name: Deploy to Pages
69+ id: deployment
70+ uses: actions/deploy-pages@v4
README.mdadded+52−0View file
@@ -0,0 +1,52 @@
1+# Cymbal2
2+
3+A time-domain physical simulation of a struck cymbal, after Samejima (2021), "Nonlinear physical modeling sound synthesis of cymbals involving dynamics of washers and sticks/mallets" (Acoust. Sci. & Tech. 42, 6). The body is a shallow spherical shell governed by the von Karman equations, so the crash is not modeled but computed: the geometric nonlinearity that couples bending to mid-surface stretching is integrated directly, and the slow migration of energy into the high frequencies after a hard stroke, the pitch glides, and the noise-like bifurcations of a crash emerge from the PDE. The support is not a rigid clamp but a washer of finite impedance, entering as the shear boundary condition at the mounting circle, which is what lets the whole shell bounce and tilt on its felt the way a mounted cymbal visibly does. The striking mechanism is carried in full: a mallet head with a Hertz-type nonlinear contact, on a stick modeled as an elastic beam, held by a hand that grips either rigidly or through a rotational spring.
4+
5+The sibling [cymbal](../cymbal) repository models the same instrument by an entirely different route: a modal eigensolve of a profiled shell of revolution, with the nonlinearity carried statistically as a cascade of modal energies and the sound radiated to a listener point. The two models make opposite trade-offs. The modal cymbal renders in seconds, carries the hammered irregularity and a radiation model, but replaces the true nonlinear dynamics with an energy random walk that discards phase. Cymbal2 integrates the true dynamics, so intermodulation, pitch glide, and the attack's nonlinear character are real, the trade-off being that a render takes minutes, the plate is ideally axisymmetric, and the output is the shell's motion at a point rather than a pressure at a listener.
6+
7+## Contents
8+
9+| Path | What it is |
10+|---|---|
11+| `matlab/` | The algorithm in MATLAB syntax, run by [numbl](https://numbl.org) |
12+| `webapp/` | Interactive browser app that executes the MATLAB code client-side |
13+
14+The MATLAB code runs from the command line with `cd matlab && numbl run cymbal2_demo.m`, or `cymbal2_validate.m` to run the checks. The web app is a Vite project; see `webapp/README.md`.
15+
16+## The model
17+
18+The signal chain, with each stage in its own `.m` file:
19+
20+```
21+polar grid, azimuthal transform cymbal2_grid.m
22+ -> radial operators per azimuthal wavenumber cymbal2_perm.m
23+ -> implicit update + boundary rows, inverted cymbal2_system.m
24+stick beam FEM, grip, head mass cymbal2_mallet.m
25+ -> hand stroke -> Hertz contact -> von Karman
26+ shell integration, step by step cymbal2_sim.m
27+ -> pickup acceleration -> decimation -> PCM cymbal2_output.m
28+```
29+
30+Five points are worth stating plainly, because they are where the model commits itself.
31+
32+**The shell is a uniform shallow spherical cap.** A real cymbal has a raised bell and a tapered bow whose curvatures differ; the paper (following Bilbao 2010) replaces both by one spherical cap whose rise sets a single radius of curvature, here 66 cm from a 4 cm rise over a 23 cm radius, with uniform 1 mm thickness. That is the price of the von Karman formulation, whose stress-function form is what makes the nonlinear terms computable: the displacement u and an Airy stress function couple through the bilinear operator L, the stretching enters as one Laplacian term in 1/Rs, and both elliptic solves separate over azimuthal wavenumber. The consequence is audible and known: the linear spectrum is that of an idealized cap, not of a hammered instrument, and there is no lathing or hammering irregularity anywhere in the model. The sibling cymbal repository is the converse commitment.
33+
34+**The support is an impedance, not a clamp.** The inner circle at the washer radius carries a moment-free condition and a shear balance against a mass-damper-spring per unit length of the mounting circumference (the paper's Eq. 15), with the measured washer stiffness 1.65e5 N/m per meter. Solved to convergence, this places the bounce of the whole shell on its felt near 13 Hz and the tilt near 4 Hz, both of which this discretization reproduces against the closed-form rigid-body values in the validation. These wobble modes are subsonic, but they matter twice: they set the boundary condition the audible modes actually see (the paper's point, visible in its Fig. 8 as a shifted low-frequency response), and under a hard stroke the shell swings away from the stick by centimeters, which is what makes the contact of a follow-through stroke a shove rather than a strike. `support_clamped = 1` restores the rigid clamp of the earlier literature and is worth hearing once: the low end stiffens and the wobble disappears.
35+
36+**The nonlinearity is integrated explicitly, and that sets a hardness ceiling.** Following the paper's scheme (Eq. 20), the stiff linear terms are implicit with the free parameter g = 0.2498, unconditionally stable on their own, while L(phi, u) and the membrane coupling are evaluated at the current step. An explicit quadratic term is only conditionally stable, and the condition tightens with amplitude: strokes whose tip arrives above roughly 10 m/s drive peak displacements past forty plate thicknesses and the integration diverges within a few hundred milliseconds, at any rate we tried. This is also the regime where the shallow-shell expansion itself has stopped being valid, so the model refuses more or less where its physics runs out; the renderer detects the divergence and reports it rather than emitting garbage. The practical range of `strike_velocity` is about 0.25 to 7 m/s of hand speed, which spans a fingertip tap to a violent crash. Two discretization choices govern the headroom, and not the ones intuition suggests. The azimuthal resolution matters more than the time step: the cascade carries energy up the nodal-diameter ladder, and if `n_theta` truncates that ladder short, energy piles up at the grid's cutoff and the integration blows up where the full grid rides it out, so a cheap render should drop `oversample` and never `n_theta`. The simulation runs at `oversample` times the output rate; 2 matches the paper's 2x rate and adds margin, but at the paper's grid even 1 carries the default stroke.
37+
38+**The striking mechanism is the paper's, with three repairs it needs in practice.** The stick is a 39 cm Euler-Bernoulli beam in cubic Hermite elements with consistent mass, damped in proportion to its mass, gripped at one end either rigidly or against a 140 N m/rad rotational spring, with the 28 g head at the tip and the Hertz-type contact f = K delta^2.54 against the shell. First, the paper's separate head SDOF coupled to the beam by a rigid constraint (its Eqs. 31-39) mixes an explicit and an implicit update through that constraint and the discrete pair is unstable; since the connection is rigid, the head mass is lumped onto the tip node, which is the same physics under substitution and inherits the beam scheme's stability. Second, the implicit weighting printed for the washer and beam schemes (the (1-phi)/2 split with phi = 0.7) is unstable whenever the mass term does not dominate, massless washer included; the weight belongs on the (k+1, k-1) average, and with it phi = 0.7 is unconditionally stable. Third, the paper's hand accelerates to a terminal speed and never reverses, so it cannot play a second stroke and eventually drags the stick through geometry no linear beam describes; here the hand accelerates at an arm-like 250 m/s^2, coasts, and on contact snaps back at 2500 m/s^2, returning to rest so a figure can repeat. Because the stick is compliant and the grip rotates, the tip arrives faster than the hand by a factor that settles around 1.3 when the clearance is sized to the stroke, which `strike_gap_m <= 0` (the default) does automatically. The multiple contacts of a real mallet bounce, which the paper shows as double force peaks, appear here the same way.
39+
40+**The output is the shell's motion at a point, not a microphone.** As in the paper, the rendered signal is the normal acceleration of the plate at a pickup point (default at 0.6 of the radius, opposite the strike), DC-blocked and decimated. This is the model's honest output: a radiation computation like the sibling cymbal's dipole sheet would need the full surface velocity history and is not part of the paper's method. The audible consequence is that low-order modes that radiate poorly from a real plate are present here at full strength, so the result is closer to a contact pickup than to a room microphone.
41+
42+Numerically, the solver departs from the paper in one structural way: instead of assembling the Kronecker-product operators over the full grid and solving one large system per step, both elliptic solves are block-diagonalized over azimuthal wavenumber (every operator and boundary row of the model is axisymmetric, so this is exact), and each step performs two batches of small precomputed radial solves plus dense transforms. The update it computes is the paper's, at about a hundredth the cost, which is the difference between a render taking minutes and taking hours. On grids coarse enough that a wavenumber's inner boundary layer cannot be resolved, the discretized boundary conditions can hand that wavenumber a negative stiffness; the system builder detects such columns by their spectral radius and truncates them, and at the default grid none are.
43+
44+## Validation
45+
46+`matlab/cymbal2_validate.m` (also `npm run selftest` in `webapp/`) checks the flat-plate limit against the closed-form free circular plate for one through four nodal diameters, the washer bounce and tilt against rigid-body closed forms in the stiff-plate limit, the clamped-center eigenfrequencies against the linear-FEM values published in the paper's Fig. 9 for the modes both discretizations resolve the same way, the exact invertibility of the azimuthal transform and the agreement of the physical-space derivative pipeline with the per-wavenumber operators, the non-amplification (spectral radius at most one) of every azimuthal column of the assembled update, the stick's static tip deflection under both grips against beam theory, exact reproducibility of a full render, pinned contact and output diagnostics from a reduced configuration, and that the crash decays. The nonlinear trajectory is sensitive to rounding, so the acoustic pins carry percent-level tolerances where the contact pins are tight.
47+
48+## Defaults
49+
50+The paper's measured instrument: an 18-inch medium crash, 23 cm radius, 4 cm rise, 1 mm bronze at 8860 kg/m^3, on its measured washer stiffness, struck by the paper's mallet (25 mm, 28 g head, on the measured stick) at 0.6 of the radius. The grid is the paper's, 24 radial intervals by 60 azimuthal points, and the damping constants are the paper's fitted sigma0 = 1.34 s^-1 and sigma1 = 1.2e-3, which put T60 near 5 s at the bottom of the spectrum and near 2 s at 5 kHz; the default tail of 6 s carries the ring to the end. The default stroke (5.75 m/s hand, the tip arriving at 8.2 m/s, 311 N peak contact) was chosen against the authors' published demonstration recording: it matches that recording's band balance to within about 2 dB from 100 Hz to 12 kHz, its high-frequency buildup over the first half second, and its decay envelope, the residual being a few dB of extra energy below 700 Hz, where the pickup hears the washer wobble. The same stroke's peak pickup acceleration, 2.4e4 m/s^2, sits at the scale of the paper's Fig. 7 waveforms. A default render integrates about 580,000 steps at 96 kHz and takes a few minutes in MATLAB or under the numbl CLI; `oversample = 1` halves that at some cost in headroom for hard strokes. The mesh resolves the dynamics reliably to roughly 2 kHz by the paper's own estimate; the crash extends far above that, but its top octaves are dispersion-warped and weak, which is the main respect in which this model is duller than a real crash. Raising `n_radial` and `n_theta` together sharpens it at quadratic cost.
51+
52+`strike_velocity` is the one control that matters most: 0.5 m/s is a tap that stays nearly linear and rings like the bell of a gong, 2 m/s begins to wash, and 4 m/s and above crash, with the high-frequency energy arriving over hundreds of milliseconds after the attack. `strike_count` with `strike_rate_hz` plays a repeated figure; each stroke's contact speed varies with `velocity_spread`. The washer parameters, the grip (`grip_clamped`, `grip_spring`), and the stick and head properties are all exposed and all audible in the attack.
export-src/cymbal2/export.jsonadded+23−0View file
@@ -0,0 +1,23 @@
1+{
2+ "id": "cymbal2",
3+ "name": "Cymbal2",
4+ "description": "Time-domain nonlinear shell cymbal (Samejima 2021); pitch scales the shell around MIDI 60 (48..72), velocity onto stroke speed. Renders are slow: the crash is computed, not modeled.",
5+ "source": "simulated_instruments/cymbal2",
6+ "knobs": [
7+ {
8+ "key": "strike_r",
9+ "label": "strike radius",
10+ "min": 0.05,
11+ "max": 0.99,
12+ "default": 0.6
13+ },
14+ {
15+ "key": "grip_spring",
16+ "label": "grip stiffness (N m/rad)",
17+ "min": 20,
18+ "max": 2000,
19+ "default": 140,
20+ "scale": "log"
21+ }
22+ ]
23+}
export-src/cymbal2/render_note.madded+37−0View file
@@ -0,0 +1,37 @@
1+function y = render_note(note)
2+% MIDI pitch scales the shell geometry: every length by
3+% s = 2^-((pitch - 60)/12), which shifts the whole spectrum by the
4+% semitone ratio exactly (shell modes go as thickness / radius^2, and the
5+% curvature and washer frequencies scale the same way when the washer
6+% circle scales too). The mallet stays fixed, so this is one player
7+% striking a family of cymbal sizes; pitch 60 is the default 18-inch
8+% shell, playable 48..72. The sound is inharmonic, so the note number is
9+% a handle on register, not a tuned fundamental. Velocity drives the hand
10+% speed of the stroke, and the crash character follows from the physics:
11+% soft notes stay nearly linear, hard notes wash. Fixed: a single stroke,
12+% a fixed seed, normalize off so velocity survives. The preview tier
13+% coarsens the azimuthal grid and shortens the tail; the full tier
14+% carries the ring to the end and takes minutes per distinct note.
15+P = cymbal2_defaults();
16+P.sample_rate = note.sample_rate;
17+P.strike_count = 1;
18+P.seed = 1;
19+P.pre_delay_ms = 0;
20+P.normalize = 0;
21+% The preview tier drops the oversampling and shortens the tail but keeps
22+% the azimuthal grid: hard strokes need the full wavenumber band to carry
23+% the cascade, and coarsening it is what destabilizes the integration.
24+if note.preview ~= 0
25+ P.oversample = 1;
26+ P.tail_s = 1.5;
27+end
28+s = 2 ^ (-(max(48, min(72, note.pitch)) - 60) / 12);
29+P.radius_m = P.radius_m * s;
30+P.rise_m = P.rise_m * s;
31+P.thickness_m = P.thickness_m * s;
32+P.washer_radius_m = P.washer_radius_m * s;
33+P.strike_velocity = 0.3 + (6.5 - 0.3) * note.velocity;
34+P.strike_r = note.params.strike_r;
35+P.grip_spring = note.params.grip_spring;
36+y = cymbal2_render(P);
37+end
exports/cymbal2.jsonadded+74−0View file
@@ -0,0 +1,74 @@
1+{
2+ "format": "numbl-instrument-v1",
3+ "id": "cymbal2",
4+ "name": "Cymbal2",
5+ "description": "Time-domain nonlinear shell cymbal (Samejima 2021); pitch scales the shell around MIDI 60 (48..72), velocity onto stroke speed. Renders are slow: the crash is computed, not modeled.",
6+ "source": "simulated_instruments/cymbal2",
7+ "knobs": [
8+ {
9+ "key": "strike_r",
10+ "label": "strike radius",
11+ "min": 0.05,
12+ "max": 0.99,
13+ "default": 0.6
14+ },
15+ {
16+ "key": "grip_spring",
17+ "label": "grip stiffness (N m/rad)",
18+ "min": 20,
19+ "max": 2000,
20+ "default": 140,
21+ "scale": "log"
22+ }
23+ ],
24+ "files": [
25+ {
26+ "name": "render_note.m",
27+ "content": "function y = render_note(note)\n% MIDI pitch scales the shell geometry: every length by\n% s = 2^-((pitch - 60)/12), which shifts the whole spectrum by the\n% semitone ratio exactly (shell modes go as thickness / radius^2, and the\n% curvature and washer frequencies scale the same way when the washer\n% circle scales too). The mallet stays fixed, so this is one player\n% striking a family of cymbal sizes; pitch 60 is the default 18-inch\n% shell, playable 48..72. The sound is inharmonic, so the note number is\n% a handle on register, not a tuned fundamental. Velocity drives the hand\n% speed of the stroke, and the crash character follows from the physics:\n% soft notes stay nearly linear, hard notes wash. Fixed: a single stroke,\n% a fixed seed, normalize off so velocity survives. The preview tier\n% coarsens the azimuthal grid and shortens the tail; the full tier\n% carries the ring to the end and takes minutes per distinct note.\nP = cymbal2_defaults();\nP.sample_rate = note.sample_rate;\nP.strike_count = 1;\nP.seed = 1;\nP.pre_delay_ms = 0;\nP.normalize = 0;\n% The preview tier drops the oversampling and shortens the tail but keeps\n% the azimuthal grid: hard strokes need the full wavenumber band to carry\n% the cascade, and coarsening it is what destabilizes the integration.\nif note.preview ~= 0\n P.oversample = 1;\n P.tail_s = 1.5;\nend\ns = 2 ^ (-(max(48, min(72, note.pitch)) - 60) / 12);\nP.radius_m = P.radius_m * s;\nP.rise_m = P.rise_m * s;\nP.thickness_m = P.thickness_m * s;\nP.washer_radius_m = P.washer_radius_m * s;\nP.strike_velocity = 0.3 + (6.5 - 0.3) * note.velocity;\nP.strike_r = note.params.strike_r;\nP.grip_spring = note.params.grip_spring;\ny = cymbal2_render(P);\nend\n"
28+ },
29+ {
30+ "name": "cymbal2_defaults.m",
31+ "content": "function P = cymbal2_defaults()\n% CYMBAL2_DEFAULTS Default parameters, following Samejima (2021) Table 1:\n% an 18-inch medium crash (45 cm diameter) modeled as a shallow spherical\n% shell, supported at the center by a washer of finite impedance and struck\n% by a mallet whose stick is an elastic beam held in a compliant grip.\n% The output duration is not a parameter: it is\n% (strike_count - 1) / strike_rate_hz + tail_s.\n\n% Shell geometry (a shallow spherical cap; rise sets the radius of curvature\n% Rs = radius^2 / (2 * rise))\nP.radius_m = 0.230;\nP.rise_m = 0.040;\nP.thickness_m = 0.001;\n\n% Material (bronze, per the paper's measured cymbal)\nP.youngs_modulus = 1.1e11;\nP.poisson_ratio = 0.33;\nP.density = 8860.0;\n\n% Damping: decay rate sigma0 + sigma1 * a^2 * k^2 at flexural wavenumber k\nP.sigma0 = 1.34; % s^-1, frequency-independent\nP.sigma1 = 1.2e-3; % s^-1 per unit nondimensional k^2\n\n% Support at the center: an SDOF washer per unit length of the inner\n% circumference, entering as the shear boundary condition at r = washer\n% radius. support_clamped = 1 replaces it with the rigid clamp of the\n% earlier literature (Bilbao 2010), which is worth hearing once.\nP.washer_radius_m = 0.029;\nP.washer_stiffness = 1.65e5; % N/m per m\nP.washer_mass = 0.0; % kg per m (paper sets 0)\nP.washer_damping = 0.0; % N s/m per m (paper sets 0)\nP.support_clamped = 0;\n\n% Mallet head (SDOF, Hertz-type nonlinear stiffness f = K * delta^alpha)\nP.head_radius_m = 0.025;\nP.head_mass_kg = 0.028;\nP.head_stiffness = 1.6e8; % N m^-alpha\nP.head_alpha = 2.54;\n\n% Mallet stick (Euler-Bernoulli beam, FEM)\nP.stick_length_m = 0.390;\nP.stick_area_m2 = 1.43139e-4;\nP.stick_youngs = 8.8e9;\nP.stick_density = 400.0;\nP.stick_inertia_m4 = 1.63044e-9;\nP.stick_damping_rate = 10.0; % mass-proportional, R = rate * M\nP.n_beam = 20; % beam elements (21 nodes)\n\n% Grip: 0 = rotational spring hand (loose grip), 1 = clamped hand (tight)\nP.grip_clamped = 0;\nP.grip_spring = 140.0; % N m/rad\n\n% Strike gesture: the hand accelerates downward to strike_velocity, coasts\n% until the head touches, holds briefly, and retracts to its rest height.\n% The paper prescribes constant acceleration to the same terminal speed\n% but never reverses, so it cannot play a second stroke. Because the stick\n% is compliant, the tip arrives somewhat slower than the hand.\nP.strike_velocity = 5.75; % peak hand speed, m/s; reproduces the\n % character of the authors' published demo\nP.strike_accel = 250.0; % downstroke hand acceleration, m/s^2\nP.retract_accel = 2500.0; % wrist-snap retract acceleration, m/s^2\nP.strike_hold_ms = 0.5; % dwell between contact and retract\nP.strike_gap_m = -1; % rest clearance; <= 0 sizes it to the stroke\nP.strike_r = 0.6; % contact radius, 0 at the washer, 1 at the rim\nP.strike_theta_deg = 180.0;\nP.strike_count = 1;\nP.strike_rate_hz = 4.0;\nP.velocity_spread = 0.0; % relative velocity scatter of repeated strokes\nP.seed = 1;\n\n% Pickup: normal acceleration of the shell at one point, as in the paper\nP.pickup_r = 0.6;\nP.pickup_theta_deg = 0.0;\n\n% Discretization\nP.n_radial = 24; % radial grid intervals (the paper's N_r)\nP.n_theta = 60; % azimuthal grid points (the paper's N_theta)\nP.g_param = 0.2498; % free parameter of the implicit shell scheme\nP.phi_param = 0.7; % free parameter of the washer and beam schemes\n\n% Rendering / output\nP.sample_rate = 48000;\nP.oversample = 2; % simulation runs at oversample * sample_rate\nP.tail_s = 6.0; % the low modes carry T60 of about 5 s\nP.pre_delay_ms = 3.0;\nP.normalize = 1;\nP.peak_dbfs = -3.0;\nP.full_scale_accel = 4.0e4; % m/s^2 mapped to full scale when normalize = 0\nend\n"
32+ },
33+ {
34+ "name": "cymbal2_demo.m",
35+ "content": "% Render a cymbal crash with default parameters and print diagnostics.\n% Set P.strike_velocity = 0.5 for a tap on the same instrument, or\n% P.support_clamped = 1 to hear the rigid center clamp of the earlier\n% literature instead of the washer. P.strike_count = 4 plays a repeated\n% figure; the hand retracts between strokes, which the paper's\n% constant-acceleration hand cannot do.\n% Run with: numbl run cymbal2_demo.m\n\nP = cymbal2_defaults();\ntic;\n[y, dinfo] = cymbal2_render(P);\nelapsed = toc;\n\nfprintf('grid : %d radial x %d azimuthal, %.1f kHz simulation rate\\n', ...\n P.n_radial, P.n_theta, dinfo.simRate / 1000);\nfprintf('shell : a = %.0f cm, rise %.0f mm -> Rs = %.2f m, washer at %.1f cm\\n', ...\n P.radius_m * 100, P.rise_m * 1000, dinfo.RsM, P.washer_radius_m * 100);\nf0 = cymbal2_modes(P, 0); f0 = f0(f0 > 1);\nf1 = cymbal2_modes(P, 1); f1 = f1(f1 > 1);\nf2 = cymbal2_modes(P, 2); f2 = f2(f2 > 1);\nfprintf('low modes : %.1f Hz (m=1 tilt), %.1f Hz (m=2), %.1f Hz (m=0 bounce)\\n', ...\n f1(1), f2(1), f0(1));\nfprintf('stroke : hand %.1f m/s -> tip %.2f m/s at contact\\n', ...\n dinfo.meanV, abs(dinfo.contactVel(1)));\nfprintf('contact : %.0f N peak, %.2f mm compression, %.1f ms over %d event(s)\\n', ...\n dinfo.peakForceN, dinfo.peakXiMm, dinfo.contactMs, dinfo.contactEvents);\nfprintf('shell peak : %.2f mm displacement = %.1f thicknesses\\n', ...\n dinfo.peakDispM * 1000, dinfo.peakDispThickness);\nfprintf('pickup : %.2e m/s^2 peak acceleration at r = %.1f cm\\n', ...\n dinfo.peakAccel, dinfo.pickupRadiusM * 100);\nfprintf('length : %.2f s @ %d Hz\\n', dinfo.lengthS, P.sample_rate);\nfprintf('render time: %.1f s\\n', elapsed);\n"
36+ },
37+ {
38+ "name": "cymbal2_grid.m",
39+ "content": "function G = cymbal2_grid(P)\n% CYMBAL2_GRID Polar grid, radial difference matrices, and the azimuthal\n% transform. The radial grid carries l = 1..Nr+2 lines: real nodes at\n% l = 3..Nr spanning [washer radius, rim radius], and two virtual lines\n% beyond each boundary, eliminated by the boundary-condition rows of the\n% assembled systems (Samejima 2021, Fig. 4). Azimuthal operators are\n% circulant, so every linear solve block-diagonalizes over the real\n% cos/sin basis; F maps physical theta samples to coefficients (U * F)\n% and Gi maps back (Uhat * Gi), with mcol giving each column's wavenumber.\n\na = P.radius_m;\nc = P.washer_radius_m;\nNr = P.n_radial;\nnth = P.n_theta;\nif mod(nth, 2) ~= 0\n error('cymbal2: n_theta must be even');\nend\nif Nr < 8\n error('cymbal2: n_radial must be at least 8');\nend\nnrt = Nr + 2;\nhr = (a - c) / (Nr - 3);\nif c <= hr\n error('cymbal2: washer_radius too small for n_radial (need c > (a-c)/(n_radial-3))');\nend\n\nr = c + ((1:nrt)' - 3) * hr; % r(1) = c - 2hr may be tiny on coarse grids;\nr(1) = max(r(1), hr / 2); % row 1 of every operator is never used\n\n% Radial stencils (dense; nrt is small). Rows outside the valid range are\n% left zero and are overwritten by boundary-condition rows downstream.\nD1 = zeros(nrt, nrt);\nD2 = zeros(nrt, nrt);\nD3 = zeros(nrt, nrt);\nfor l = 2:nrt - 1\n D1(l, l - 1) = -1 / (2 * hr); D1(l, l + 1) = 1 / (2 * hr);\n D2(l, l - 1) = 1 / hr^2; D2(l, l) = -2 / hr^2; D2(l, l + 1) = 1 / hr^2;\nend\nfor l = 3:nrt - 2\n D3(l, l - 2) = -1 / (2 * hr^3);\n D3(l, l - 1) = 2 / (2 * hr^3);\n D3(l, l + 1) = -2 / (2 * hr^3);\n D3(l, l + 2) = 1 / (2 * hr^3);\nend\n\n% Azimuthal stencils as dense circulant matrices, applied on the right\n% (U * T2 etc., so they carry a transpose; T2 is symmetric, T1 is\n% antisymmetric and T1t below is arranged so Uth = U * T1t).\nhth = 2 * pi / nth;\nT2 = zeros(nth, nth);\nT1t = zeros(nth, nth);\nfor j = 1:nth\n jm = 1 + mod(j - 2, nth);\n jp = 1 + mod(j, nth);\n T2(j, j) = -2 / hth^2; T2(j, jm) = 1 / hth^2; T2(j, jp) = 1 / hth^2;\n % row j of T1 has -1 at jm, +1 at jp (scaled); transpose: column j\n T1t(jm, j) = -1 / (2 * hth);\n T1t(jp, j) = 1 / (2 * hth);\nend\n\n% Real transform: columns ordered m = 0, (1c,1s), (2c,2s), ..., nth/2\nth = ((1:nth)' - 1) * hth;\nF = zeros(nth, nth);\nGi = zeros(nth, nth);\nmcol = zeros(1, nth);\nF(:, 1) = 1 / nth;\nGi(1, :) = 1;\ncol = 2;\nfor m = 1:nth / 2 - 1\n F(:, col) = 2 * cos(m * th) / nth; Gi(col, :) = cos(m * th)'; mcol(col) = m; col = col + 1;\n F(:, col) = 2 * sin(m * th) / nth; Gi(col, :) = sin(m * th)'; mcol(col) = m; col = col + 1;\nend\nF(:, col) = cos((nth / 2) * th) / nth;\nGi(col, :) = cos((nth / 2) * th)';\nmcol(col) = nth / 2;\n\n% Eigenvalue of the circulant second-difference for wavenumber m\nlam2 = (2 * cos(mcol * hth) - 2) / hth^2;\n\nG.a = a; G.c = c; G.Nr = Nr; G.nrt = nrt; G.nth = nth;\nG.hr = hr; G.hth = hth;\nG.r = r; G.th = th;\nG.D1 = D1; G.D2 = D2; G.D3 = D3;\nG.T1t = T1t; G.T2 = T2;\nG.F = F; G.Gi = Gi; G.mcol = mcol; G.lam2 = lam2;\nG.Rd = 1.0 ./ r; % 1/r column, broadcast over theta\nG.R2d = 1.0 ./ r.^2;\nend\n"
40+ },
41+ {
42+ "name": "cymbal2_mallet.m",
43+ "content": "function MB = cymbal2_mallet(P)\n% CYMBAL2_MALLET Precompute the striking mechanism: the mallet stick as an\n% Euler-Bernoulli beam (Hermite cubic elements, consistent mass,\n% mass-proportional damping), held at one end by the hand, with the mallet\n% head attached at the other end as a separate SDOF body (Samejima 2021,\n% Sec. 3). The hand grips either rigidly (grip_clamped = 1: transverse\n% displacement and rotation both prescribed) or through a rotational\n% spring (grip_clamped = 0: displacement prescribed, rotation free against\n% grip_spring). The head couples to the beam tip by the rigid connection\n% u_M = u_S. The paper integrates the head as its own SDOF (Eq. 34) and\n% solves the connection force from the constraint; because that couples an\n% explicit head update to the implicit beam, the discrete pair is\n% unstable, so the head mass is lumped onto the tip node instead, which is\n% the same physics under the rigid connection (substituting u_M = u_S and\n% f_S = -f_M2 into Eq. (31) puts M_M on the tip and f_M1 on the beam).\n\nne = P.n_beam;\nnn = ne + 1;\nndof = 2 * nn;\nle = P.stick_length_m / ne;\nEI = P.stick_youngs * P.stick_inertia_m4;\nrA = P.stick_density * P.stick_area_m2;\n\nKe = EI / le^3 * [12, 6 * le, -12, 6 * le; ...\n 6 * le, 4 * le^2, -6 * le, 2 * le^2; ...\n -12, -6 * le, 12, -6 * le; ...\n 6 * le, 2 * le^2, -6 * le, 4 * le^2];\nMe = rA * le / 420 * [156, 22 * le, 54, -13 * le; ...\n 22 * le, 4 * le^2, 13 * le, -3 * le^2; ...\n 54, 13 * le, 156, -22 * le; ...\n -13 * le, -3 * le^2, -22 * le, 4 * le^2];\n\nK = zeros(ndof, ndof);\nM = zeros(ndof, ndof);\nfor e = 1:ne\n ix = (2 * e - 1):(2 * e + 2);\n K(ix, ix) = K(ix, ix) + Ke;\n M(ix, ix) = M(ix, ix) + Me;\nend\nR = P.stick_damping_rate * M;\nM(ndof - 1, ndof - 1) = M(ndof - 1, ndof - 1) + P.head_mass_kg;\n\n% Node 1 is the grip, node nn the tip. Dof 1 (grip displacement) is always\n% prescribed to the hand motion; dof 2 (grip rotation) is prescribed to 0\n% for the clamped hand, or free with the rotational spring added.\nif P.grip_clamped ~= 0\n fi = 3:ndof;\nelse\n K(2, 2) = K(2, 2) + P.grip_spring;\n fi = 2:ndof;\nend\nKf = K(fi, fi); Mf = M(fi, fi); Rf = R(fi, fi);\nKc = K(fi, 1); Mc = M(fi, 1); Rc = R(fi, 1);\n\nht = 1 / (P.sample_rate * P.oversample);\nps = P.phi_param;\n% Stiffness split with the implicit weight ps/2 on the (k+1, k-1) average,\n% as in the shell's washer row; the complementary split printed in the\n% paper is unstable for the beam's highest elements at audio rates.\nAS = Mf + (ht / 2) * Rf + (ps / 2) * ht^2 * Kf;\nBS = 2 * Mf - (1 - ps) * ht^2 * Kf;\nCS = -Mf + (ht / 2) * Rf - (ps / 2) * ht^2 * Kf;\nASi = inv(AS);\nMB.ASiB = ASi * BS;\nMB.ASiC = ASi * CS;\n\n% Hand coupling, folded so the per-step load is\n% vM * hA(k) + vR * hB(k) + vK * hC(k)\n% with hA = -(uH+ - 2 uH + uH-), hB = -(ht/2)(uH+ - uH-),\n% hC = -ht^2 (ps uH + ((1-ps)/2)(uH+ + uH-)).\nMB.vM = ASi * Mc;\nMB.vR = ASi * Rc;\nMB.vK = ASi * Kc;\n\ntipi = find(fi == ndof - 1);\nMB.tipi = tipi;\nMB.tipcol = ASi(:, tipi);\nMB.gtip = ASi(tipi, tipi);\nMB.nfree = length(fi);\nMB.ht = ht;\nMB.ps = ps;\nend\n"
44+ },
45+ {
46+ "name": "cymbal2_modes.m",
47+ "content": "function [f, W] = cymbal2_modes(P, m)\n% CYMBAL2_MODES Eigenfrequencies (Hz) of the linearized shell for one\n% azimuthal wavenumber m, with the membrane (curvature) stiffening and the\n% support boundary condition included; the washer is taken as its static\n% stiffness (its mass and damping are ignored here, as in the paper's\n% linear FEM reference). Used by the demo and the validation, not by the\n% renderer. W returns the mode shapes on the real radial nodes l = 3..Nr.\n\nG = cymbal2_grid(P);\nnrt = G.nrt;\nNr = G.Nr;\nhth = G.hth;\nlam2 = (2 * cos(m * hth) - 2) / hth^2;\nop = cymbal2_perm(P, G, lam2);\n\nE = P.youngs_modulus;\nh = P.thickness_m;\nrhoh = P.density * h;\nRs = P.radius_m^2 / (2 * P.rise_m);\n\n% Airy solve operator: Fm * phi = Z * (E h / Rs) * Lap * u, with\n% phi = phi_r = 0 at both circles and identity rows on the outermost\n% virtual lines.\nFm = zeros(nrt, nrt);\nFm(1, 1) = 1;\nFm(2, :) = G.D1(3, :);\nFm(3, 3) = 1;\nfor l = 4:Nr - 1\n Fm(l, :) = op.Bih(l, :);\nend\nFm(Nr, Nr) = 1;\nFm(Nr + 1, :) = G.D1(Nr, :);\nFm(Nr + 2, Nr + 2) = 1;\nZl = zeros(nrt, nrt);\nfor l = 4:Nr - 1\n Zl(l, :) = op.Lap(l, :);\nend\nPhimap = Fm \\ ((E * h / Rs) * Zl);\n\n% Full stiffness on the equation-of-motion rows l = 3..Nr\nKfull = zeros(nrt, nrt);\nKterm = op.Dflex * op.Bih + (1 / Rs) * op.Lap * Phimap;\nfor l = 3:Nr\n Kfull(l, :) = Kterm(l, :);\nend\n\n% Boundary rows, static. With the washer the shell's equation of motion\n% holds on l = 3..Nr and the four virtual lines are eliminated; with the\n% rigid clamp u and u_r vanish at l = 3, so u(3) itself joins the\n% eliminated set and the equations of motion start at l = 4 (the innermost\n% virtual line is then never engaged).\ne3 = zeros(1, nrt); e3(3) = 1;\nif P.support_clamped ~= 0\n Cmat = [e3; G.D1(3, :); op.MomE; op.ShE];\n keep = 4:Nr;\n elim = [2, 3, Nr + 1, Nr + 2];\nelse\n Cmat = [op.MomW; op.ShW - P.washer_stiffness * e3; op.MomE; op.ShE];\n keep = 3:Nr;\n elim = [1, 2, Nr + 1, Nr + 2];\nend\nX = -Cmat(:, elim) \\ Cmat(:, keep);\nKred = Kfull(keep, keep) + Kfull(keep, elim) * X;\n\n[V, Dm] = eig(Kred / rhoh);\nom2 = real(diag(Dm));\n[om2, idx] = sort(om2);\nom2 = max(om2, 0);\nf = sqrt(om2) / (2 * pi);\nW = real(V(:, idx));\nend\n"
48+ },
49+ {
50+ "name": "cymbal2_output.m",
51+ "content": "function [y, peakAccel] = cymbal2_output(P, acc)\n% CYMBAL2_OUTPUT Output conditioning: DC blocker on the raw pickup\n% acceleration, Kaiser-windowed-sinc anti-alias decimation from the\n% simulation rate to the output rate, tail fade, and normalization. Same\n% structure as the sibling instruments' output stages.\n\npeakAccel = max(abs(acc));\nR = P.oversample;\nT = length(acc);\n\n% One-pole DC blocker (the pickup is an acceleration, but the contact\n% leaves a small quasi-static offset while the washer relaxes)\nfc = 15.0;\na1 = exp(-2 * pi * fc / (P.sample_rate * R));\np = zeros(1, T);\nxp = 0; yp = 0;\nfor k = 1:T\n yv = acc(k) - xp + a1 * yp;\n xp = acc(k);\n yp = yv;\n p(k) = yv;\nend\n\nif R > 1\n taps = 12 * R + 1;\n M2 = (taps - 1) / 2;\n fcNorm = 0.45 / R;\n m = (0:taps - 1) - M2;\n hh = sin(2 * pi * fcNorm * m) ./ (pi * m);\n hh(M2 + 1) = 2 * fcNorm;\n r = m / M2;\n beta = 8;\n hh = hh .* (c2_i0(beta * sqrt(max(0, 1 - r.^2))) / c2_i0(beta));\n hh = hh / sum(hh);\n pf = conv(p, hh);\n outLen = floor(T / R);\n y = reshape(pf((0:outLen - 1) * R + M2 + 1), 1, outLen);\nelse\n y = p;\nend\n\n% Squared fade on the tail\nfade = min(2048, floor(length(y) / 8));\nif fade > 1\n y(end - fade + 1:end) = y(end - fade + 1:end) .* linspace(1, 0, fade).^2;\nend\n\nif P.normalize ~= 0\n pk = max(abs(y));\n if pk > 0\n y = y / pk * 10^(P.peak_dbfs / 20);\n end\nelse\n y = y / P.full_scale_accel;\nend\nend\n\nfunction v = c2_i0(x)\n% Modified Bessel function of the first kind, order zero (series)\nv = ones(size(x));\nterm = ones(size(x));\nhx = x / 2;\nfor k = 1:40\n term = term .* (hx / k).^2;\n v = v + term;\nend\nend\n"
52+ },
53+ {
54+ "name": "cymbal2_perm.m",
55+ "content": "function op = cymbal2_perm(P, G, lam2)\n% CYMBAL2_PERM Radial operators for one azimuthal wavenumber, dimensional.\n% lam2 is the (negative) eigenvalue of the circulant second-difference in\n% theta for this wavenumber, so every d^2/dtheta^2 becomes multiplication\n% by lam2. Returns the Laplacian, the biharmonic (as the exact discrete\n% composition Lap * Lap), and the radial-moment and Kelvin-Kirchhoff shear\n% rows used as boundary conditions at the washer (l = 3) and at the free\n% edge (l = Nr). The moment and shear expressions are those of Samejima\n% (2021), Eqs. (6)-(7) and (14)-(15).\n\nnrt = G.nrt;\nnu = P.poisson_ratio;\nDfl = P.youngs_modulus * P.thickness_m^3 / (12 * (1 - nu^2));\n\nRm = diag(G.Rd);\nR2 = diag(G.R2d);\nR3 = diag(G.Rd.^3);\n\nLap = G.D2 + Rm * G.D1 + lam2 * R2;\nBih = Lap * Lap;\n\n% Radial bending moment (per unit length, divided by -D):\n% u_rr + nu * (u_r / r + u_thth / r^2)\nMom = G.D2 + nu * (Rm * G.D1 + lam2 * R2);\n\n% Kelvin-Kirchhoff shear: V = -D * [u_rrr + u_rr / r - u_r / r^2\n% + (nu - 3) / r^3 * u_thth + (2 - nu) / r^2 * u_r_thth]\nSh = -Dfl * (G.D3 + Rm * G.D2 - R2 * G.D1 ...\n + (nu - 3) * lam2 * R3 + (2 - nu) * lam2 * R2 * G.D1);\n\nop.Lap = Lap;\nop.Bih = Bih;\nop.MomW = Mom(3, :); % moment row at the washer circle\nop.MomE = Mom(nrt - 2, :); % moment row at the free edge (l = Nr)\nop.ShW = Sh(3, :);\nop.ShE = Sh(nrt - 2, :);\nop.Dflex = Dfl;\nend\n"
56+ },
57+ {
58+ "name": "cymbal2_render.m",
59+ "content": "function [y, dinfo] = cymbal2_render(P)\n% CYMBAL2_RENDER Full render: precompute -> time integration -> output.\n\nS = cymbal2_system(P);\nMB = cymbal2_mallet(P);\n[acc, dinfo] = cymbal2_sim(P, S, MB);\n[y, peakAccel] = cymbal2_output(P, acc);\n\ndinfo.peakAccel = peakAccel;\ndinfo.strikeRadiusM = S.strikeRadiusM;\ndinfo.pickupRadiusM = S.pickupRadiusM;\ndinfo.RsM = S.Rs;\ndinfo.lengthS = length(y) / P.sample_rate;\nend\n"
60+ },
61+ {
62+ "name": "cymbal2_sim.m",
63+ "content": "function [acc, dinfo, fhist] = cymbal2_sim(P, S, MB)\n% CYMBAL2_SIM Time-domain integration of the full system: von Karman\n% shallow spherical shell (implicit in the stiff linear terms, explicit in\n% the nonlinearity), washer support built into the boundary rows, mallet\n% head + stick coupled through the Hertz contact. Returns the normal\n% acceleration of the shell at the pickup node at the simulation rate,\n% diagnostics, and the contact force history.\n\nG = S.G;\nnrt = G.nrt;\nNr = G.Nr;\nnth = G.nth;\nht = S.ht;\nrhoh = S.rhoh;\nRs = S.Rs;\nEhRs = S.Eh / Rs;\nEh2 = S.Eh / 2;\nls = S.ls; js = S.js; lp = S.lp; jp = S.jp;\nfsc = S.forceScale;\nnlr0 = S.eomRow0;\n\nK = max(1, round(P.strike_count));\n\n% Stroke schedule; per-stroke peak hand speeds with a seeded spread. Each\n% stroke is launched early by its predicted hand travel time (acceleration\n% plus coast over the rest clearance), so the contacts land near the\n% nominal beats instead of drifting with velocity. A nonpositive\n% strike_gap_m asks for the clearance to be sized to the stroke: contact\n% arrives just after the hand reaches full speed, which keeps the tip's\n% whip phase, and with it the contact speed, consistent across strokes.\ntstr = P.pre_delay_ms * 1e-3 + (0:K - 1) / max(1e-3, P.strike_rate_hz);\nvstr = zeros(1, K);\nstate = P.seed * 2654435761;\nfor i = 1:K\n state = mod(state * 1664525 + 1013904223, 4294967296);\n x = state / 4294967296;\n vstr(i) = max(0.05, P.strike_velocity * (1 + P.velocity_spread * (2 * x - 1)));\nend\ngap = P.strike_gap_m;\nif gap <= 0\n gap = 1.05 * max(vstr)^2 / (2 * P.strike_accel) + 0.005;\nend\ntravel = gap ./ vstr + vstr / (2 * P.strike_accel);\ntlaunch = tstr - travel;\noff = -min(0, min(tlaunch));\ntlaunch = tlaunch + off;\n\nduration = off + (K - 1) / max(1e-3, P.strike_rate_hz) + P.tail_s;\nnsteps = floor(duration * P.sample_rate * P.oversample);\n\nKM = P.head_stiffness;\nalM = P.head_alpha;\ntipi = MB.tipi;\nps = MB.ps;\n\n% Hand state machine: the hand accelerates downward, coasts at the stroke\n% speed until the head makes contact (or a miss depth), holds briefly,\n% then retracts at half speed and brakes so it arrives back at rest,\n% ready for the next scheduled stroke. This replaces the paper's\n% constant-acceleration hand, which never reverses. The machine runs one\n% step ahead because the beam coupling needs u_H at k + 1.\naH = P.strike_accel;\naHr = P.retract_accel;\nholdSteps = max(1, round(P.strike_hold_ms * 1e-3 / ht));\nmissDepth = gap + 0.06;\nuHm = 0; uH0 = 0; uHp = 0;\nvH = 0;\nphase = 0;\nstrokeIdx = 1;\nholdLeft = 0;\nvtar = 0;\n\nD1p = G.D1; D2p = G.D2;\nT1t = G.T1t; T2 = G.T2;\nF = G.F; Gi = G.Gi;\nRd = G.Rd; R2d = G.R2d;\nAiB = S.AiB; AiC = S.AiC; Ai = S.Ai; Fi = S.Fi;\nsNL = ht^2 / rhoh;\n\nU = zeros(nrt, nth);\nU1 = zeros(nrt, nth);\nUh3 = zeros(nrt, 1, nth);\nUh13 = zeros(nrt, 1, nth);\nd = zeros(MB.nfree, 1);\nd1 = zeros(MB.nfree, 1);\n\nacc = zeros(1, nsteps);\nfhist = zeros(1, nsteps);\nshellActive = 0;\npeakDisp = 0;\npeakXi = 0;\ninContact = 0;\ncvel = zeros(1, K);\n\nfor k = 1:nsteps\n % ---- mallet: contact force from the state at step k ----\n xi = U(ls, js) - d(tipi) - gap;\n if xi > 0\n fM1 = KM * xi^alM;\n shellActive = 1;\n if inContact == 0 && strokeIdx <= K && cvel(strokeIdx) == 0\n cvel(strokeIdx) = (d(tipi) - d1(tipi)) / ht; % tip speed at contact\n end\n inContact = 1;\n if xi > peakXi\n peakXi = xi;\n end\n else\n fM1 = 0;\n inContact = 0;\n end\n fhist(k) = fM1;\n hA = -(uHp - 2 * uH0 + uHm);\n hB = -(ht / 2) * (uHp - uHm);\n hC = -ht^2 * (ps * uH0 + ((1 - ps) / 2) * (uHp + uHm));\n dnew = MB.ASiB * d + MB.ASiC * d1 + MB.vM * hA + MB.vR * hB + MB.vK * hC ...\n + (ht^2 * fM1) * MB.tipcol;\n\n % ---- advance the hand one step ----\n if phase == 0\n if strokeIdx <= K && k * ht >= tlaunch(strokeIdx)\n phase = 1;\n vtar = vstr(strokeIdx);\n end\n elseif phase == 1\n vH = max(vH - aH * ht, -vtar);\n if inContact ~= 0\n phase = 3;\n holdLeft = holdSteps;\n elseif vH <= -vtar\n phase = 2;\n end\n elseif phase == 2\n if inContact ~= 0\n phase = 3;\n holdLeft = holdSteps;\n elseif uHp <= -missDepth\n phase = 4;\n end\n elseif phase == 3\n holdLeft = holdLeft - 1;\n if holdLeft <= 0\n phase = 4;\n end\n elseif phase == 4\n vH = min(vH + aHr * ht, 0.5 * vtar);\n if vH > 0 && uHp >= -vH^2 / (2 * aHr)\n phase = 5;\n end\n else\n vH = vH - aHr * ht;\n if vH <= 0 || uHp >= 0\n vH = 0;\n phase = 0;\n strokeIdx = strokeIdx + 1;\n end\n end\n uHm = uH0;\n uH0 = uHp;\n uHp = uHp + vH * ht;\n\n if shellActive ~= 0\n % ---- shell derivatives and von Karman brackets at step k ----\n Ur = D1p * U;\n Urr = D2p * U;\n v1 = Rd .* Ur + R2d .* (U * T2);\n v2 = Rd .* (Ur * T1t) - R2d .* (U * T1t);\n LapU = Urr + v1;\n rhsA = EhRs * LapU - Eh2 * 2 * (Urr .* v1 - v2 .* v2);\n rhsA(1:3, :) = 0;\n rhsA(Nr:nrt, :) = 0;\n Ph3 = pagemtimes(Fi, reshape(rhsA * F, nrt, 1, nth));\n Phi = reshape(Ph3, nrt, nth) * Gi;\n Pr = D1p * Phi;\n Prr = D2p * Phi;\n w1 = Rd .* Pr + R2d .* (Phi * T2);\n w2 = Rd .* (Pr * T1t) - R2d .* (Phi * T1t);\n NL = Prr .* v1 + Urr .* w1 - 2 * w2 .* v2 - (Prr + w1) / Rs;\n NL(ls, js) = NL(ls, js) - fM1 * fsc;\n NL(1:nlr0 - 1, :) = 0;\n NL(Nr + 1:nrt, :) = 0;\n NLh = pagemtimes(Ai, reshape((sNL * NL) * F, nrt, 1, nth));\n Uh3n = pagemtimes(AiB, Uh3) + pagemtimes(AiC, Uh13) + NLh;\n Unew = reshape(Uh3n, nrt, nth) * Gi;\n\n acc(k) = (Unew(lp, jp) - 2 * U(lp, jp) + U1(lp, jp)) / ht^2;\n U1 = U;\n U = Unew;\n Uh13 = Uh3;\n Uh3 = Uh3n;\n if mod(k, 512) == 0\n pd = max(abs(U(:, js)));\n if pd > peakDisp\n peakDisp = pd;\n end\n if ~isfinite(pd)\n error('cymbal2: simulation diverged at t = %.4f s (reduce strike_velocity, or raise n_theta / oversample)', k * ht);\n end\n end\n end\n\n d1 = d;\n d = dnew;\nend\n\n% Contact bookkeeping from the force history\ninC = fhist > 0;\ndinfo.peakForceN = max(fhist);\ndinfo.contactSteps = sum(inC);\ndinfo.contactMs = 1000 * ht * dinfo.contactSteps;\ndinfo.contactEvents = sum(diff([0, inC]) > 0);\ndinfo.peakXiMm = peakXi * 1000;\ndinfo.peakDispM = peakDisp;\ndinfo.peakDispThickness = peakDisp / P.thickness_m;\ndinfo.peakAccel = max(abs(acc));\ndinfo.nsteps = nsteps;\ndinfo.simRate = P.sample_rate * P.oversample;\ndinfo.tstr = tstr;\ndinfo.vstr = vstr;\ndinfo.meanV = mean(vstr);\ndinfo.contactVel = cvel; % tip speed at first contact of each stroke\ndinfo.strokesPlayed = min(strokeIdx - 1 + (phase > 0), K);\nend\n"
64+ },
65+ {
66+ "name": "cymbal2_system.m",
67+ "content": "function S = cymbal2_system(P)\n% CYMBAL2_SYSTEM Precompute everything the time loop needs. The linear\n% updates block-diagonalize over azimuthal wavenumber, so the implicit\n% shell solve and the Airy (stress function) solve are stored as one small\n% radial matrix per transform column: the time step becomes\n% Uhat+ = (A^-1 B) Uhat + (A^-1 C) Uhat- + A^-1 (ht^2/rho h) NLhat\n% with A, B, C the matrices of Samejima (2021) Eq. (20) restricted to one\n% wavenumber, the boundary-condition rows of Eqs. (22)-(30) built in, and\n% NL the explicit nonlinear + forcing terms evaluated in physical space.\n% Rows are equilibrated before inversion; the stored inverses absorb the\n% scaling so the loop applies them to raw right-hand sides.\n\nG = cymbal2_grid(P);\nnrt = G.nrt;\nNr = G.Nr;\nnth = G.nth;\n\nh = P.thickness_m;\nrhoh = P.density * h;\nEh = P.youngs_modulus * h;\nRs = P.radius_m^2 / (2 * P.rise_m);\nht = 1 / (P.sample_rate * P.oversample);\ng = P.g_param;\nphiw = P.phi_param;\nkap = 0; % set below from cymbal2_perm (needs Dflex)\ns1 = P.sigma1 * P.radius_m^2;\ns0 = P.sigma0;\n\nMW = P.washer_mass;\nRW = P.washer_damping;\nKW = P.washer_stiffness;\n% Spring term split as phiw/2 on the (k+1, k-1) average and (1 - phiw) on\n% step k. Note the paper prints the complementary split, which is unstable\n% for a massless washer (the roots of (1-phi)/2 (z^2+1) + phi z include\n% |z| > 1 whenever phi > 1/2); with the implicit weight on the average,\n% phi = 0.7 is unconditionally stable.\nAW = MW + (ht / 2) * RW + (phiw / 2) * ht^2 * KW;\nBW = 2 * MW - (1 - phiw) * ht^2 * KW;\nCW = -MW + (ht / 2) * RW - (phiw / 2) * ht^2 * KW;\n\nAiB = zeros(nrt, nrt, nth);\nAiC = zeros(nrt, nrt, nth);\nAi = zeros(nrt, nrt, nth);\nFi = zeros(nrt, nrt, nth);\nnDropped = 0;\n\ne1 = zeros(1, nrt); e1(1) = 1;\ne3 = zeros(1, nrt); e3(3) = 1;\nI = eye(nrt);\n\nfor m = 0:nth / 2\n lam2 = (2 * cos(m * G.hth) - 2) / G.hth^2;\n op = cymbal2_perm(P, G, lam2);\n kap = op.Dflex / rhoh;\n\n Am = zeros(nrt, nrt);\n Bm = zeros(nrt, nrt);\n Cm = zeros(nrt, nrt);\n if P.support_clamped ~= 0\n l0 = 4;\n else\n l0 = 3;\n end\n for l = l0:Nr\n Am(l, :) = (1 + s0 * ht) * I(l, :) + g * kap * ht^2 * op.Bih(l, :) - s1 * ht * op.Lap(l, :);\n Bm(l, :) = 2 * I(l, :) - (1 - 2 * g) * kap * ht^2 * op.Bih(l, :);\n Cm(l, :) = (s0 * ht - 1) * I(l, :) - g * kap * ht^2 * op.Bih(l, :) - s1 * ht * op.Lap(l, :);\n end\n if P.support_clamped ~= 0\n Am(1, :) = e1;\n Am(2, :) = e3;\n Am(3, :) = G.D1(3, :);\n else\n % The shear side of the washer balance gets the same implicit split as\n % the spring; the paper evaluates it at step k only, which feeds a\n % spurious mode on the virtual lines when the washer is massless.\n Am(1, :) = op.MomW;\n Am(2, :) = AW * e3 - (phiw / 2) * ht^2 * op.ShW;\n Bm(2, :) = BW * e3 + (1 - phiw) * ht^2 * op.ShW;\n Cm(2, :) = CW * e3 + (phiw / 2) * ht^2 * op.ShW;\n end\n Am(Nr + 1, :) = op.MomE;\n Am(Nr + 2, :) = op.ShE;\n\n % Row equilibration, absorbed into the stored inverses\n sc = 1.0 ./ max(abs(Am), [], 2);\n Asc = sc .* Am;\n Aisc = inv(Asc);\n AiBm = Aisc * (sc .* Bm);\n AiCm = Aisc * (sc .* Cm);\n Aim = Aisc .* sc'; % inv(Am) applied to raw rhs\n\n % On grids too coarse to resolve a wavenumber's inner boundary layer,\n % the discretized boundary conditions can hand that wavenumber a\n % negative stiffness; such columns are truncated rather than integrated.\n % At the default grid none are.\n sr = max(abs(eig([AiBm, AiCm; I, zeros(nrt)])));\n if sr > 1 + 1e-7\n AiBm = zeros(nrt, nrt);\n AiCm = zeros(nrt, nrt);\n Aim = zeros(nrt, nrt);\n nDropped = nDropped + 1 + (m > 0 && m < nth / 2);\n end\n\n % Airy solve: biharmonic rows inside, phi = phi_r = 0 at both circles\n Fm = zeros(nrt, nrt);\n Fm(1, 1) = 1;\n Fm(2, :) = G.D1(3, :);\n Fm(3, 3) = 1;\n for l = 4:Nr - 1\n Fm(l, :) = op.Bih(l, :);\n end\n Fm(Nr, Nr) = 1;\n Fm(Nr + 1, :) = G.D1(Nr, :);\n Fm(Nr + 2, Nr + 2) = 1;\n scf = 1.0 ./ max(abs(Fm), [], 2);\n Fim = inv(scf .* Fm) .* scf';\n\n if m == 0\n cols = 1;\n elseif m == nth / 2\n cols = nth;\n else\n cols = [2 * m, 2 * m + 1];\n end\n for cc = 1:length(cols)\n AiB(:, :, cols(cc)) = AiBm;\n AiC(:, :, cols(cc)) = AiCm;\n Ai(:, :, cols(cc)) = Aim;\n Fi(:, :, cols(cc)) = Fim;\n end\nend\n\nif nDropped > nth / 2\n error(['cymbal2: the radial grid cannot resolve the inner boundary layer ' ...\n '(need roughly washer_radius > 3 (a - c) / (n_radial - 3)); raise n_radial']);\nend\n\n% Strike and pickup nodes (radial fraction 0 at the washer, 1 at the rim)\nls = 3 + round(min(1, max(0, P.strike_r)) * (Nr - 3));\nlp = 3 + round(min(1, max(0, P.pickup_r)) * (Nr - 3));\njs = 1 + mod(round(P.strike_theta_deg / 360 * nth), nth);\njp = 1 + mod(round(P.pickup_theta_deg / 360 * nth), nth);\n\nS.G = G;\nS.AiB = AiB; S.AiC = AiC; S.Ai = Ai; S.Fi = Fi;\nS.rhoh = rhoh; S.Eh = Eh; S.Rs = Rs; S.ht = ht;\nS.ls = ls; S.lp = lp; S.js = js; S.jp = jp;\nS.eomRow0 = 3 + (P.support_clamped ~= 0); % first equation-of-motion row\nS.nDropped = nDropped; % truncated theta columns\nS.strikeRadiusM = G.r(ls);\nS.pickupRadiusM = G.r(lp);\nS.forceScale = 1 / (G.r(ls) * G.hr * G.hth); % point force -> pressure\nend\n"
68+ },
69+ {
70+ "name": "cymbal2_validate.m",
71+ "content": "% Validate the cymbal2 model. There is no external reference\n% implementation, so the checks are (a) closed-form limits the\n% discretization must reproduce (the free circular plate, rigid-body\n% motion on the washer, cantilever statics of the stick), (b) agreement\n% with the linear-FEM eigenfrequencies published in Samejima (2021) for\n% the modes that model resolves the same way, (c) structural identities\n% of the transform and of the assembled update, (d) exact reproducibility,\n% and (e) pinned end-to-end diagnostics from a reduced configuration,\n% recorded when the model was built and treated as a regression reference.\n% Run with: numbl run cymbal2_validate.m\n\nfails = 0;\nP = cymbal2_defaults();\n\n% ---- 1. flat limit: the completely free circular plate ----\n% With the curvature removed and the washer stiffness zeroed, the shell\n% must reproduce the free plate's frequency equation (the 2 x 2 J_n / I_n\n% determinant under zero radial moment and zero Kelvin-Kirchhoff shear).\n% m = 0 is excluded: its limit as the free hole closes is logarithmic, so\n% a finite hole shifts it far more than the m >= 1 modes.\nQ = P;\nQ.rise_m = 1e-9;\nQ.washer_stiffness = 0;\nQ.washer_radius_m = 0.010;\nQ.n_radial = 96;\nDp = Q.youngs_modulus * Q.thickness_m^3 / (12 * (1 - Q.poisson_ratio^2));\nscal = sqrt(Dp / (Q.density * Q.thickness_m)) / Q.radius_m^2 / (2 * pi);\nerrP = 0;\nfor m = 1:4\n f = cymbal2_modes(Q, m);\n fe = f(f > 0.5);\n lam = c2v_plate_roots(m, Q.poisson_ratio);\n fc = lam.^2 * scal;\n nc = min(3, min(length(fe), length(fc)));\n errP = max(errP, max(abs(fe(1:nc) - fc(1:nc)') ./ fc(1:nc)'));\nend\nok = errP < 0.04;\nfails = fails + (~ok);\nfprintf('flat plate : max freq err %.1e over m = 1..4 %s\\n', errP, c2v_verdict(ok));\n\n% ---- 2. rigid limit: bounce and tilt on the washer ----\n% With the plate stiffened a hundredfold the lowest m = 0 and m = 1 modes\n% must approach the closed-form rigid-body values on the line spring:\n% bounce f^2 = KW 2 pi c / (rho h pi (a^2 - c^2)) / (2 pi)^2 and tilt\n% f^2 = 4 KW c^3 / (rho h (a^4 - c^4)) / (2 pi)^2.\nQ = P;\nQ.rise_m = 1e-9;\nQ.youngs_modulus = P.youngs_modulus * 100;\nQ.n_radial = 96;\na = Q.radius_m; c = Q.washer_radius_m;\nrhoh = Q.density * Q.thickness_m;\nKW = Q.washer_stiffness;\nf0r = sqrt(KW * 2 * pi * c / (rhoh * pi * (a^2 - c^2))) / (2 * pi);\nf1r = sqrt(4 * KW * c^3 / (rhoh * (a^4 - c^4))) / (2 * pi);\nfa = cymbal2_modes(Q, 0); fa = fa(fa > 0.01);\nfb = cymbal2_modes(Q, 1); fb = fb(fb > 0.01);\ne0 = abs(fa(1) / f0r - 1);\ne1 = abs(fb(1) / f1r - 1);\nok = e0 < 0.03 && e1 < 0.03;\nfails = fails + (~ok);\nfprintf('washer : bounce err %.1e, tilt err %.1e vs rigid closed form %s\\n', e0, e1, c2v_verdict(ok));\n\n% ---- 3. the paper's linear FEM reference, where comparable ----\n% Samejima's Fig. 9 gives linear-FEM eigenfrequencies of the same shell.\n% For the clamped center, the (2,1), (3,1), (4,1) modes must agree to a\n% few percent ((1,1) is excluded: it is controlled by the support model,\n% where this shell's moment-free hole and their 2-D elements differ).\nQ = P;\nQ.support_clamped = 1;\nref = [20.1, 43.3, 77.7];\ngot = zeros(1, 3);\nfor m = 2:4\n f = cymbal2_modes(Q, m);\n f = f(f > 1);\n got(m - 1) = f(1);\nend\nerrF = max(abs(got - ref) ./ ref);\nok = errF < 0.06;\nfails = fails + (~ok);\nfprintf('paper FEM : clamped (2,1)(3,1)(4,1) = %.1f %.1f %.1f Hz (ref %.1f %.1f %.1f), err %.1e %s\\n', ...\n got(1), got(2), got(3), ref(1), ref(2), ref(3), errF, c2v_verdict(ok));\n\n% ---- 4. transform and update-structure identities ----\n% The azimuthal transform must be exactly invertible, and the physical\n% derivative pipeline must agree with the per-wavenumber operators: for a\n% field cos(m theta) g(r), the physical-space Laplacian must equal the\n% per-m radial Laplacian applied to g.\nG = cymbal2_grid(P);\neT = max(max(abs(G.F * G.Gi - eye(G.nth))));\ngr = exp(-((G.r - 0.12) / 0.05).^2);\nm = 7;\nUf = gr * cos(m * G.th');\nLapPhys = G.D2 * Uf + G.Rd .* (G.D1 * Uf) + G.R2d .* (Uf * G.T2);\nlam2 = (2 * cos(m * G.hth) - 2) / G.hth^2;\nop = cymbal2_perm(P, G, lam2);\nLapM = (op.Lap * gr) * cos(m * G.th');\nr2 = 2:G.nrt - 1;\neL = max(max(abs(LapPhys(r2, :) - LapM(r2, :)))) / max(max(abs(LapM(r2, :))));\nok = eT < 1e-10 && eL < 1e-10;\nfails = fails + (~ok);\nfprintf('operators : transform inverse %.1e, physical vs per-m Laplacian %.1e %s\\n', ...\n eT, eL, c2v_verdict(ok));\n\n% ---- 5. the assembled update is not amplifying ----\n% Every azimuthal column's two-step update must have spectral radius at\n% most 1 (the damping makes physical modes strictly contractive; boundary\n% constraints contribute zeros). The system builder truncates columns\n% that fail; at the defaults none may be.\nS = cymbal2_system(P);\nnrt = S.G.nrt;\nsrMax = 0;\nfor col = [1, 2, S.G.nth / 2, S.G.nth]\n T = [S.AiB(:, :, col), S.AiC(:, :, col); eye(nrt), zeros(nrt)];\n srMax = max(srMax, max(abs(eig(T))));\nend\nok = S.nDropped == 0 && srMax <= 1 + 1e-9;\nfails = fails + (~ok);\nfprintf('stability : %d truncated columns, sampled spectral radius %.9f %s\\n', ...\n S.nDropped, srMax, c2v_verdict(ok));\n\n% ---- 6. stick statics ----\n% A static tip load on the clamped-hand stick must deflect by\n% F L^3 / (3 E I); the rotational-spring hand adds F L^2 / K_R times L.\nMBc = cymbal2_mallet(setfield(P, 'grip_clamped', 1));\nMBs = cymbal2_mallet(P);\nEI = P.stick_youngs * P.stick_inertia_m4;\nL = P.stick_length_m;\nFl = 1.0;\n% static balance of the update: at a fixed point, (I - A^-1 B - A^-1 C) d\n% equals A^-1 ht^2 K d, so the deflection under a constant tip load F is\n% (I - ASiB - ASiC) \\ (tipcol * ht^2 * F)\nhtv = MBc.ht;\ndc = ((eye(size(MBc.ASiB, 1)) - MBc.ASiB - MBc.ASiC) \\ (MBc.tipcol * htv^2)) * Fl;\nds = ((eye(size(MBs.ASiB, 1)) - MBs.ASiB - MBs.ASiC) \\ (MBs.tipcol * htv^2)) * Fl;\nwc = dc(MBc.tipi);\nws = ds(MBs.tipi);\nwcRef = Fl * L^3 / (3 * EI);\nwsRef = wcRef + (Fl * L / P.grip_spring) * L; % grip rotation F L / K_R times the arm\nec = abs(wc / wcRef - 1);\nes = abs(ws / wsRef - 1);\nok = ec < 1e-6 && es < 1e-6;\nfails = fails + (~ok);\nfprintf('stick : clamped-hand tip %.3e (ref %.3e), spring-hand %.3e (ref %.3e) %s\\n', ...\n wc, wcRef, ws, wsRef, c2v_verdict(ok));\n\n% ---- 7. reproducibility and pinned diagnostics ----\n% Two identical renders of a reduced configuration must agree exactly,\n% and its diagnostics are pinned. The trajectory of a nonlinear shell is\n% sensitive, so the acoustic pins carry percent-level tolerances while\n% the contact pins are tight.\nQ = P;\nQ.sample_rate = 44100;\nQ.oversample = 1;\nQ.tail_s = 1.5;\n[y1, d1] = cymbal2_render(Q);\n[y2, d2] = cymbal2_render(Q);\nrms1 = sqrt(mean(y1.^2));\nrms2 = sqrt(mean(y2.^2));\nsame = max(abs(y1 - y2)) == 0;\nok = same ...\n && d1.contactEvents == d2.contactEvents ...\n && abs(d1.peakForceN - 311.14) < 4 ...\n && abs(d1.contactVel(1) + 8.223) < 0.08 ...\n && abs(rms1 - 0.0973) < 0.010;\nfails = fails + (~ok);\nfprintf('pinned : repro %d, peakF %.2f N (ref 311.14), cvel %.3f (ref -8.223), rms %.4f (ref 0.0973) %s\\n', ...\n same, d1.peakForceN, d1.contactVel(1), rms1, c2v_verdict(ok));\n\n% ---- 8. the crash decays ----\nn = length(y1);\nfs = Q.sample_rate;\nhead = y1(1:round(0.2 * fs));\ntail = y1(n - round(0.2 * fs):n);\ndropE = 10 * log10(sum(tail.^2) / sum(head.^2));\nok = dropE < -3 && dropE > -50;\nfails = fails + (~ok);\nfprintf('decay : last/first 0.2 s %.1f dB %s\\n', dropE, c2v_verdict(ok));\n\nif fails > 0\n fprintf('%d CHECK(S) FAILED\\n', fails);\nelse\n fprintf('all checks passed\\n');\nend\n\nfunction lam = c2v_plate_roots(n, nu)\n% Roots of the free-edge frequency equation of a circular plate\nlam = [];\nprev = c2v_plate_det(n, nu, 0.30);\nfor x = 0.30:0.01:20\n cur = c2v_plate_det(n, nu, x + 0.01);\n if prev * cur < 0\n lo = x; hi = x + 0.01;\n for it = 1:60\n mid = (lo + hi) / 2;\n if c2v_plate_det(n, nu, lo) * c2v_plate_det(n, nu, mid) < 0\n hi = mid;\n else\n lo = mid;\n end\n end\n lam = [lam, (lo + hi) / 2];\n end\n prev = cur;\nend\nend\n\nfunction d = c2v_plate_det(n, nu, x)\nJ = besselj(n, x); I = besseli(n, x);\nJp = besselj(n - 1, x) - (n / x) * J;\nIp = besseli(n - 1, x) - (n / x) * I;\nJpp = -Jp / x + (n^2 / x^2 - 1) * J;\nIpp = -Ip / x + (n^2 / x^2 + 1) * I;\na11 = Jpp + nu * (Jp / x - n^2 * J / x^2);\na12 = Ipp + nu * (Ip / x - n^2 * I / x^2);\na21 = -Jp - (1 - nu) * (n^2 / x^2) * (Jp - J / x);\na22 = Ip - (1 - nu) * (n^2 / x^2) * (Ip - I / x);\nd = a11 * a22 - a12 * a21;\nend\n\nfunction s = c2v_verdict(ok)\nif ok\n s = 'ok';\nelse\n s = 'MISMATCH';\nend\nend\n"
72+ }
73+ ]
74+}
matlab/cymbal2_defaults.madded+89−0View file
@@ -0,0 +1,89 @@
1+function P = cymbal2_defaults()
2+% CYMBAL2_DEFAULTS Default parameters, following Samejima (2021) Table 1:
3+% an 18-inch medium crash (45 cm diameter) modeled as a shallow spherical
4+% shell, supported at the center by a washer of finite impedance and struck
5+% by a mallet whose stick is an elastic beam held in a compliant grip.
6+% The output duration is not a parameter: it is
7+% (strike_count - 1) / strike_rate_hz + tail_s.
8+
9+% Shell geometry (a shallow spherical cap; rise sets the radius of curvature
10+% Rs = radius^2 / (2 * rise))
11+P.radius_m = 0.230;
12+P.rise_m = 0.040;
13+P.thickness_m = 0.001;
14+
15+% Material (bronze, per the paper's measured cymbal)
16+P.youngs_modulus = 1.1e11;
17+P.poisson_ratio = 0.33;
18+P.density = 8860.0;
19+
20+% Damping: decay rate sigma0 + sigma1 * a^2 * k^2 at flexural wavenumber k
21+P.sigma0 = 1.34; % s^-1, frequency-independent
22+P.sigma1 = 1.2e-3; % s^-1 per unit nondimensional k^2
23+
24+% Support at the center: an SDOF washer per unit length of the inner
25+% circumference, entering as the shear boundary condition at r = washer
26+% radius. support_clamped = 1 replaces it with the rigid clamp of the
27+% earlier literature (Bilbao 2010), which is worth hearing once.
28+P.washer_radius_m = 0.029;
29+P.washer_stiffness = 1.65e5; % N/m per m
30+P.washer_mass = 0.0; % kg per m (paper sets 0)
31+P.washer_damping = 0.0; % N s/m per m (paper sets 0)
32+P.support_clamped = 0;
33+
34+% Mallet head (SDOF, Hertz-type nonlinear stiffness f = K * delta^alpha)
35+P.head_radius_m = 0.025;
36+P.head_mass_kg = 0.028;
37+P.head_stiffness = 1.6e8; % N m^-alpha
38+P.head_alpha = 2.54;
39+
40+% Mallet stick (Euler-Bernoulli beam, FEM)
41+P.stick_length_m = 0.390;
42+P.stick_area_m2 = 1.43139e-4;
43+P.stick_youngs = 8.8e9;
44+P.stick_density = 400.0;
45+P.stick_inertia_m4 = 1.63044e-9;
46+P.stick_damping_rate = 10.0; % mass-proportional, R = rate * M
47+P.n_beam = 20; % beam elements (21 nodes)
48+
49+% Grip: 0 = rotational spring hand (loose grip), 1 = clamped hand (tight)
50+P.grip_clamped = 0;
51+P.grip_spring = 140.0; % N m/rad
52+
53+% Strike gesture: the hand accelerates downward to strike_velocity, coasts
54+% until the head touches, holds briefly, and retracts to its rest height.
55+% The paper prescribes constant acceleration to the same terminal speed
56+% but never reverses, so it cannot play a second stroke. Because the stick
57+% is compliant, the tip arrives somewhat slower than the hand.
58+P.strike_velocity = 5.75; % peak hand speed, m/s; reproduces the
59+ % character of the authors' published demo
60+P.strike_accel = 250.0; % downstroke hand acceleration, m/s^2
61+P.retract_accel = 2500.0; % wrist-snap retract acceleration, m/s^2
62+P.strike_hold_ms = 0.5; % dwell between contact and retract
63+P.strike_gap_m = -1; % rest clearance; <= 0 sizes it to the stroke
64+P.strike_r = 0.6; % contact radius, 0 at the washer, 1 at the rim
65+P.strike_theta_deg = 180.0;
66+P.strike_count = 1;
67+P.strike_rate_hz = 4.0;
68+P.velocity_spread = 0.0; % relative velocity scatter of repeated strokes
69+P.seed = 1;
70+
71+% Pickup: normal acceleration of the shell at one point, as in the paper
72+P.pickup_r = 0.6;
73+P.pickup_theta_deg = 0.0;
74+
75+% Discretization
76+P.n_radial = 24; % radial grid intervals (the paper's N_r)
77+P.n_theta = 60; % azimuthal grid points (the paper's N_theta)
78+P.g_param = 0.2498; % free parameter of the implicit shell scheme
79+P.phi_param = 0.7; % free parameter of the washer and beam schemes
80+
81+% Rendering / output
82+P.sample_rate = 48000;
83+P.oversample = 2; % simulation runs at oversample * sample_rate
84+P.tail_s = 6.0; % the low modes carry T60 of about 5 s
85+P.pre_delay_ms = 3.0;
86+P.normalize = 1;
87+P.peak_dbfs = -3.0;
88+P.full_scale_accel = 4.0e4; % m/s^2 mapped to full scale when normalize = 0
89+end
matlab/cymbal2_demo.madded+32−0View file
@@ -0,0 +1,32 @@
1+% Render a cymbal crash with default parameters and print diagnostics.
2+% Set P.strike_velocity = 0.5 for a tap on the same instrument, or
3+% P.support_clamped = 1 to hear the rigid center clamp of the earlier
4+% literature instead of the washer. P.strike_count = 4 plays a repeated
5+% figure; the hand retracts between strokes, which the paper's
6+% constant-acceleration hand cannot do.
7+% Run with: numbl run cymbal2_demo.m
8+
9+P = cymbal2_defaults();
10+tic;
11+[y, dinfo] = cymbal2_render(P);
12+elapsed = toc;
13+
14+fprintf('grid : %d radial x %d azimuthal, %.1f kHz simulation rate\n', ...
15+ P.n_radial, P.n_theta, dinfo.simRate / 1000);
16+fprintf('shell : a = %.0f cm, rise %.0f mm -> Rs = %.2f m, washer at %.1f cm\n', ...
17+ P.radius_m * 100, P.rise_m * 1000, dinfo.RsM, P.washer_radius_m * 100);
18+f0 = cymbal2_modes(P, 0); f0 = f0(f0 > 1);
19+f1 = cymbal2_modes(P, 1); f1 = f1(f1 > 1);
20+f2 = cymbal2_modes(P, 2); f2 = f2(f2 > 1);
21+fprintf('low modes : %.1f Hz (m=1 tilt), %.1f Hz (m=2), %.1f Hz (m=0 bounce)\n', ...
22+ f1(1), f2(1), f0(1));
23+fprintf('stroke : hand %.1f m/s -> tip %.2f m/s at contact\n', ...
24+ dinfo.meanV, abs(dinfo.contactVel(1)));
25+fprintf('contact : %.0f N peak, %.2f mm compression, %.1f ms over %d event(s)\n', ...
26+ dinfo.peakForceN, dinfo.peakXiMm, dinfo.contactMs, dinfo.contactEvents);
27+fprintf('shell peak : %.2f mm displacement = %.1f thicknesses\n', ...
28+ dinfo.peakDispM * 1000, dinfo.peakDispThickness);
29+fprintf('pickup : %.2e m/s^2 peak acceleration at r = %.1f cm\n', ...
30+ dinfo.peakAccel, dinfo.pickupRadiusM * 100);
31+fprintf('length : %.2f s @ %d Hz\n', dinfo.lengthS, P.sample_rate);
32+fprintf('render time: %.1f s\n', elapsed);
matlab/cymbal2_grid.madded+88−0View file
@@ -0,0 +1,88 @@
1+function G = cymbal2_grid(P)
2+% CYMBAL2_GRID Polar grid, radial difference matrices, and the azimuthal
3+% transform. The radial grid carries l = 1..Nr+2 lines: real nodes at
4+% l = 3..Nr spanning [washer radius, rim radius], and two virtual lines
5+% beyond each boundary, eliminated by the boundary-condition rows of the
6+% assembled systems (Samejima 2021, Fig. 4). Azimuthal operators are
7+% circulant, so every linear solve block-diagonalizes over the real
8+% cos/sin basis; F maps physical theta samples to coefficients (U * F)
9+% and Gi maps back (Uhat * Gi), with mcol giving each column's wavenumber.
10+
11+a = P.radius_m;
12+c = P.washer_radius_m;
13+Nr = P.n_radial;
14+nth = P.n_theta;
15+if mod(nth, 2) ~= 0
16+ error('cymbal2: n_theta must be even');
17+end
18+if Nr < 8
19+ error('cymbal2: n_radial must be at least 8');
20+end
21+nrt = Nr + 2;
22+hr = (a - c) / (Nr - 3);
23+if c <= hr
24+ error('cymbal2: washer_radius too small for n_radial (need c > (a-c)/(n_radial-3))');
25+end
26+
27+r = c + ((1:nrt)' - 3) * hr; % r(1) = c - 2hr may be tiny on coarse grids;
28+r(1) = max(r(1), hr / 2); % row 1 of every operator is never used
29+
30+% Radial stencils (dense; nrt is small). Rows outside the valid range are
31+% left zero and are overwritten by boundary-condition rows downstream.
32+D1 = zeros(nrt, nrt);
33+D2 = zeros(nrt, nrt);
34+D3 = zeros(nrt, nrt);
35+for l = 2:nrt - 1
36+ D1(l, l - 1) = -1 / (2 * hr); D1(l, l + 1) = 1 / (2 * hr);
37+ D2(l, l - 1) = 1 / hr^2; D2(l, l) = -2 / hr^2; D2(l, l + 1) = 1 / hr^2;
38+end
39+for l = 3:nrt - 2
40+ D3(l, l - 2) = -1 / (2 * hr^3);
41+ D3(l, l - 1) = 2 / (2 * hr^3);
42+ D3(l, l + 1) = -2 / (2 * hr^3);
43+ D3(l, l + 2) = 1 / (2 * hr^3);
44+end
45+
46+% Azimuthal stencils as dense circulant matrices, applied on the right
47+% (U * T2 etc., so they carry a transpose; T2 is symmetric, T1 is
48+% antisymmetric and T1t below is arranged so Uth = U * T1t).
49+hth = 2 * pi / nth;
50+T2 = zeros(nth, nth);
51+T1t = zeros(nth, nth);
52+for j = 1:nth
53+ jm = 1 + mod(j - 2, nth);
54+ jp = 1 + mod(j, nth);
55+ T2(j, j) = -2 / hth^2; T2(j, jm) = 1 / hth^2; T2(j, jp) = 1 / hth^2;
56+ % row j of T1 has -1 at jm, +1 at jp (scaled); transpose: column j
57+ T1t(jm, j) = -1 / (2 * hth);
58+ T1t(jp, j) = 1 / (2 * hth);
59+end
60+
61+% Real transform: columns ordered m = 0, (1c,1s), (2c,2s), ..., nth/2
62+th = ((1:nth)' - 1) * hth;
63+F = zeros(nth, nth);
64+Gi = zeros(nth, nth);
65+mcol = zeros(1, nth);
66+F(:, 1) = 1 / nth;
67+Gi(1, :) = 1;
68+col = 2;
69+for m = 1:nth / 2 - 1
70+ F(:, col) = 2 * cos(m * th) / nth; Gi(col, :) = cos(m * th)'; mcol(col) = m; col = col + 1;
71+ F(:, col) = 2 * sin(m * th) / nth; Gi(col, :) = sin(m * th)'; mcol(col) = m; col = col + 1;
72+end
73+F(:, col) = cos((nth / 2) * th) / nth;
74+Gi(col, :) = cos((nth / 2) * th)';
75+mcol(col) = nth / 2;
76+
77+% Eigenvalue of the circulant second-difference for wavenumber m
78+lam2 = (2 * cos(mcol * hth) - 2) / hth^2;
79+
80+G.a = a; G.c = c; G.Nr = Nr; G.nrt = nrt; G.nth = nth;
81+G.hr = hr; G.hth = hth;
82+G.r = r; G.th = th;
83+G.D1 = D1; G.D2 = D2; G.D3 = D3;
84+G.T1t = T1t; G.T2 = T2;
85+G.F = F; G.Gi = Gi; G.mcol = mcol; G.lam2 = lam2;
86+G.Rd = 1.0 ./ r; % 1/r column, broadcast over theta
87+G.R2d = 1.0 ./ r.^2;
88+end
matlab/cymbal2_mallet.madded+82−0View file
@@ -0,0 +1,82 @@
1+function MB = cymbal2_mallet(P)
2+% CYMBAL2_MALLET Precompute the striking mechanism: the mallet stick as an
3+% Euler-Bernoulli beam (Hermite cubic elements, consistent mass,
4+% mass-proportional damping), held at one end by the hand, with the mallet
5+% head attached at the other end as a separate SDOF body (Samejima 2021,
6+% Sec. 3). The hand grips either rigidly (grip_clamped = 1: transverse
7+% displacement and rotation both prescribed) or through a rotational
8+% spring (grip_clamped = 0: displacement prescribed, rotation free against
9+% grip_spring). The head couples to the beam tip by the rigid connection
10+% u_M = u_S. The paper integrates the head as its own SDOF (Eq. 34) and
11+% solves the connection force from the constraint; because that couples an
12+% explicit head update to the implicit beam, the discrete pair is
13+% unstable, so the head mass is lumped onto the tip node instead, which is
14+% the same physics under the rigid connection (substituting u_M = u_S and
15+% f_S = -f_M2 into Eq. (31) puts M_M on the tip and f_M1 on the beam).
16+
17+ne = P.n_beam;
18+nn = ne + 1;
19+ndof = 2 * nn;
20+le = P.stick_length_m / ne;
21+EI = P.stick_youngs * P.stick_inertia_m4;
22+rA = P.stick_density * P.stick_area_m2;
23+
24+Ke = EI / le^3 * [12, 6 * le, -12, 6 * le; ...
25+ 6 * le, 4 * le^2, -6 * le, 2 * le^2; ...
26+ -12, -6 * le, 12, -6 * le; ...
27+ 6 * le, 2 * le^2, -6 * le, 4 * le^2];
28+Me = rA * le / 420 * [156, 22 * le, 54, -13 * le; ...
29+ 22 * le, 4 * le^2, 13 * le, -3 * le^2; ...
30+ 54, 13 * le, 156, -22 * le; ...
31+ -13 * le, -3 * le^2, -22 * le, 4 * le^2];
32+
33+K = zeros(ndof, ndof);
34+M = zeros(ndof, ndof);
35+for e = 1:ne
36+ ix = (2 * e - 1):(2 * e + 2);
37+ K(ix, ix) = K(ix, ix) + Ke;
38+ M(ix, ix) = M(ix, ix) + Me;
39+end
40+R = P.stick_damping_rate * M;
41+M(ndof - 1, ndof - 1) = M(ndof - 1, ndof - 1) + P.head_mass_kg;
42+
43+% Node 1 is the grip, node nn the tip. Dof 1 (grip displacement) is always
44+% prescribed to the hand motion; dof 2 (grip rotation) is prescribed to 0
45+% for the clamped hand, or free with the rotational spring added.
46+if P.grip_clamped ~= 0
47+ fi = 3:ndof;
48+else
49+ K(2, 2) = K(2, 2) + P.grip_spring;
50+ fi = 2:ndof;
51+end
52+Kf = K(fi, fi); Mf = M(fi, fi); Rf = R(fi, fi);
53+Kc = K(fi, 1); Mc = M(fi, 1); Rc = R(fi, 1);
54+
55+ht = 1 / (P.sample_rate * P.oversample);
56+ps = P.phi_param;
57+% Stiffness split with the implicit weight ps/2 on the (k+1, k-1) average,
58+% as in the shell's washer row; the complementary split printed in the
59+% paper is unstable for the beam's highest elements at audio rates.
60+AS = Mf + (ht / 2) * Rf + (ps / 2) * ht^2 * Kf;
61+BS = 2 * Mf - (1 - ps) * ht^2 * Kf;
62+CS = -Mf + (ht / 2) * Rf - (ps / 2) * ht^2 * Kf;
63+ASi = inv(AS);
64+MB.ASiB = ASi * BS;
65+MB.ASiC = ASi * CS;
66+
67+% Hand coupling, folded so the per-step load is
68+% vM * hA(k) + vR * hB(k) + vK * hC(k)
69+% with hA = -(uH+ - 2 uH + uH-), hB = -(ht/2)(uH+ - uH-),
70+% hC = -ht^2 (ps uH + ((1-ps)/2)(uH+ + uH-)).
71+MB.vM = ASi * Mc;
72+MB.vR = ASi * Rc;
73+MB.vK = ASi * Kc;
74+
75+tipi = find(fi == ndof - 1);
76+MB.tipi = tipi;
77+MB.tipcol = ASi(:, tipi);
78+MB.gtip = ASi(tipi, tipi);
79+MB.nfree = length(fi);
80+MB.ht = ht;
81+MB.ps = ps;
82+end
matlab/cymbal2_modes.madded+71−0View file
@@ -0,0 +1,71 @@
1+function [f, W] = cymbal2_modes(P, m)
2+% CYMBAL2_MODES Eigenfrequencies (Hz) of the linearized shell for one
3+% azimuthal wavenumber m, with the membrane (curvature) stiffening and the
4+% support boundary condition included; the washer is taken as its static
5+% stiffness (its mass and damping are ignored here, as in the paper's
6+% linear FEM reference). Used by the demo and the validation, not by the
7+% renderer. W returns the mode shapes on the real radial nodes l = 3..Nr.
8+
9+G = cymbal2_grid(P);
10+nrt = G.nrt;
11+Nr = G.Nr;
12+hth = G.hth;
13+lam2 = (2 * cos(m * hth) - 2) / hth^2;
14+op = cymbal2_perm(P, G, lam2);
15+
16+E = P.youngs_modulus;
17+h = P.thickness_m;
18+rhoh = P.density * h;
19+Rs = P.radius_m^2 / (2 * P.rise_m);
20+
21+% Airy solve operator: Fm * phi = Z * (E h / Rs) * Lap * u, with
22+% phi = phi_r = 0 at both circles and identity rows on the outermost
23+% virtual lines.
24+Fm = zeros(nrt, nrt);
25+Fm(1, 1) = 1;
26+Fm(2, :) = G.D1(3, :);
27+Fm(3, 3) = 1;
28+for l = 4:Nr - 1
29+ Fm(l, :) = op.Bih(l, :);
30+end
31+Fm(Nr, Nr) = 1;
32+Fm(Nr + 1, :) = G.D1(Nr, :);
33+Fm(Nr + 2, Nr + 2) = 1;
34+Zl = zeros(nrt, nrt);
35+for l = 4:Nr - 1
36+ Zl(l, :) = op.Lap(l, :);
37+end
38+Phimap = Fm \ ((E * h / Rs) * Zl);
39+
40+% Full stiffness on the equation-of-motion rows l = 3..Nr
41+Kfull = zeros(nrt, nrt);
42+Kterm = op.Dflex * op.Bih + (1 / Rs) * op.Lap * Phimap;
43+for l = 3:Nr
44+ Kfull(l, :) = Kterm(l, :);
45+end
46+
47+% Boundary rows, static. With the washer the shell's equation of motion
48+% holds on l = 3..Nr and the four virtual lines are eliminated; with the
49+% rigid clamp u and u_r vanish at l = 3, so u(3) itself joins the
50+% eliminated set and the equations of motion start at l = 4 (the innermost
51+% virtual line is then never engaged).
52+e3 = zeros(1, nrt); e3(3) = 1;
53+if P.support_clamped ~= 0
54+ Cmat = [e3; G.D1(3, :); op.MomE; op.ShE];
55+ keep = 4:Nr;
56+ elim = [2, 3, Nr + 1, Nr + 2];
57+else
58+ Cmat = [op.MomW; op.ShW - P.washer_stiffness * e3; op.MomE; op.ShE];
59+ keep = 3:Nr;
60+ elim = [1, 2, Nr + 1, Nr + 2];
61+end
62+X = -Cmat(:, elim) \ Cmat(:, keep);
63+Kred = Kfull(keep, keep) + Kfull(keep, elim) * X;
64+
65+[V, Dm] = eig(Kred / rhoh);
66+om2 = real(diag(Dm));
67+[om2, idx] = sort(om2);
68+om2 = max(om2, 0);
69+f = sqrt(om2) / (2 * pi);
70+W = real(V(:, idx));
71+end
matlab/cymbal2_output.madded+67−0View file
@@ -0,0 +1,67 @@
1+function [y, peakAccel] = cymbal2_output(P, acc)
2+% CYMBAL2_OUTPUT Output conditioning: DC blocker on the raw pickup
3+% acceleration, Kaiser-windowed-sinc anti-alias decimation from the
4+% simulation rate to the output rate, tail fade, and normalization. Same
5+% structure as the sibling instruments' output stages.
6+
7+peakAccel = max(abs(acc));
8+R = P.oversample;
9+T = length(acc);
10+
11+% One-pole DC blocker (the pickup is an acceleration, but the contact
12+% leaves a small quasi-static offset while the washer relaxes)
13+fc = 15.0;
14+a1 = exp(-2 * pi * fc / (P.sample_rate * R));
15+p = zeros(1, T);
16+xp = 0; yp = 0;
17+for k = 1:T
18+ yv = acc(k) - xp + a1 * yp;
19+ xp = acc(k);
20+ yp = yv;
21+ p(k) = yv;
22+end
23+
24+if R > 1
25+ taps = 12 * R + 1;
26+ M2 = (taps - 1) / 2;
27+ fcNorm = 0.45 / R;
28+ m = (0:taps - 1) - M2;
29+ hh = sin(2 * pi * fcNorm * m) ./ (pi * m);
30+ hh(M2 + 1) = 2 * fcNorm;
31+ r = m / M2;
32+ beta = 8;
33+ hh = hh .* (c2_i0(beta * sqrt(max(0, 1 - r.^2))) / c2_i0(beta));
34+ hh = hh / sum(hh);
35+ pf = conv(p, hh);
36+ outLen = floor(T / R);
37+ y = reshape(pf((0:outLen - 1) * R + M2 + 1), 1, outLen);
38+else
39+ y = p;
40+end
41+
42+% Squared fade on the tail
43+fade = min(2048, floor(length(y) / 8));
44+if fade > 1
45+ y(end - fade + 1:end) = y(end - fade + 1:end) .* linspace(1, 0, fade).^2;
46+end
47+
48+if P.normalize ~= 0
49+ pk = max(abs(y));
50+ if pk > 0
51+ y = y / pk * 10^(P.peak_dbfs / 20);
52+ end
53+else
54+ y = y / P.full_scale_accel;
55+end
56+end
57+
58+function v = c2_i0(x)
59+% Modified Bessel function of the first kind, order zero (series)
60+v = ones(size(x));
61+term = ones(size(x));
62+hx = x / 2;
63+for k = 1:40
64+ term = term .* (hx / k).^2;
65+ v = v + term;
66+end
67+end
matlab/cymbal2_perm.madded+38−0View file
@@ -0,0 +1,38 @@
1+function op = cymbal2_perm(P, G, lam2)
2+% CYMBAL2_PERM Radial operators for one azimuthal wavenumber, dimensional.
3+% lam2 is the (negative) eigenvalue of the circulant second-difference in
4+% theta for this wavenumber, so every d^2/dtheta^2 becomes multiplication
5+% by lam2. Returns the Laplacian, the biharmonic (as the exact discrete
6+% composition Lap * Lap), and the radial-moment and Kelvin-Kirchhoff shear
7+% rows used as boundary conditions at the washer (l = 3) and at the free
8+% edge (l = Nr). The moment and shear expressions are those of Samejima
9+% (2021), Eqs. (6)-(7) and (14)-(15).
10+
11+nrt = G.nrt;
12+nu = P.poisson_ratio;
13+Dfl = P.youngs_modulus * P.thickness_m^3 / (12 * (1 - nu^2));
14+
15+Rm = diag(G.Rd);
16+R2 = diag(G.R2d);
17+R3 = diag(G.Rd.^3);
18+
19+Lap = G.D2 + Rm * G.D1 + lam2 * R2;
20+Bih = Lap * Lap;
21+
22+% Radial bending moment (per unit length, divided by -D):
23+% u_rr + nu * (u_r / r + u_thth / r^2)
24+Mom = G.D2 + nu * (Rm * G.D1 + lam2 * R2);
25+
26+% Kelvin-Kirchhoff shear: V = -D * [u_rrr + u_rr / r - u_r / r^2
27+% + (nu - 3) / r^3 * u_thth + (2 - nu) / r^2 * u_r_thth]
28+Sh = -Dfl * (G.D3 + Rm * G.D2 - R2 * G.D1 ...
29+ + (nu - 3) * lam2 * R3 + (2 - nu) * lam2 * R2 * G.D1);
30+
31+op.Lap = Lap;
32+op.Bih = Bih;
33+op.MomW = Mom(3, :); % moment row at the washer circle
34+op.MomE = Mom(nrt - 2, :); % moment row at the free edge (l = Nr)
35+op.ShW = Sh(3, :);
36+op.ShE = Sh(nrt - 2, :);
37+op.Dflex = Dfl;
38+end
matlab/cymbal2_render.madded+14−0View file
@@ -0,0 +1,14 @@
1+function [y, dinfo] = cymbal2_render(P)
2+% CYMBAL2_RENDER Full render: precompute -> time integration -> output.
3+
4+S = cymbal2_system(P);
5+MB = cymbal2_mallet(P);
6+[acc, dinfo] = cymbal2_sim(P, S, MB);
7+[y, peakAccel] = cymbal2_output(P, acc);
8+
9+dinfo.peakAccel = peakAccel;
10+dinfo.strikeRadiusM = S.strikeRadiusM;
11+dinfo.pickupRadiusM = S.pickupRadiusM;
12+dinfo.RsM = S.Rs;
13+dinfo.lengthS = length(y) / P.sample_rate;
14+end
matlab/cymbal2_sim.madded+223−0View file
@@ -0,0 +1,223 @@
1+function [acc, dinfo, fhist] = cymbal2_sim(P, S, MB)
2+% CYMBAL2_SIM Time-domain integration of the full system: von Karman
3+% shallow spherical shell (implicit in the stiff linear terms, explicit in
4+% the nonlinearity), washer support built into the boundary rows, mallet
5+% head + stick coupled through the Hertz contact. Returns the normal
6+% acceleration of the shell at the pickup node at the simulation rate,
7+% diagnostics, and the contact force history.
8+
9+G = S.G;
10+nrt = G.nrt;
11+Nr = G.Nr;
12+nth = G.nth;
13+ht = S.ht;
14+rhoh = S.rhoh;
15+Rs = S.Rs;
16+EhRs = S.Eh / Rs;
17+Eh2 = S.Eh / 2;
18+ls = S.ls; js = S.js; lp = S.lp; jp = S.jp;
19+fsc = S.forceScale;
20+nlr0 = S.eomRow0;
21+
22+K = max(1, round(P.strike_count));
23+
24+% Stroke schedule; per-stroke peak hand speeds with a seeded spread. Each
25+% stroke is launched early by its predicted hand travel time (acceleration
26+% plus coast over the rest clearance), so the contacts land near the
27+% nominal beats instead of drifting with velocity. A nonpositive
28+% strike_gap_m asks for the clearance to be sized to the stroke: contact
29+% arrives just after the hand reaches full speed, which keeps the tip's
30+% whip phase, and with it the contact speed, consistent across strokes.
31+tstr = P.pre_delay_ms * 1e-3 + (0:K - 1) / max(1e-3, P.strike_rate_hz);
32+vstr = zeros(1, K);
33+state = P.seed * 2654435761;
34+for i = 1:K
35+ state = mod(state * 1664525 + 1013904223, 4294967296);
36+ x = state / 4294967296;
37+ vstr(i) = max(0.05, P.strike_velocity * (1 + P.velocity_spread * (2 * x - 1)));
38+end
39+gap = P.strike_gap_m;
40+if gap <= 0
41+ gap = 1.05 * max(vstr)^2 / (2 * P.strike_accel) + 0.005;
42+end
43+travel = gap ./ vstr + vstr / (2 * P.strike_accel);
44+tlaunch = tstr - travel;
45+off = -min(0, min(tlaunch));
46+tlaunch = tlaunch + off;
47+
48+duration = off + (K - 1) / max(1e-3, P.strike_rate_hz) + P.tail_s;
49+nsteps = floor(duration * P.sample_rate * P.oversample);
50+
51+KM = P.head_stiffness;
52+alM = P.head_alpha;
53+tipi = MB.tipi;
54+ps = MB.ps;
55+
56+% Hand state machine: the hand accelerates downward, coasts at the stroke
57+% speed until the head makes contact (or a miss depth), holds briefly,
58+% then retracts at half speed and brakes so it arrives back at rest,
59+% ready for the next scheduled stroke. This replaces the paper's
60+% constant-acceleration hand, which never reverses. The machine runs one
61+% step ahead because the beam coupling needs u_H at k + 1.
62+aH = P.strike_accel;
63+aHr = P.retract_accel;
64+holdSteps = max(1, round(P.strike_hold_ms * 1e-3 / ht));
65+missDepth = gap + 0.06;
66+uHm = 0; uH0 = 0; uHp = 0;
67+vH = 0;
68+phase = 0;
69+strokeIdx = 1;
70+holdLeft = 0;
71+vtar = 0;
72+
73+D1p = G.D1; D2p = G.D2;
74+T1t = G.T1t; T2 = G.T2;
75+F = G.F; Gi = G.Gi;
76+Rd = G.Rd; R2d = G.R2d;
77+AiB = S.AiB; AiC = S.AiC; Ai = S.Ai; Fi = S.Fi;
78+sNL = ht^2 / rhoh;
79+
80+U = zeros(nrt, nth);
81+U1 = zeros(nrt, nth);
82+Uh3 = zeros(nrt, 1, nth);
83+Uh13 = zeros(nrt, 1, nth);
84+d = zeros(MB.nfree, 1);
85+d1 = zeros(MB.nfree, 1);
86+
87+acc = zeros(1, nsteps);
88+fhist = zeros(1, nsteps);
89+shellActive = 0;
90+peakDisp = 0;
91+peakXi = 0;
92+inContact = 0;
93+cvel = zeros(1, K);
94+
95+for k = 1:nsteps
96+ % ---- mallet: contact force from the state at step k ----
97+ xi = U(ls, js) - d(tipi) - gap;
98+ if xi > 0
99+ fM1 = KM * xi^alM;
100+ shellActive = 1;
101+ if inContact == 0 && strokeIdx <= K && cvel(strokeIdx) == 0
102+ cvel(strokeIdx) = (d(tipi) - d1(tipi)) / ht; % tip speed at contact
103+ end
104+ inContact = 1;
105+ if xi > peakXi
106+ peakXi = xi;
107+ end
108+ else
109+ fM1 = 0;
110+ inContact = 0;
111+ end
112+ fhist(k) = fM1;
113+ hA = -(uHp - 2 * uH0 + uHm);
114+ hB = -(ht / 2) * (uHp - uHm);
115+ hC = -ht^2 * (ps * uH0 + ((1 - ps) / 2) * (uHp + uHm));
116+ dnew = MB.ASiB * d + MB.ASiC * d1 + MB.vM * hA + MB.vR * hB + MB.vK * hC ...
117+ + (ht^2 * fM1) * MB.tipcol;
118+
119+ % ---- advance the hand one step ----
120+ if phase == 0
121+ if strokeIdx <= K && k * ht >= tlaunch(strokeIdx)
122+ phase = 1;
123+ vtar = vstr(strokeIdx);
124+ end
125+ elseif phase == 1
126+ vH = max(vH - aH * ht, -vtar);
127+ if inContact ~= 0
128+ phase = 3;
129+ holdLeft = holdSteps;
130+ elseif vH <= -vtar
131+ phase = 2;
132+ end
133+ elseif phase == 2
134+ if inContact ~= 0
135+ phase = 3;
136+ holdLeft = holdSteps;
137+ elseif uHp <= -missDepth
138+ phase = 4;
139+ end
140+ elseif phase == 3
141+ holdLeft = holdLeft - 1;
142+ if holdLeft <= 0
143+ phase = 4;
144+ end
145+ elseif phase == 4
146+ vH = min(vH + aHr * ht, 0.5 * vtar);
147+ if vH > 0 && uHp >= -vH^2 / (2 * aHr)
148+ phase = 5;
149+ end
150+ else
151+ vH = vH - aHr * ht;
152+ if vH <= 0 || uHp >= 0
153+ vH = 0;
154+ phase = 0;
155+ strokeIdx = strokeIdx + 1;
156+ end
157+ end
158+ uHm = uH0;
159+ uH0 = uHp;
160+ uHp = uHp + vH * ht;
161+
162+ if shellActive ~= 0
163+ % ---- shell derivatives and von Karman brackets at step k ----
164+ Ur = D1p * U;
165+ Urr = D2p * U;
166+ v1 = Rd .* Ur + R2d .* (U * T2);
167+ v2 = Rd .* (Ur * T1t) - R2d .* (U * T1t);
168+ LapU = Urr + v1;
169+ rhsA = EhRs * LapU - Eh2 * 2 * (Urr .* v1 - v2 .* v2);
170+ rhsA(1:3, :) = 0;
171+ rhsA(Nr:nrt, :) = 0;
172+ Ph3 = pagemtimes(Fi, reshape(rhsA * F, nrt, 1, nth));
173+ Phi = reshape(Ph3, nrt, nth) * Gi;
174+ Pr = D1p * Phi;
175+ Prr = D2p * Phi;
176+ w1 = Rd .* Pr + R2d .* (Phi * T2);
177+ w2 = Rd .* (Pr * T1t) - R2d .* (Phi * T1t);
178+ NL = Prr .* v1 + Urr .* w1 - 2 * w2 .* v2 - (Prr + w1) / Rs;
179+ NL(ls, js) = NL(ls, js) - fM1 * fsc;
180+ NL(1:nlr0 - 1, :) = 0;
181+ NL(Nr + 1:nrt, :) = 0;
182+ NLh = pagemtimes(Ai, reshape((sNL * NL) * F, nrt, 1, nth));
183+ Uh3n = pagemtimes(AiB, Uh3) + pagemtimes(AiC, Uh13) + NLh;
184+ Unew = reshape(Uh3n, nrt, nth) * Gi;
185+
186+ acc(k) = (Unew(lp, jp) - 2 * U(lp, jp) + U1(lp, jp)) / ht^2;
187+ U1 = U;
188+ U = Unew;
189+ Uh13 = Uh3;
190+ Uh3 = Uh3n;
191+ if mod(k, 512) == 0
192+ pd = max(abs(U(:, js)));
193+ if pd > peakDisp
194+ peakDisp = pd;
195+ end
196+ if ~isfinite(pd)
197+ error('cymbal2: simulation diverged at t = %.4f s (reduce strike_velocity, or raise n_theta / oversample)', k * ht);
198+ end
199+ end
200+ end
201+
202+ d1 = d;
203+ d = dnew;
204+end
205+
206+% Contact bookkeeping from the force history
207+inC = fhist > 0;
208+dinfo.peakForceN = max(fhist);
209+dinfo.contactSteps = sum(inC);
210+dinfo.contactMs = 1000 * ht * dinfo.contactSteps;
211+dinfo.contactEvents = sum(diff([0, inC]) > 0);
212+dinfo.peakXiMm = peakXi * 1000;
213+dinfo.peakDispM = peakDisp;
214+dinfo.peakDispThickness = peakDisp / P.thickness_m;
215+dinfo.peakAccel = max(abs(acc));
216+dinfo.nsteps = nsteps;
217+dinfo.simRate = P.sample_rate * P.oversample;
218+dinfo.tstr = tstr;
219+dinfo.vstr = vstr;
220+dinfo.meanV = mean(vstr);
221+dinfo.contactVel = cvel; % tip speed at first contact of each stroke
222+dinfo.strokesPlayed = min(strokeIdx - 1 + (phase > 0), K);
223+end
matlab/cymbal2_system.madded+154−0View file
@@ -0,0 +1,154 @@
1+function S = cymbal2_system(P)
2+% CYMBAL2_SYSTEM Precompute everything the time loop needs. The linear
3+% updates block-diagonalize over azimuthal wavenumber, so the implicit
4+% shell solve and the Airy (stress function) solve are stored as one small
5+% radial matrix per transform column: the time step becomes
6+% Uhat+ = (A^-1 B) Uhat + (A^-1 C) Uhat- + A^-1 (ht^2/rho h) NLhat
7+% with A, B, C the matrices of Samejima (2021) Eq. (20) restricted to one
8+% wavenumber, the boundary-condition rows of Eqs. (22)-(30) built in, and
9+% NL the explicit nonlinear + forcing terms evaluated in physical space.
10+% Rows are equilibrated before inversion; the stored inverses absorb the
11+% scaling so the loop applies them to raw right-hand sides.
12+
13+G = cymbal2_grid(P);
14+nrt = G.nrt;
15+Nr = G.Nr;
16+nth = G.nth;
17+
18+h = P.thickness_m;
19+rhoh = P.density * h;
20+Eh = P.youngs_modulus * h;
21+Rs = P.radius_m^2 / (2 * P.rise_m);
22+ht = 1 / (P.sample_rate * P.oversample);
23+g = P.g_param;
24+phiw = P.phi_param;
25+kap = 0; % set below from cymbal2_perm (needs Dflex)
26+s1 = P.sigma1 * P.radius_m^2;
27+s0 = P.sigma0;
28+
29+MW = P.washer_mass;
30+RW = P.washer_damping;
31+KW = P.washer_stiffness;
32+% Spring term split as phiw/2 on the (k+1, k-1) average and (1 - phiw) on
33+% step k. Note the paper prints the complementary split, which is unstable
34+% for a massless washer (the roots of (1-phi)/2 (z^2+1) + phi z include
35+% |z| > 1 whenever phi > 1/2); with the implicit weight on the average,
36+% phi = 0.7 is unconditionally stable.
37+AW = MW + (ht / 2) * RW + (phiw / 2) * ht^2 * KW;
38+BW = 2 * MW - (1 - phiw) * ht^2 * KW;
39+CW = -MW + (ht / 2) * RW - (phiw / 2) * ht^2 * KW;
40+
41+AiB = zeros(nrt, nrt, nth);
42+AiC = zeros(nrt, nrt, nth);
43+Ai = zeros(nrt, nrt, nth);
44+Fi = zeros(nrt, nrt, nth);
45+nDropped = 0;
46+
47+e1 = zeros(1, nrt); e1(1) = 1;
48+e3 = zeros(1, nrt); e3(3) = 1;
49+I = eye(nrt);
50+
51+for m = 0:nth / 2
52+ lam2 = (2 * cos(m * G.hth) - 2) / G.hth^2;
53+ op = cymbal2_perm(P, G, lam2);
54+ kap = op.Dflex / rhoh;
55+
56+ Am = zeros(nrt, nrt);
57+ Bm = zeros(nrt, nrt);
58+ Cm = zeros(nrt, nrt);
59+ if P.support_clamped ~= 0
60+ l0 = 4;
61+ else
62+ l0 = 3;
63+ end
64+ for l = l0:Nr
65+ Am(l, :) = (1 + s0 * ht) * I(l, :) + g * kap * ht^2 * op.Bih(l, :) - s1 * ht * op.Lap(l, :);
66+ Bm(l, :) = 2 * I(l, :) - (1 - 2 * g) * kap * ht^2 * op.Bih(l, :);
67+ Cm(l, :) = (s0 * ht - 1) * I(l, :) - g * kap * ht^2 * op.Bih(l, :) - s1 * ht * op.Lap(l, :);
68+ end
69+ if P.support_clamped ~= 0
70+ Am(1, :) = e1;
71+ Am(2, :) = e3;
72+ Am(3, :) = G.D1(3, :);
73+ else
74+ % The shear side of the washer balance gets the same implicit split as
75+ % the spring; the paper evaluates it at step k only, which feeds a
76+ % spurious mode on the virtual lines when the washer is massless.
77+ Am(1, :) = op.MomW;
78+ Am(2, :) = AW * e3 - (phiw / 2) * ht^2 * op.ShW;
79+ Bm(2, :) = BW * e3 + (1 - phiw) * ht^2 * op.ShW;
80+ Cm(2, :) = CW * e3 + (phiw / 2) * ht^2 * op.ShW;
81+ end
82+ Am(Nr + 1, :) = op.MomE;
83+ Am(Nr + 2, :) = op.ShE;
84+
85+ % Row equilibration, absorbed into the stored inverses
86+ sc = 1.0 ./ max(abs(Am), [], 2);
87+ Asc = sc .* Am;
88+ Aisc = inv(Asc);
89+ AiBm = Aisc * (sc .* Bm);
90+ AiCm = Aisc * (sc .* Cm);
91+ Aim = Aisc .* sc'; % inv(Am) applied to raw rhs
92+
93+ % On grids too coarse to resolve a wavenumber's inner boundary layer,
94+ % the discretized boundary conditions can hand that wavenumber a
95+ % negative stiffness; such columns are truncated rather than integrated.
96+ % At the default grid none are.
97+ sr = max(abs(eig([AiBm, AiCm; I, zeros(nrt)])));
98+ if sr > 1 + 1e-7
99+ AiBm = zeros(nrt, nrt);
100+ AiCm = zeros(nrt, nrt);
101+ Aim = zeros(nrt, nrt);
102+ nDropped = nDropped + 1 + (m > 0 && m < nth / 2);
103+ end
104+
105+ % Airy solve: biharmonic rows inside, phi = phi_r = 0 at both circles
106+ Fm = zeros(nrt, nrt);
107+ Fm(1, 1) = 1;
108+ Fm(2, :) = G.D1(3, :);
109+ Fm(3, 3) = 1;
110+ for l = 4:Nr - 1
111+ Fm(l, :) = op.Bih(l, :);
112+ end
113+ Fm(Nr, Nr) = 1;
114+ Fm(Nr + 1, :) = G.D1(Nr, :);
115+ Fm(Nr + 2, Nr + 2) = 1;
116+ scf = 1.0 ./ max(abs(Fm), [], 2);
117+ Fim = inv(scf .* Fm) .* scf';
118+
119+ if m == 0
120+ cols = 1;
121+ elseif m == nth / 2
122+ cols = nth;
123+ else
124+ cols = [2 * m, 2 * m + 1];
125+ end
126+ for cc = 1:length(cols)
127+ AiB(:, :, cols(cc)) = AiBm;
128+ AiC(:, :, cols(cc)) = AiCm;
129+ Ai(:, :, cols(cc)) = Aim;
130+ Fi(:, :, cols(cc)) = Fim;
131+ end
132+end
133+
134+if nDropped > nth / 2
135+ error(['cymbal2: the radial grid cannot resolve the inner boundary layer ' ...
136+ '(need roughly washer_radius > 3 (a - c) / (n_radial - 3)); raise n_radial']);
137+end
138+
139+% Strike and pickup nodes (radial fraction 0 at the washer, 1 at the rim)
140+ls = 3 + round(min(1, max(0, P.strike_r)) * (Nr - 3));
141+lp = 3 + round(min(1, max(0, P.pickup_r)) * (Nr - 3));
142+js = 1 + mod(round(P.strike_theta_deg / 360 * nth), nth);
143+jp = 1 + mod(round(P.pickup_theta_deg / 360 * nth), nth);
144+
145+S.G = G;
146+S.AiB = AiB; S.AiC = AiC; S.Ai = Ai; S.Fi = Fi;
147+S.rhoh = rhoh; S.Eh = Eh; S.Rs = Rs; S.ht = ht;
148+S.ls = ls; S.lp = lp; S.js = js; S.jp = jp;
149+S.eomRow0 = 3 + (P.support_clamped ~= 0); % first equation-of-motion row
150+S.nDropped = nDropped; % truncated theta columns
151+S.strikeRadiusM = G.r(ls);
152+S.pickupRadiusM = G.r(lp);
153+S.forceScale = 1 / (G.r(ls) * G.hr * G.hth); % point force -> pressure
154+end
matlab/cymbal2_validate.madded+226−0View file
@@ -0,0 +1,226 @@
1+% Validate the cymbal2 model. There is no external reference
2+% implementation, so the checks are (a) closed-form limits the
3+% discretization must reproduce (the free circular plate, rigid-body
4+% motion on the washer, cantilever statics of the stick), (b) agreement
5+% with the linear-FEM eigenfrequencies published in Samejima (2021) for
6+% the modes that model resolves the same way, (c) structural identities
7+% of the transform and of the assembled update, (d) exact reproducibility,
8+% and (e) pinned end-to-end diagnostics from a reduced configuration,
9+% recorded when the model was built and treated as a regression reference.
10+% Run with: numbl run cymbal2_validate.m
11+
12+fails = 0;
13+P = cymbal2_defaults();
14+
15+% ---- 1. flat limit: the completely free circular plate ----
16+% With the curvature removed and the washer stiffness zeroed, the shell
17+% must reproduce the free plate's frequency equation (the 2 x 2 J_n / I_n
18+% determinant under zero radial moment and zero Kelvin-Kirchhoff shear).
19+% m = 0 is excluded: its limit as the free hole closes is logarithmic, so
20+% a finite hole shifts it far more than the m >= 1 modes.
21+Q = P;
22+Q.rise_m = 1e-9;
23+Q.washer_stiffness = 0;
24+Q.washer_radius_m = 0.010;
25+Q.n_radial = 96;
26+Dp = Q.youngs_modulus * Q.thickness_m^3 / (12 * (1 - Q.poisson_ratio^2));
27+scal = sqrt(Dp / (Q.density * Q.thickness_m)) / Q.radius_m^2 / (2 * pi);
28+errP = 0;
29+for m = 1:4
30+ f = cymbal2_modes(Q, m);
31+ fe = f(f > 0.5);
32+ lam = c2v_plate_roots(m, Q.poisson_ratio);
33+ fc = lam.^2 * scal;
34+ nc = min(3, min(length(fe), length(fc)));
35+ errP = max(errP, max(abs(fe(1:nc) - fc(1:nc)') ./ fc(1:nc)'));
36+end
37+ok = errP < 0.04;
38+fails = fails + (~ok);
39+fprintf('flat plate : max freq err %.1e over m = 1..4 %s\n', errP, c2v_verdict(ok));
40+
41+% ---- 2. rigid limit: bounce and tilt on the washer ----
42+% With the plate stiffened a hundredfold the lowest m = 0 and m = 1 modes
43+% must approach the closed-form rigid-body values on the line spring:
44+% bounce f^2 = KW 2 pi c / (rho h pi (a^2 - c^2)) / (2 pi)^2 and tilt
45+% f^2 = 4 KW c^3 / (rho h (a^4 - c^4)) / (2 pi)^2.
46+Q = P;
47+Q.rise_m = 1e-9;
48+Q.youngs_modulus = P.youngs_modulus * 100;
49+Q.n_radial = 96;
50+a = Q.radius_m; c = Q.washer_radius_m;
51+rhoh = Q.density * Q.thickness_m;
52+KW = Q.washer_stiffness;
53+f0r = sqrt(KW * 2 * pi * c / (rhoh * pi * (a^2 - c^2))) / (2 * pi);
54+f1r = sqrt(4 * KW * c^3 / (rhoh * (a^4 - c^4))) / (2 * pi);
55+fa = cymbal2_modes(Q, 0); fa = fa(fa > 0.01);
56+fb = cymbal2_modes(Q, 1); fb = fb(fb > 0.01);
57+e0 = abs(fa(1) / f0r - 1);
58+e1 = abs(fb(1) / f1r - 1);
59+ok = e0 < 0.03 && e1 < 0.03;
60+fails = fails + (~ok);
61+fprintf('washer : bounce err %.1e, tilt err %.1e vs rigid closed form %s\n', e0, e1, c2v_verdict(ok));
62+
63+% ---- 3. the paper's linear FEM reference, where comparable ----
64+% Samejima's Fig. 9 gives linear-FEM eigenfrequencies of the same shell.
65+% For the clamped center, the (2,1), (3,1), (4,1) modes must agree to a
66+% few percent ((1,1) is excluded: it is controlled by the support model,
67+% where this shell's moment-free hole and their 2-D elements differ).
68+Q = P;
69+Q.support_clamped = 1;
70+ref = [20.1, 43.3, 77.7];
71+got = zeros(1, 3);
72+for m = 2:4
73+ f = cymbal2_modes(Q, m);
74+ f = f(f > 1);
75+ got(m - 1) = f(1);
76+end
77+errF = max(abs(got - ref) ./ ref);
78+ok = errF < 0.06;
79+fails = fails + (~ok);
80+fprintf('paper FEM : clamped (2,1)(3,1)(4,1) = %.1f %.1f %.1f Hz (ref %.1f %.1f %.1f), err %.1e %s\n', ...
81+ got(1), got(2), got(3), ref(1), ref(2), ref(3), errF, c2v_verdict(ok));
82+
83+% ---- 4. transform and update-structure identities ----
84+% The azimuthal transform must be exactly invertible, and the physical
85+% derivative pipeline must agree with the per-wavenumber operators: for a
86+% field cos(m theta) g(r), the physical-space Laplacian must equal the
87+% per-m radial Laplacian applied to g.
88+G = cymbal2_grid(P);
89+eT = max(max(abs(G.F * G.Gi - eye(G.nth))));
90+gr = exp(-((G.r - 0.12) / 0.05).^2);
91+m = 7;
92+Uf = gr * cos(m * G.th');
93+LapPhys = G.D2 * Uf + G.Rd .* (G.D1 * Uf) + G.R2d .* (Uf * G.T2);
94+lam2 = (2 * cos(m * G.hth) - 2) / G.hth^2;
95+op = cymbal2_perm(P, G, lam2);
96+LapM = (op.Lap * gr) * cos(m * G.th');
97+r2 = 2:G.nrt - 1;
98+eL = max(max(abs(LapPhys(r2, :) - LapM(r2, :)))) / max(max(abs(LapM(r2, :))));
99+ok = eT < 1e-10 && eL < 1e-10;
100+fails = fails + (~ok);
101+fprintf('operators : transform inverse %.1e, physical vs per-m Laplacian %.1e %s\n', ...
102+ eT, eL, c2v_verdict(ok));
103+
104+% ---- 5. the assembled update is not amplifying ----
105+% Every azimuthal column's two-step update must have spectral radius at
106+% most 1 (the damping makes physical modes strictly contractive; boundary
107+% constraints contribute zeros). The system builder truncates columns
108+% that fail; at the defaults none may be.
109+S = cymbal2_system(P);
110+nrt = S.G.nrt;
111+srMax = 0;
112+for col = [1, 2, S.G.nth / 2, S.G.nth]
113+ T = [S.AiB(:, :, col), S.AiC(:, :, col); eye(nrt), zeros(nrt)];
114+ srMax = max(srMax, max(abs(eig(T))));
115+end
116+ok = S.nDropped == 0 && srMax <= 1 + 1e-9;
117+fails = fails + (~ok);
118+fprintf('stability : %d truncated columns, sampled spectral radius %.9f %s\n', ...
119+ S.nDropped, srMax, c2v_verdict(ok));
120+
121+% ---- 6. stick statics ----
122+% A static tip load on the clamped-hand stick must deflect by
123+% F L^3 / (3 E I); the rotational-spring hand adds F L^2 / K_R times L.
124+MBc = cymbal2_mallet(setfield(P, 'grip_clamped', 1));
125+MBs = cymbal2_mallet(P);
126+EI = P.stick_youngs * P.stick_inertia_m4;
127+L = P.stick_length_m;
128+Fl = 1.0;
129+% static balance of the update: at a fixed point, (I - A^-1 B - A^-1 C) d
130+% equals A^-1 ht^2 K d, so the deflection under a constant tip load F is
131+% (I - ASiB - ASiC) \ (tipcol * ht^2 * F)
132+htv = MBc.ht;
133+dc = ((eye(size(MBc.ASiB, 1)) - MBc.ASiB - MBc.ASiC) \ (MBc.tipcol * htv^2)) * Fl;
134+ds = ((eye(size(MBs.ASiB, 1)) - MBs.ASiB - MBs.ASiC) \ (MBs.tipcol * htv^2)) * Fl;
135+wc = dc(MBc.tipi);
136+ws = ds(MBs.tipi);
137+wcRef = Fl * L^3 / (3 * EI);
138+wsRef = wcRef + (Fl * L / P.grip_spring) * L; % grip rotation F L / K_R times the arm
139+ec = abs(wc / wcRef - 1);
140+es = abs(ws / wsRef - 1);
141+ok = ec < 1e-6 && es < 1e-6;
142+fails = fails + (~ok);
143+fprintf('stick : clamped-hand tip %.3e (ref %.3e), spring-hand %.3e (ref %.3e) %s\n', ...
144+ wc, wcRef, ws, wsRef, c2v_verdict(ok));
145+
146+% ---- 7. reproducibility and pinned diagnostics ----
147+% Two identical renders of a reduced configuration must agree exactly,
148+% and its diagnostics are pinned. The trajectory of a nonlinear shell is
149+% sensitive, so the acoustic pins carry percent-level tolerances while
150+% the contact pins are tight.
151+Q = P;
152+Q.sample_rate = 44100;
153+Q.oversample = 1;
154+Q.tail_s = 1.5;
155+[y1, d1] = cymbal2_render(Q);
156+[y2, d2] = cymbal2_render(Q);
157+rms1 = sqrt(mean(y1.^2));
158+rms2 = sqrt(mean(y2.^2));
159+same = max(abs(y1 - y2)) == 0;
160+ok = same ...
161+ && d1.contactEvents == d2.contactEvents ...
162+ && abs(d1.peakForceN - 311.14) < 4 ...
163+ && abs(d1.contactVel(1) + 8.223) < 0.08 ...
164+ && abs(rms1 - 0.0973) < 0.010;
165+fails = fails + (~ok);
166+fprintf('pinned : repro %d, peakF %.2f N (ref 311.14), cvel %.3f (ref -8.223), rms %.4f (ref 0.0973) %s\n', ...
167+ same, d1.peakForceN, d1.contactVel(1), rms1, c2v_verdict(ok));
168+
169+% ---- 8. the crash decays ----
170+n = length(y1);
171+fs = Q.sample_rate;
172+head = y1(1:round(0.2 * fs));
173+tail = y1(n - round(0.2 * fs):n);
174+dropE = 10 * log10(sum(tail.^2) / sum(head.^2));
175+ok = dropE < -3 && dropE > -50;
176+fails = fails + (~ok);
177+fprintf('decay : last/first 0.2 s %.1f dB %s\n', dropE, c2v_verdict(ok));
178+
179+if fails > 0
180+ fprintf('%d CHECK(S) FAILED\n', fails);
181+else
182+ fprintf('all checks passed\n');
183+end
184+
185+function lam = c2v_plate_roots(n, nu)
186+% Roots of the free-edge frequency equation of a circular plate
187+lam = [];
188+prev = c2v_plate_det(n, nu, 0.30);
189+for x = 0.30:0.01:20
190+ cur = c2v_plate_det(n, nu, x + 0.01);
191+ if prev * cur < 0
192+ lo = x; hi = x + 0.01;
193+ for it = 1:60
194+ mid = (lo + hi) / 2;
195+ if c2v_plate_det(n, nu, lo) * c2v_plate_det(n, nu, mid) < 0
196+ hi = mid;
197+ else
198+ lo = mid;
199+ end
200+ end
201+ lam = [lam, (lo + hi) / 2];
202+ end
203+ prev = cur;
204+end
205+end
206+
207+function d = c2v_plate_det(n, nu, x)
208+J = besselj(n, x); I = besseli(n, x);
209+Jp = besselj(n - 1, x) - (n / x) * J;
210+Ip = besseli(n - 1, x) - (n / x) * I;
211+Jpp = -Jp / x + (n^2 / x^2 - 1) * J;
212+Ipp = -Ip / x + (n^2 / x^2 + 1) * I;
213+a11 = Jpp + nu * (Jp / x - n^2 * J / x^2);
214+a12 = Ipp + nu * (Ip / x - n^2 * I / x^2);
215+a21 = -Jp - (1 - nu) * (n^2 / x^2) * (Jp - J / x);
216+a22 = Ip - (1 - nu) * (n^2 / x^2) * (Ip - I / x);
217+d = a11 * a22 - a12 * a21;
218+end
219+
220+function s = c2v_verdict(ok)
221+if ok
222+ s = 'ok';
223+else
224+ s = 'MISMATCH';
225+end
226+end
scripts/build-exports.mjsadded+79−0View file
@@ -0,0 +1,79 @@
1+#!/usr/bin/env node
2+// Build exported instrument bundles (format numbl-instrument-v1) for one
3+// instrument repository. Run from the repo root:
4+//
5+// node scripts/build-exports.mjs
6+//
7+// For each export-src/<name>/export.json, the bundle's files are that
8+// directory's .m sources (the drivers, render_note.m required) plus every
9+// matlab/*.m of the repo, and the result is written to exports/<id>.json.
10+// The canonical copy of this script lives in
11+// simulated_instruments/instrument-format; instrument repos carry a copy so
12+// they build standalone.
13+
14+import { mkdirSync, readdirSync, readFileSync, writeFileSync, existsSync } from 'node:fs';
15+import { join } from 'node:path';
16+
17+const ID_RE = /^[a-z0-9][a-z0-9-]*$/;
18+const SRC = 'export-src';
19+const OUT = 'exports';
20+
21+function fail(msg) {
22+ console.error(`build-exports: ${msg}`);
23+ process.exit(1);
24+}
25+
26+function checkKnobs(id, knobs) {
27+ if (!Array.isArray(knobs)) fail(`${id}: 'knobs' must be an array`);
28+ for (const k of knobs) {
29+ if (typeof k.key !== 'string' || typeof k.label !== 'string') fail(`${id}: knobs need key and label`);
30+ if (typeof k.min !== 'number' || typeof k.max !== 'number' || typeof k.default !== 'number' || !(k.min < k.max)) {
31+ fail(`${id}: knob '${k.key}' needs min < max and a default`);
32+ }
33+ if (k.default < k.min || k.default > k.max) fail(`${id}: knob '${k.key}' default outside its range`);
34+ if (k.step !== undefined && !(k.step > 0)) fail(`${id}: knob '${k.key}' step must be > 0`);
35+ if (k.scale !== undefined && k.scale !== 'log') fail(`${id}: knob '${k.key}' has unknown scale`);
36+ if (k.scale === 'log' && k.min <= 0) fail(`${id}: knob '${k.key}' log scale needs min > 0`);
37+ }
38+}
39+
40+function mfilesIn(dir) {
41+ if (!existsSync(dir)) return [];
42+ return readdirSync(dir)
43+ .filter((f) => f.endsWith('.m'))
44+ .sort()
45+ .map((name) => ({ name, content: readFileSync(join(dir, name), 'utf8') }));
46+}
47+
48+if (!existsSync(SRC)) fail(`no ${SRC}/ directory here; run from the instrument repo root`);
49+const shared = mfilesIn('matlab');
50+const dirs = readdirSync(SRC, { withFileTypes: true }).filter((d) => d.isDirectory()).map((d) => d.name).sort();
51+if (dirs.length === 0) fail(`${SRC}/ has no instrument directories`);
52+mkdirSync(OUT, { recursive: true });
53+
54+for (const dir of dirs) {
55+ const metaPath = join(SRC, dir, 'export.json');
56+ if (!existsSync(metaPath)) fail(`${metaPath} is missing`);
57+ const meta = JSON.parse(readFileSync(metaPath, 'utf8'));
58+ const id = meta.id;
59+ if (typeof id !== 'string' || !ID_RE.test(id)) fail(`${metaPath}: 'id' must match ${ID_RE}`);
60+ if (typeof meta.name !== 'string' || meta.name === '') fail(`${id}: missing 'name'`);
61+ checkKnobs(id, meta.knobs ?? []);
62+ const drivers = mfilesIn(join(SRC, dir));
63+ if (!drivers.some((f) => f.name === 'render_note.m')) fail(`${id}: no render_note.m in ${SRC}/${dir}/`);
64+ const names = new Set(drivers.map((f) => f.name));
65+ const clash = shared.find((f) => names.has(f.name));
66+ if (clash) fail(`${id}: driver file '${clash.name}' collides with matlab/${clash.name}`);
67+ const bundle = {
68+ format: 'numbl-instrument-v1',
69+ id,
70+ name: meta.name,
71+ ...(meta.description ? { description: meta.description } : {}),
72+ ...(meta.source ? { source: meta.source } : {}),
73+ knobs: meta.knobs ?? [],
74+ files: [...drivers, ...shared],
75+ };
76+ const outPath = join(OUT, `${id}.json`);
77+ writeFileSync(outPath, JSON.stringify(bundle, null, 2) + '\n');
78+ console.log(`wrote ${outPath} (${bundle.files.length} files, ${bundle.knobs.length} knobs)`);
79+}
webapp/.gitignoreadded+2−0View file
@@ -0,0 +1,2 @@
1+node_modules/
2+dist/
webapp/README.mdadded+26−0View file
@@ -0,0 +1,26 @@
1+# Cymbal2 Simulator — web app
2+
3+Interactive browser implementation of the physical model described in `../README.md`. The algorithm itself is written in MATLAB syntax, in `../matlab/*.m`, and executed client-side by [numbl](https://numbl.org) inside a Web Worker; the TypeScript here is only the UI, the plots, the WAV encoder, and the marshalling of parameters and tensors in and out of the numbl runtime.
4+
5+The signal chain: a prescribed hand stroke drives an Euler-Bernoulli stick with a Hertz-contact mallet head against a von Karman shallow spherical shell on a washer support, integrated step by step in the time domain; the output is the shell's normal acceleration at a pickup point, decimated and normalized.
6+
7+## Usage
8+
9+```
10+npm install
11+npm run dev # development server
12+npm run build # production build in dist/
13+npm run selftest # run ../matlab/cymbal2_validate.m through numbl in Node
14+```
15+
16+The self-test runs the closed-form and reproducibility checks and the pinned diagnostics; the two full renders inside it make it take a few minutes. The MATLAB code can also be run directly: `cd ../matlab && numbl run cymbal2_demo.m` (or `cymbal2_validate.m`), and runs unchanged in MATLAB proper, where it is several times faster.
17+
18+## Where the time goes
19+
20+Everywhere. Unlike the sibling `cymbal` app, whose expensive stage (the eigensolves) could be cached and whose synthesis was an oscillator bank, here the sound is a nonlinear time-domain integration of the whole shell, several hundred thousand steps of coupled solves, and every parameter feeds it. There is therefore no debounced auto-render: the Render preview button runs a cheap pass (no oversampling, tail clamped to 1.5 s; the azimuthal grid is kept, because hard strokes diverge without the full wavenumber ladder), Render full runs the real thing, and both report elapsed time and can be cancelled. Expect tens of seconds for a preview and minutes for a full render in the browser; the only piece the worker caches is the cheap linear mode table that the spectrum plot overlays.
21+
22+## Reading and editing the code in the app
23+
24+The Code button opens the MATLAB source in a CodeMirror editor with MATLAB syntax highlighting. Every `.m` file in `../matlab` is listed, tagged by whether the app executes it or it is a CLI entry point. The sources are editable and an edit takes effect on the next render, so the physics can be changed in place: setting the curvature term `Rs` large in `cymbal2_system.m` flattens the shell into a plate, and zeroing the `Luu` term in `cymbal2_sim.m` removes the stretching nonlinearity, leaving a linear shell that never crashes however hard it is hit.
25+
26+Edits live in memory only. They are deliberately not persisted, so a reload restores the code the app shipped with; "Revert file" and "Revert all" undo them sooner, and a dot marks the files that differ from the shipped source. A fault in edited code is reported with its file and line, the last successful render stays on screen, and reverting recovers.
webapp/index.htmladded+13−0View file
@@ -0,0 +1,13 @@
1+<!doctype html>
2+<html lang="en">
3+ <head>
4+ <meta charset="UTF-8" />
5+ <meta name="viewport" content="width=device-width, initial-scale=1.0" />
6+ <link rel="icon" href="data:image/svg+xml,<svg xmlns=%22http://www.w3.org/2000/svg%22 viewBox=%220 0 100 100%22><text y=%22.9em%22 font-size=%2290%22>🛎️</text></svg>" />
7+ <title>Cymbal2 Simulator</title>
8+ </head>
9+ <body>
10+ <div id="root"></div>
11+ <script type="module" src="/src/main.tsx"></script>
12+ </body>
13+</html>
webapp/package-lock.jsonadded+4172−0View file
This diff is 4,177 lines long and is not shown.
webapp/package.jsonadded+32−0View file
@@ -0,0 +1,32 @@
1+{
2+ "name": "cymbal2-webapp",
3+ "private": true,
4+ "version": "0.1.0",
5+ "type": "module",
6+ "scripts": {
7+ "dev": "vite",
8+ "build": "tsc --noEmit && vite build",
9+ "preview": "vite preview",
10+ "selftest": "tsx src/selftest.ts"
11+ },
12+ "dependencies": {
13+ "@codemirror/commands": "^6.11.0",
14+ "@codemirror/language": "^6.12.4",
15+ "@codemirror/legacy-modes": "^6.5.3",
16+ "@codemirror/state": "^6.7.1",
17+ "@codemirror/theme-one-dark": "^6.1.3",
18+ "@codemirror/view": "^6.43.9",
19+ "codemirror": "^6.0.2",
20+ "numbl": "^0.4.19",
21+ "react": "^18.3.1",
22+ "react-dom": "^18.3.1"
23+ },
24+ "devDependencies": {
25+ "@types/react": "^18.3.5",
26+ "@types/react-dom": "^18.3.0",
27+ "@vitejs/plugin-react": "^4.3.1",
28+ "tsx": "^4.19.0",
29+ "typescript": "^5.5.4",
30+ "vite": "^5.4.3"
31+ }
32+}
webapp/src/App.tsxadded+259−0View file
@@ -0,0 +1,259 @@
1+import { useCallback, useEffect, useRef, useState } from 'react';
2+import { DEFAULTS, loadParams, saveParams, type Cymbal2Params } from './params';
3+import { DEFAULT_MFILES, type MFile } from './mfiles';
4+import { useRenderer } from './hooks/useRenderer';
5+import { encodeWav } from './lib/wav';
6+import ParamPanel from './components/ParamPanel';
7+import CodePanel from './components/CodePanel';
8+import StrikePicker from './components/StrikePicker';
9+import DiagnosticsPanel from './components/DiagnosticsPanel';
10+import WaveformPlot from './components/plots/WaveformPlot';
11+import SpectrumPlot from './components/plots/SpectrumPlot';
12+import SpectrogramPlot from './components/plots/SpectrogramPlot';
13+import ContactPlot from './components/plots/ContactPlot';
14+
15+export default function App() {
16+ const [params, setParams] = useState<Cymbal2Params>(() => loadParams());
17+ // Editable copy of the MATLAB implementation. Deliberately not persisted:
18+ // a reload restores the shipped code.
19+ const [mfiles, setMfiles] = useState<MFile[]>(() => DEFAULT_MFILES.map((f) => ({ ...f })));
20+ const [codeOpen, setCodeOpen] = useState(false);
21+ const { result, error, busy, stage, elapsedS, renderPreview, renderFull, cancel } = useRenderer(
22+ params,
23+ mfiles,
24+ );
25+ const audioCtxRef = useRef<AudioContext | null>(null);
26+
27+ useEffect(() => saveParams(params), [params]);
28+
29+ const onChange = useCallback((key: keyof Cymbal2Params, value: number | boolean) => {
30+ setParams((p) => ({ ...p, [key]: value }));
31+ }, []);
32+
33+ const reset = useCallback(() => setParams({ ...DEFAULTS }), []);
34+
35+ const onChangeFile = useCallback((name: string, source: string) => {
36+ setMfiles((fs) => fs.map((f) => (f.name === name ? { ...f, source } : f)));
37+ }, []);
38+
39+ const onResetFile = useCallback((name: string) => {
40+ setMfiles((fs) =>
41+ fs.map((f) =>
42+ f.name === name
43+ ? { ...f, source: DEFAULT_MFILES.find((d) => d.name === name)?.source ?? f.source }
44+ : f,
45+ ),
46+ );
47+ }, []);
48+
49+ const onResetAllFiles = useCallback(() => setMfiles(DEFAULT_MFILES.map((f) => ({ ...f }))), []);
50+ const codeEdited = mfiles.some(
51+ (f) => f.source !== DEFAULT_MFILES.find((d) => d.name === f.name)?.source,
52+ );
53+
54+ const [looping, setLooping] = useState(false);
55+ const loopSrcRef = useRef<AudioBufferSourceNode | null>(null);
56+
57+ const getAudioCtx = useCallback((sampleRate: number) => {
58+ let ctx = audioCtxRef.current;
59+ if (!ctx || ctx.state === 'closed') {
60+ ctx = new AudioContext({ sampleRate });
61+ audioCtxRef.current = ctx;
62+ }
63+ void ctx.resume();
64+ return ctx;
65+ }, []);
66+
67+ const makeBuffer = useCallback((ctx: AudioContext, r: NonNullable<typeof result>) => {
68+ const buf = ctx.createBuffer(1, r.y.length, r.sampleRate);
69+ buf.copyToChannel(new Float32Array(r.y), 0);
70+ return buf;
71+ }, []);
72+
73+ const play = useCallback(() => {
74+ if (!result) return;
75+ const ctx = getAudioCtx(result.sampleRate);
76+ const src = ctx.createBufferSource();
77+ src.buffer = makeBuffer(ctx, result);
78+ src.connect(ctx.destination);
79+ src.start();
80+ }, [result, getAudioCtx, makeBuffer]);
81+
82+ useEffect(() => {
83+ if (!looping || !result) return;
84+ const ctx = getAudioCtx(result.sampleRate);
85+ const src = ctx.createBufferSource();
86+ src.buffer = makeBuffer(ctx, result);
87+ src.loop = true;
88+ src.connect(ctx.destination);
89+ src.start();
90+ loopSrcRef.current = src;
91+ return () => {
92+ loopSrcRef.current = null;
93+ try {
94+ src.stop();
95+ } catch {
96+ // already stopped
97+ }
98+ src.disconnect();
99+ };
100+ }, [looping, result, getAudioCtx, makeBuffer]);
101+
102+ const download = useCallback(() => {
103+ if (!result) return;
104+ const blob = encodeWav(result.y, result.sampleRate, params.bit_depth === 16 ? 16 : 24);
105+ const url = URL.createObjectURL(blob);
106+ const a = document.createElement('a');
107+ a.href = url;
108+ a.download = result.quality === 'full' ? 'cymbal2.wav' : 'cymbal2_preview.wav';
109+ a.click();
110+ setTimeout(() => URL.revokeObjectURL(url), 5000);
111+ }, [result, params.bit_depth]);
112+
113+ return (
114+ <div className="app">
115+ <div className="topbar">
116+ <h1>Cymbal2 Simulator</h1>
117+ <button
118+ className="primary"
119+ onClick={renderPreview}
120+ disabled={busy}
121+ title="Coarse azimuthal grid, no oversampling, short tail — tens of seconds"
122+ >
123+ Render preview
124+ </button>
125+ <button
126+ className="primary"
127+ onClick={renderFull}
128+ disabled={busy}
129+ title="Full time-domain integration at the set quality — takes minutes in the browser"
130+ >
131+ Render full
132+ </button>
133+ {busy && (
134+ <button onClick={cancel} title="Stop the integration">
135+ ✕ Cancel
136+ </button>
137+ )}
138+ <button
139+ className={`play ${result?.quality === 'preview' ? 'preview' : ''}`}
140+ onClick={play}
141+ disabled={!result}
142+ >
143+ ▶ Play{result?.quality === 'preview' ? ' (preview)' : ''}
144+ </button>
145+ <button
146+ className={`loop ${looping ? 'active' : ''}`}
147+ onClick={() => setLooping((v) => !v)}
148+ disabled={!result && !looping}
149+ title="Repeat the last render continuously"
150+ >
151+ ⟳ Loop{looping ? ' on' : ''}
152+ </button>
153+ <button onClick={download} disabled={!result}>
154+ Download .wav{result?.quality === 'preview' ? ' (preview)' : ''}
155+ </button>
156+ <button onClick={reset}>Reset to defaults</button>
157+ <button
158+ className={codeOpen ? 'code active' : 'code'}
159+ onClick={() => setCodeOpen((v) => !v)}
160+ title="View and edit the MATLAB source that numbl executes"
161+ >
162+ {'{ }'} Code{codeEdited ? ' •' : ''}
163+ </button>
164+ <div className="status">
165+ {error && <span className="badge err">code error</span>}
166+ {result && !busy && (
167+ <span className={`badge ${result.quality}`}>
168+ {result.quality === 'full' ? 'full quality' : 'preview quality'}
169+ </span>
170+ )}
171+ {busy ? (
172+ <span>
173+ {stage}… {elapsedS.toFixed(0)} s (a full render is a real time-domain
174+ integration; expect minutes)
175+ </span>
176+ ) : (
177+ <span>ready — renders are explicit, nothing runs on its own</span>
178+ )}
179+ </div>
180+ </div>
181+ <div className="main">
182+ <div className="sidebar">
183+ <ParamPanel params={params} onChange={onChange} />
184+ </div>
185+ <div className="content">
186+ {error && (
187+ <div className="card wide error-card">
188+ <h2>MATLAB error</h2>
189+ <div>
190+ <strong>
191+ {error.file ? `${error.file}${error.line ? `:${error.line}` : ''}` : 'error'}
192+ </strong>{' '}
193+ {error.message}
194+ </div>
195+ <div className="code-hint">
196+ Showing the last successful render. Open Code to fix the source, or revert it.
197+ </div>
198+ </div>
199+ )}
200+ <div className="card wide">
201+ <h2>Cross-section — drag the strike and pickup points</h2>
202+ <StrikePicker params={params} onChange={onChange} />
203+ </div>
204+ {result && (
205+ <>
206+ <div className="card">
207+ <h2>Diagnostics</h2>
208+ <DiagnosticsPanel result={result} />
209+ </div>
210+ <div className="card">
211+ <h2>Contact force</h2>
212+ <ContactPlot forceT={result.forceT} forceN={result.forceN} />
213+ </div>
214+ <div className="card">
215+ <h2>Waveform</h2>
216+ <WaveformPlot y={result.y} sampleRate={result.sampleRate} />
217+ </div>
218+ <div className="card">
219+ <h2>Spectrum</h2>
220+ <SpectrumPlot
221+ y={result.y}
222+ sampleRate={result.sampleRate}
223+ modeF={result.modeF}
224+ modeN={result.modeM}
225+ />
226+ </div>
227+ <div className="card wide">
228+ <h2>Spectrogram — watch the crash build after the attack</h2>
229+ <SpectrogramPlot y={result.y} sampleRate={result.sampleRate} />
230+ </div>
231+ </>
232+ )}
233+ {!result && !busy && (
234+ <div className="card wide">
235+ <h2>No render yet</h2>
236+ <div className="code-hint">
237+ Press Render preview for a quick coarse pass, or Render full for the real thing.
238+ Unlike the sibling modal-cymbal app, every render here is a nonlinear time-domain
239+ integration of the shell, so nothing renders automatically.
240+ </div>
241+ </div>
242+ )}
243+ </div>
244+ </div>
245+ {codeOpen && (
246+ <CodePanel
247+ mfiles={mfiles}
248+ onChangeFile={onChangeFile}
249+ onResetFile={onResetFile}
250+ onResetAll={onResetAllFiles}
251+ onClose={() => setCodeOpen(false)}
252+ onRenderFull={renderFull}
253+ error={error}
254+ busy={busy}
255+ />
256+ )}
257+ </div>
258+ );
259+}
webapp/src/components/CodeEditor.tsxadded+91−0View file
@@ -0,0 +1,91 @@
1+// CodeMirror 6 editor for one .m file. The parent gives each instance a key
2+// of the file name, so switching files remounts rather than syncing docs.
3+
4+import { useEffect, useRef } from 'react';
5+import { EditorState } from '@codemirror/state';
6+import { EditorView, keymap, lineNumbers, highlightActiveLine } from '@codemirror/view';
7+import { defaultKeymap, history, historyKeymap, indentWithTab } from '@codemirror/commands';
8+import { StreamLanguage, syntaxHighlighting, HighlightStyle } from '@codemirror/language';
9+import type { StreamParser } from '@codemirror/language';
10+import { octave } from '@codemirror/legacy-modes/mode/octave';
11+import { tags as t } from '@lezer/highlight';
12+
13+// The octave mode's keyword list is missing these MATLAB keywords, so promote
14+// them after it has tokenized the word.
15+const EXTRA_KEYWORDS = new Set(['function', 'break', 'parfor', 'spmd']);
16+
17+const matlabMode: StreamParser<unknown> = {
18+ ...octave,
19+ token(stream, state) {
20+ const start = stream.pos;
21+ const tok = octave.token(stream, state);
22+ if (EXTRA_KEYWORDS.has(stream.string.slice(start, stream.pos))) return 'keyword';
23+ return tok;
24+ },
25+};
26+
27+const matlabHighlight = HighlightStyle.define([
28+ { tag: t.keyword, color: '#e8a33d' },
29+ { tag: t.comment, color: '#6b7a86', fontStyle: 'italic' },
30+ { tag: t.string, color: '#a5d6a7' },
31+ { tag: t.number, color: '#4dd0e1' },
32+ { tag: t.operator, color: '#c9a0dc' },
33+ { tag: t.variableName, color: '#d6dae2' },
34+ { tag: t.standard(t.variableName), color: '#82b1ff' }, // builtins: zeros, sqrt, ...
35+ { tag: t.function(t.variableName), color: '#82b1ff' },
36+ { tag: t.atom, color: '#4dd0e1' },
37+ { tag: t.meta, color: '#c9a0dc' },
38+ { tag: t.invalid, color: '#e57373' },
39+]);
40+
41+const editorTheme = EditorView.theme(
42+ {
43+ '&': { color: '#d6dae2', backgroundColor: '#181b20', height: '100%' },
44+ '.cm-content': { fontFamily: 'ui-monospace, SFMono-Regular, Menlo, monospace', fontSize: '12.5px' },
45+ '.cm-gutters': { backgroundColor: '#181b20', color: '#4d5560', border: 'none' },
46+ '.cm-activeLine': { backgroundColor: '#20242b' },
47+ '.cm-cursor': { borderLeftColor: '#e8a33d' },
48+ '&.cm-focused': { outline: 'none' },
49+ '.cm-scroller': { overflow: 'auto' },
50+ '.cm-selectionBackground, ::selection': { backgroundColor: '#31465a' },
51+ '&.cm-focused .cm-selectionBackground': { backgroundColor: '#31465a' },
52+ },
53+ { dark: true },
54+);
55+
56+interface Props {
57+ source: string;
58+ onChange: (source: string) => void;
59+}
60+
61+export default function CodeEditor({ source, onChange }: Props) {
62+ const hostRef = useRef<HTMLDivElement | null>(null);
63+ const onChangeRef = useRef(onChange);
64+ onChangeRef.current = onChange;
65+
66+ useEffect(() => {
67+ const host = hostRef.current!;
68+ const view = new EditorView({
69+ parent: host,
70+ state: EditorState.create({
71+ doc: source,
72+ extensions: [
73+ lineNumbers(),
74+ highlightActiveLine(),
75+ history(),
76+ keymap.of([...defaultKeymap, ...historyKeymap, indentWithTab]),
77+ StreamLanguage.define(matlabMode),
78+ syntaxHighlighting(matlabHighlight),
79+ editorTheme,
80+ EditorView.updateListener.of((u) => {
81+ if (u.docChanged) onChangeRef.current(u.state.doc.toString());
82+ }),
83+ ],
84+ }),
85+ });
86+ return () => view.destroy();
87+ // eslint-disable-next-line react-hooks/exhaustive-deps
88+ }, []);
89+
90+ return <div className="code-editor" ref={hostRef} />;
91+}
webapp/src/components/CodePanel.tsxadded+109−0View file
@@ -0,0 +1,109 @@
1+// Overlay for reading and editing the MATLAB implementation. Edits take
2+// effect on the next debounced render and live only as long as the page.
3+
4+import { useState } from 'react';
5+import { DEFAULT_MFILES, EXECUTED_BY_APP, type MFile } from '../mfiles';
6+import type { RenderError } from '../hooks/useRenderer';
7+import CodeEditor from './CodeEditor';
8+
9+interface Props {
10+ mfiles: MFile[];
11+ onChangeFile: (name: string, source: string) => void;
12+ onResetFile: (name: string) => void;
13+ onResetAll: () => void;
14+ onClose: () => void;
15+ onRenderFull: () => void;
16+ error: RenderError | null;
17+ busy: boolean;
18+}
19+
20+function defaultSourceOf(name: string): string {
21+ return DEFAULT_MFILES.find((f) => f.name === name)?.source ?? '';
22+}
23+
24+export default function CodePanel({
25+ mfiles,
26+ onChangeFile,
27+ onResetFile,
28+ onResetAll,
29+ onClose,
30+ onRenderFull,
31+ error,
32+ busy,
33+}: Props) {
34+ const [selected, setSelected] = useState(
35+ mfiles.find((f) => EXECUTED_BY_APP.has(f.name))?.name ?? mfiles[0]?.name ?? '',
36+ );
37+ const file = mfiles.find((f) => f.name === selected) ?? mfiles[0];
38+ const dirty = (f: MFile) => f.source !== defaultSourceOf(f.name);
39+ const anyDirty = mfiles.some(dirty);
40+
41+ return (
42+ <div className="modal-backdrop" onClick={onClose}>
43+ <div className="code-panel" onClick={(e) => e.stopPropagation()}>
44+ <div className="code-header">
45+ <h2>MATLAB source</h2>
46+ <span className="code-note">
47+ Executed by numbl. Edits re-render as you type and are discarded on reload.
48+ </span>
49+ <button onClick={onRenderFull} disabled={busy}>
50+ Render full quality
51+ </button>
52+ <button onClick={onResetAll} disabled={!anyDirty}>
53+ Revert all
54+ </button>
55+ <button onClick={onClose}>Close</button>
56+ </div>
57+ <div className="code-body">
58+ <div className="code-files">
59+ {mfiles.map((f) => (
60+ <button
61+ key={f.name}
62+ className={f.name === file?.name ? 'active' : ''}
63+ onClick={() => setSelected(f.name)}
64+ >
65+ <span className="fname">
66+ {f.name}
67+ {dirty(f) ? <span className="dot" title="edited" /> : null}
68+ </span>
69+ <span className="ftag">{EXECUTED_BY_APP.has(f.name) ? 'in app' : 'CLI'}</span>
70+ </button>
71+ ))}
72+ </div>
73+ <div className="code-main">
74+ {file && (
75+ <>
76+ <div className="code-toolbar">
77+ <span>{file.name}</span>
78+ {!EXECUTED_BY_APP.has(file.name) && (
79+ <span className="code-hint">
80+ entry point for the numbl CLI; the app does not run this file
81+ </span>
82+ )}
83+ <button
84+ onClick={() => onResetFile(file.name)}
85+ disabled={!dirty(file)}
86+ style={{ marginLeft: 'auto' }}
87+ >
88+ Revert file
89+ </button>
90+ </div>
91+ <CodeEditor
92+ key={file.name}
93+ source={file.source}
94+ onChange={(src) => onChangeFile(file.name, src)}
95+ />
96+ </>
97+ )}
98+ </div>
99+ </div>
100+ {error && (
101+ <div className="code-error">
102+ <strong>{error.file ? `${error.file}${error.line ? `:${error.line}` : ''}` : 'error'}</strong>{' '}
103+ {error.message}
104+ </div>
105+ )}
106+ </div>
107+ </div>
108+ );
109+}
webapp/src/components/DiagnosticsPanel.tsxadded+64−0View file
@@ -0,0 +1,64 @@
1+import type { RenderResponse } from '../worker/renderWorker';
2+
3+export default function DiagnosticsPanel({ result }: { result: RenderResponse }) {
4+ const d = result.diag;
5+ return (
6+ <div className="diag">
7+ <table>
8+ <tbody>
9+ <tr>
10+ <th>Stroke</th>
11+ <td>
12+ {d.strokesPlayed === 1
13+ ? `single, hand ${d.meanV.toFixed(2)} m/s`
14+ : `${d.strokesPlayed} strokes, mean hand ${d.meanV.toFixed(2)} m/s`}
15+ {' '}→ tip {Math.abs(d.contactVel).toFixed(2)} m/s at contact
16+ </td>
17+ </tr>
18+ <tr>
19+ <th>Contact</th>
20+ <td>
21+ {d.peakForceN.toFixed(0)} N peak · {d.peakXiMm.toFixed(2)} mm compression ·{' '}
22+ {d.contactMs.toFixed(1)} ms over {d.contactEvents} event{d.contactEvents === 1 ? '' : 's'}
23+ </td>
24+ </tr>
25+ <tr>
26+ <th>Where</th>
27+ <td>
28+ strike r = {(d.strikeRadiusM * 100).toFixed(1)} cm · pickup r ={' '}
29+ {(d.pickupRadiusM * 100).toFixed(1)} cm
30+ </td>
31+ </tr>
32+ <tr>
33+ <th>Shell</th>
34+ <td>
35+ curvature radius {d.RsM.toFixed(2)} m · peak {d.peakDispMm.toFixed(1)} mm ={' '}
36+ {d.peakDispThickness.toFixed(1)} thicknesses
37+ </td>
38+ </tr>
39+ <tr>
40+ <th>Pickup</th>
41+ <td>{d.peakAccel.toExponential(2)} m/s² peak acceleration</td>
42+ </tr>
43+ <tr>
44+ <th>Integration</th>
45+ <td>
46+ {d.nsteps.toLocaleString()} steps at {(d.simRate / 1000).toFixed(1)} kHz ·{' '}
47+ {d.nTheta} azimuthal points · oversample {d.oversample}
48+ </td>
49+ </tr>
50+ <tr>
51+ <th>Length</th>
52+ <td>
53+ {d.durationS.toFixed(2)} s · {result.quality} · rendered in{' '}
54+ {(d.renderMs / 1000).toFixed(1)} s
55+ </td>
56+ </tr>
57+ </tbody>
58+ </table>
59+ {d.warnings.map((w, i) => (
60+ <div key={i} className="warning">⚠ {w}</div>
61+ ))}
62+ </div>
63+ );
64+}
webapp/src/components/ParamPanel.tsxadded+160−0View file
@@ -0,0 +1,160 @@
1+import { PARAM_GROUPS, type ParamSpec, type Cymbal2Params } from '../params';
2+
3+interface Props {
4+ params: Cymbal2Params;
5+ onChange: (key: keyof Cymbal2Params, value: number | boolean) => void;
6+}
7+
8+const SLIDER_STEPS = 1000;
9+
10+function sliderFromValue(spec: ParamSpec, v: number): number {
11+ if (spec.log) {
12+ const t = (Math.log(v) - Math.log(spec.min)) / (Math.log(spec.max) - Math.log(spec.min));
13+ return Math.round(t * SLIDER_STEPS);
14+ }
15+ return Math.round(((v - spec.min) / (spec.max - spec.min)) * SLIDER_STEPS);
16+}
17+
18+function valueFromSlider(spec: ParamSpec, s: number): number {
19+ const t = s / SLIDER_STEPS;
20+ let v: number;
21+ if (spec.log) {
22+ v = Math.exp(Math.log(spec.min) + t * (Math.log(spec.max) - Math.log(spec.min)));
23+ // keep three significant digits so the displayed number stays readable
24+ const mag = 10 ** Math.floor(Math.log10(Math.abs(v)));
25+ v = Math.round((v / mag) * 100) / 100 * mag;
26+ } else {
27+ v = spec.min + t * (spec.max - spec.min);
28+ if (spec.step) v = Math.round(v / spec.step) * spec.step;
29+ }
30+ if (spec.integer) v = Math.round(v);
31+ return Math.min(spec.max, Math.max(spec.min, v));
32+}
33+
34+function formatValue(v: number): string {
35+ if (v === 0) return '0';
36+ const a = Math.abs(v);
37+ if (a >= 1e5 || a < 1e-3) return v.toExponential(2);
38+ return String(Math.round(v * 1e6) / 1e6);
39+}
40+
41+function ParamRow({ spec, params, onChange }: Props & { spec: ParamSpec }) {
42+ const v = params[spec.key] as number;
43+ return (
44+ <div className="param-row" title={spec.help ?? ''}>
45+ <div className="param-label">
46+ <span className="name">
47+ {spec.label}
48+ {spec.unit ? <span style={{ color: 'var(--muted)' }}> ({spec.unit})</span> : null}
49+ </span>
50+ </div>
51+ <div className="param-controls">
52+ <input
53+ type="range"
54+ min={0}
55+ max={SLIDER_STEPS}
56+ value={sliderFromValue(spec, v)}
57+ onChange={(e) => onChange(spec.key, valueFromSlider(spec, Number(e.target.value)))}
58+ />
59+ <input
60+ type="number"
61+ value={formatValue(v)}
62+ step={spec.step ?? 'any'}
63+ onChange={(e) => {
64+ const nv = Number(e.target.value);
65+ if (Number.isFinite(nv)) onChange(spec.key, spec.integer ? Math.round(nv) : nv);
66+ }}
67+ />
68+ </div>
69+ </div>
70+ );
71+}
72+
73+function BoolRow({
74+ label,
75+ value,
76+ onToggle,
77+ note,
78+}: {
79+ label: string;
80+ value: boolean;
81+ onToggle: (v: boolean) => void;
82+ note?: string;
83+}) {
84+ return (
85+ <div className="param-row">
86+ <div className="param-label">
87+ <span className="name">{label}</span>
88+ </div>
89+ <div className="param-controls">
90+ <input type="checkbox" checked={value} onChange={(e) => onToggle(e.target.checked)} />
91+ {note && <span className="param-help">{note}</span>}
92+ </div>
93+ </div>
94+ );
95+}
96+
97+export default function ParamPanel({ params, onChange }: Props) {
98+ return (
99+ <div>
100+ {PARAM_GROUPS.map((g, gi) => (
101+ <details key={g.name} className="group" open={gi < 2}>
102+ <summary>{g.name}</summary>
103+ <div className="group-body">
104+ {g.specs.map((spec) => (
105+ <ParamRow key={spec.key} spec={spec} params={params} onChange={onChange} />
106+ ))}
107+ {g.name === 'Support (washer)' && (
108+ <BoolRow
109+ label="Rigid clamp"
110+ value={params.support_clamped}
111+ onToggle={(v) => onChange('support_clamped', v)}
112+ note={
113+ params.support_clamped
114+ ? 'clamped center (the earlier literature)'
115+ : 'washer impedance (the paper)'
116+ }
117+ />
118+ )}
119+ {g.name === 'Stick and grip' && (
120+ <BoolRow
121+ label="Clamped hand"
122+ value={params.grip_clamped}
123+ onToggle={(v) => onChange('grip_clamped', v)}
124+ note={params.grip_clamped ? 'tight grip' : 'loose grip (rotational spring)'}
125+ />
126+ )}
127+ {g.name === 'Output' && (
128+ <>
129+ <BoolRow
130+ label="Normalize peak"
131+ value={params.normalize}
132+ onToggle={(v) => onChange('normalize', v)}
133+ note={
134+ params.normalize
135+ ? `scaled to ${params.peak_dbfs} dBFS`
136+ : `absolute: ${params.full_scale_accel.toExponential(1)} m/s² = 0 dBFS`
137+ }
138+ />
139+ <div className="param-row">
140+ <div className="param-label">
141+ <span className="name">WAV bit depth</span>
142+ </div>
143+ <div className="param-controls">
144+ <select
145+ value={params.bit_depth}
146+ onChange={(e) => onChange('bit_depth', Number(e.target.value))}
147+ >
148+ <option value={16}>16-bit</option>
149+ <option value={24}>24-bit</option>
150+ </select>
151+ </div>
152+ </div>
153+ </>
154+ )}
155+ </div>
156+ </details>
157+ ))}
158+ </div>
159+ );
160+}
webapp/src/components/StrikePicker.tsxadded+197−0View file
@@ -0,0 +1,197 @@
1+// Where the mallet lands and where the pickup listens, set by dragging on a
2+// scale cross-section of the shell. The model addresses both by radial
3+// fraction (0 at the washer circle, 1 at the rim); this draws the actual
4+// spherical cap, the washer, the mallet above the contact, and the pickup
5+// point, and converts. The strike sits on the theta = 180 meridian and the
6+// pickup (by default) on theta = 0, so the cross-section shows the strike on
7+// the left half and the pickup on the right.
8+
9+import { useCallback, useRef } from 'react';
10+import type { Cymbal2Params } from '../params';
11+import { useCanvas, drawFrame, PLOT_COLORS } from './plots/canvasUtils';
12+
13+interface Props {
14+ params: Cymbal2Params;
15+ onChange: (key: keyof Cymbal2Params, value: number) => void;
16+}
17+
18+const HEIGHT = 190;
19+
20+export default function StrikePicker({ params, onChange }: Props) {
21+ const dragRef = useRef<'strike' | 'pickup' | null>(null);
22+ const geomRef = useRef({ cx: 0, pw: 100 });
23+
24+ const a = params.radius_m;
25+ const c = params.washer_radius_m;
26+ const rise = params.rise_m;
27+ // spherical cap: z(r) = rise * (1 - r^2 / a^2) (shallow approximation)
28+ const capZ = useCallback((r: number) => rise * (1 - (r * r) / (a * a)), [rise, a]);
29+
30+ const strikeLeft = Math.cos((params.strike_theta_deg * Math.PI) / 180) < 0;
31+
32+ const ref = useCanvas(
33+ (ctx, w, h) => {
34+ drawFrame(ctx, w, h);
35+ const padS = 16;
36+ const pw = (w - 2 * padS) / 2;
37+ const cx = w / 2;
38+ geomRef.current = { cx, pw };
39+ const px = (r: number, left: boolean) => cx + (left ? -1 : 1) * (r / a) * pw;
40+ const baseY = h - 52;
41+ const zScale = (h - 110) / Math.max(1e-4, rise);
42+ const py = (z: number) => baseY - z * zScale;
43+
44+ // centimetre ruler
45+ ctx.font = '10px system-ui';
46+ ctx.textAlign = 'center';
47+ for (let cm = -Math.floor(a * 100); cm <= a * 100 + 1e-9; cm += 5) {
48+ const x = px(Math.abs(cm) / 100, cm < 0);
49+ ctx.strokeStyle = PLOT_COLORS.grid;
50+ ctx.beginPath();
51+ ctx.moveTo(x, baseY + 8);
52+ ctx.lineTo(x, baseY + (cm % 10 === 0 ? 14 : 11));
53+ ctx.stroke();
54+ if (cm % 10 === 0) {
55+ ctx.fillStyle = PLOT_COLORS.axis;
56+ ctx.fillText(`${cm}`, x, baseY + 25);
57+ }
58+ }
59+ ctx.fillStyle = PLOT_COLORS.axis;
60+ ctx.textAlign = 'left';
61+ ctx.fillText('radius (cm)', padS, baseY + 25);
62+
63+ // the stand and washer
64+ ctx.strokeStyle = 'rgba(139, 147, 161, 0.5)';
65+ ctx.setLineDash([3, 3]);
66+ ctx.beginPath();
67+ ctx.moveTo(cx, py(capZ(0)) + 3);
68+ ctx.lineTo(cx, baseY + 6);
69+ ctx.stroke();
70+ ctx.setLineDash([]);
71+ ctx.fillStyle = 'rgba(139, 147, 161, 0.7)';
72+ const wx0 = px(c, true);
73+ const wx1 = px(c, false);
74+ ctx.fillRect(wx0, py(capZ(c)) + 2, wx1 - wx0, 5);
75+ ctx.fillStyle = PLOT_COLORS.axis;
76+ ctx.textAlign = 'center';
77+ ctx.fillText('washer', cx, py(capZ(c)) + 18);
78+
79+ // the shell cap, both halves
80+ ctx.strokeStyle = '#c9a227';
81+ ctx.lineWidth = 2.5;
82+ ctx.beginPath();
83+ const N = 80;
84+ for (let i = 0; i <= N; i++) {
85+ const r = -a + (2 * a * i) / N;
86+ const x = px(Math.abs(r), r < 0);
87+ const y = py(capZ(Math.abs(r)));
88+ if (i === 0) ctx.moveTo(x, y);
89+ else ctx.lineTo(x, y);
90+ }
91+ ctx.stroke();
92+ ctx.lineWidth = 1;
93+
94+ // strike marker with the mallet above it
95+ const rs = c + params.strike_r * (a - c);
96+ const xs = px(rs, strikeLeft);
97+ const ys = py(capZ(rs));
98+ ctx.strokeStyle = PLOT_COLORS.accent;
99+ ctx.lineWidth = 2.5;
100+ ctx.beginPath();
101+ ctx.moveTo(xs + (strikeLeft ? 34 : -34), ys - 44);
102+ ctx.lineTo(xs + (strikeLeft ? 6 : -6), ys - 16);
103+ ctx.stroke();
104+ ctx.fillStyle = PLOT_COLORS.accent;
105+ ctx.beginPath();
106+ ctx.arc(xs, ys - 11, 6, 0, 2 * Math.PI);
107+ ctx.fill();
108+ ctx.font = '11px system-ui';
109+ ctx.textAlign = strikeLeft ? 'left' : 'right';
110+ ctx.fillText(
111+ `strike r = ${(rs * 100).toFixed(1)} cm`,
112+ xs + (strikeLeft ? 12 : -12),
113+ ys - 30,
114+ );
115+
116+ // pickup marker on the other meridian
117+ const rp = c + params.pickup_r * (a - c);
118+ const xp = px(rp, !strikeLeft);
119+ const yp = py(capZ(rp));
120+ ctx.strokeStyle = '#4dd0e1';
121+ ctx.fillStyle = '#4dd0e1';
122+ ctx.beginPath();
123+ ctx.arc(xp, yp, 5, 0, 2 * Math.PI);
124+ ctx.fill();
125+ ctx.beginPath();
126+ ctx.moveTo(xp, yp + 6);
127+ ctx.lineTo(xp, yp + 16);
128+ ctx.stroke();
129+ ctx.textAlign = !strikeLeft ? 'left' : 'right';
130+ ctx.fillText(`pickup r = ${(rp * 100).toFixed(1)} cm`, xp + (!strikeLeft ? 10 : -10), yp + 26);
131+
132+ ctx.fillStyle = PLOT_COLORS.axis;
133+ ctx.font = '10px system-ui';
134+ ctx.textAlign = 'left';
135+ ctx.fillText(
136+ `rise ${(rise * 100).toFixed(1)} cm -> curvature radius ${((a * a) / (2 * rise)).toFixed(2)} m (drawn to its own vertical scale)`,
137+ padS,
138+ 16,
139+ );
140+ },
141+ [params.strike_r, params.pickup_r, a, c, rise, capZ, strikeLeft],
142+ HEIGHT,
143+ );
144+
145+ const fracFromEvent = useCallback(
146+ (e: React.PointerEvent<HTMLCanvasElement>) => {
147+ const rect = e.currentTarget.getBoundingClientRect();
148+ const { cx, pw } = geomRef.current;
149+ const x = e.clientX - rect.left;
150+ const left = x < cx;
151+ const r = (Math.abs(x - cx) / pw) * a;
152+ const frac = Math.min(0.99, Math.max(0.01, (r - c) / (a - c)));
153+ return { frac, left };
154+ },
155+ [a, c],
156+ );
157+
158+ const onDown = useCallback(
159+ (e: React.PointerEvent<HTMLCanvasElement>) => {
160+ const { frac, left } = fracFromEvent(e);
161+ const which = left === strikeLeft ? 'strike' : 'pickup';
162+ dragRef.current = which;
163+ e.currentTarget.setPointerCapture(e.pointerId);
164+ onChange(which === 'strike' ? 'strike_r' : 'pickup_r', frac);
165+ },
166+ [fracFromEvent, onChange, strikeLeft],
167+ );
168+
169+ const onMove = useCallback(
170+ (e: React.PointerEvent<HTMLCanvasElement>) => {
171+ if (!dragRef.current) return;
172+ onChange(dragRef.current === 'strike' ? 'strike_r' : 'pickup_r', fracFromEvent(e).frac);
173+ },
174+ [fracFromEvent, onChange],
175+ );
176+
177+ const onUp = useCallback((e: React.PointerEvent<HTMLCanvasElement>) => {
178+ dragRef.current = null;
179+ e.currentTarget.releasePointerCapture(e.pointerId);
180+ }, []);
181+
182+ return (
183+ <div className="strike-picker">
184+ <canvas
185+ ref={ref}
186+ style={{ cursor: 'ew-resize', touchAction: 'none' }}
187+ onPointerDown={onDown}
188+ onPointerMove={onMove}
189+ onPointerUp={onUp}
190+ onPointerCancel={onUp}
191+ />
192+ <div className="picker-hint">
193+ drag the amber marker (strike, on its meridian) or the cyan one (pickup, opposite)
194+ </div>
195+ </div>
196+ );
197+}
webapp/src/components/plots/ContactPlot.tsxadded+96−0View file
@@ -0,0 +1,96 @@
1+// Contact force against time, zoomed to the strokes: the Hertz contact
2+// between the mallet head and the shell, max-pooled by the worker. The
3+// double peaks of a single stroke are the beam's reflection returning to
4+// the head (the paper's Fig. 12); later bumps are the stick bouncing on the
5+// wobbling shell before the wrist clears it.
6+
7+import { useCanvas, drawFrame, PLOT_COLORS } from './canvasUtils';
8+
9+interface Props {
10+ forceT: Float32Array;
11+ forceN: Float32Array;
12+}
13+
14+export default function ContactPlot({ forceT, forceN }: Props) {
15+ const ref = useCanvas(
16+ (ctx, w, h) => {
17+ drawFrame(ctx, w, h);
18+ const padL = 40;
19+ const padB = 18;
20+ const pw = w - padL - 8;
21+ const ph = h - padB - 10;
22+
23+ let peak = 0;
24+ let tFirst = Infinity;
25+ let tLast = 0;
26+ for (let i = 0; i < forceN.length; i++) {
27+ if (forceN[i] > peak) peak = forceN[i];
28+ if (forceN[i] > 0) {
29+ tFirst = Math.min(tFirst, forceT[i]);
30+ tLast = Math.max(tLast, forceT[i]);
31+ }
32+ }
33+ if (!(peak > 0)) {
34+ ctx.fillStyle = PLOT_COLORS.axis;
35+ ctx.font = '11px system-ui';
36+ ctx.fillText('no contact', padL + 8, 30);
37+ return;
38+ }
39+ const t0 = Math.max(0, tFirst - 0.01);
40+ const t1 = tLast + 0.02;
41+ const fx = (t: number) => padL + ((t - t0) / (t1 - t0)) * pw;
42+
43+ ctx.font = '10px system-ui';
44+ ctx.textAlign = 'right';
45+ for (let g = 0; g <= 4; g++) {
46+ const f = (peak * g) / 4;
47+ const yy = 10 + ph - (g / 4) * ph;
48+ ctx.strokeStyle = PLOT_COLORS.grid;
49+ ctx.beginPath();
50+ ctx.moveTo(padL, yy);
51+ ctx.lineTo(padL + pw, yy);
52+ ctx.stroke();
53+ ctx.fillStyle = PLOT_COLORS.axis;
54+ ctx.fillText(f.toFixed(0), padL - 4, yy + 3);
55+ }
56+ ctx.textAlign = 'center';
57+ const span = t1 - t0;
58+ const steps = [0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1];
59+ const tStep = steps.find((s) => span / s <= 10) ?? 1;
60+ for (let t = Math.ceil(t0 / tStep) * tStep; t <= t1; t += tStep) {
61+ const x = fx(t);
62+ ctx.strokeStyle = PLOT_COLORS.grid;
63+ ctx.beginPath();
64+ ctx.moveTo(x, 10);
65+ ctx.lineTo(x, 10 + ph);
66+ ctx.stroke();
67+ ctx.fillStyle = PLOT_COLORS.axis;
68+ ctx.fillText(tStep >= 0.1 ? `${t.toFixed(1)}s` : `${(t * 1000).toFixed(0)}ms`, x, h - 4);
69+ }
70+
71+ ctx.strokeStyle = PLOT_COLORS.accent;
72+ ctx.lineWidth = 1.4;
73+ ctx.beginPath();
74+ let started = false;
75+ for (let i = 0; i < forceN.length; i++) {
76+ if (forceT[i] < t0 || forceT[i] > t1) continue;
77+ const x = fx(forceT[i]);
78+ const yy = 10 + ph - (forceN[i] / peak) * ph;
79+ if (!started) {
80+ ctx.moveTo(x, yy);
81+ started = true;
82+ } else {
83+ ctx.lineTo(x, yy);
84+ }
85+ }
86+ ctx.stroke();
87+ ctx.lineWidth = 1;
88+ ctx.fillStyle = PLOT_COLORS.axis;
89+ ctx.textAlign = 'left';
90+ ctx.fillText(`peak ${peak.toFixed(0)} N`, padL + 6, 20);
91+ },
92+ [forceT, forceN],
93+ 170,
94+ );
95+ return <canvas ref={ref} />;
96+}
webapp/src/components/plots/SpectrogramPlot.tsxadded+64−0View file
@@ -0,0 +1,64 @@
1+// Time-frequency view of the rendered strike; per-partial decay is visible as
2+// horizontal streaks of different lengths.
3+
4+import { stft } from '../../lib/fft';
5+import { useCanvas, drawFrame, PLOT_COLORS, colormap } from './canvasUtils';
6+
7+interface Props {
8+ y: Float32Array;
9+ sampleRate: number;
10+}
11+
12+export default function SpectrogramPlot({ y, sampleRate }: Props) {
13+ const ref = useCanvas(
14+ (ctx, w, h) => {
15+ drawFrame(ctx, w, h);
16+ const padL = 34;
17+ const padB = 16;
18+ const pw = w - padL - 6;
19+ const ph = h - padB - 6;
20+
21+ const { frames, bins, db, floorDb } = stft(y, 512, 128);
22+
23+ const img = ctx.createImageData(frames, bins);
24+ for (let f = 0; f < frames; f++) {
25+ for (let b = 0; b < bins; b++) {
26+ const v = db[f * bins + b];
27+ const t = Math.min(1, Math.max(0, (v - floorDb) / -floorDb));
28+ const [r, g, bl] = colormap(t);
29+ // flip vertically: low frequency at the bottom
30+ const row = bins - 1 - b;
31+ const idx = (row * frames + f) * 4;
32+ img.data[idx] = r;
33+ img.data[idx + 1] = g;
34+ img.data[idx + 2] = bl;
35+ img.data[idx + 3] = 255;
36+ }
37+ }
38+ const off = new OffscreenCanvas(frames, bins);
39+ const octx = off.getContext('2d')!;
40+ octx.putImageData(img, 0, 0);
41+ ctx.imageSmoothingEnabled = true;
42+ ctx.drawImage(off, padL, 6, pw, ph);
43+
44+ // axes
45+ ctx.fillStyle = PLOT_COLORS.axis;
46+ ctx.font = '10px system-ui';
47+ ctx.textAlign = 'right';
48+ const nyq = sampleRate / 2;
49+ for (let kf = 0; kf <= nyq; kf += 4000) {
50+ const yy = 6 + ph - (kf / nyq) * ph;
51+ ctx.fillText(kf >= 1000 ? `${kf / 1000}k` : `${kf}`, padL - 4, yy + 3);
52+ }
53+ ctx.textAlign = 'center';
54+ const durS = y.length / sampleRate;
55+ for (let t = 0; t <= durS + 1e-9; t += 0.1) {
56+ const x = padL + (t / durS) * pw;
57+ ctx.fillText(`${t.toFixed(1)}s`, x, h - 4);
58+ }
59+ },
60+ [y, sampleRate],
61+ 200,
62+ );
63+ return <canvas ref={ref} />;
64+}
webapp/src/components/plots/SpectrumPlot.tsxadded+114−0View file
@@ -0,0 +1,114 @@
1+// Magnitude spectrum of the rendered stroke, log-frequency axis, with the
2+// retained plate modes overlaid and colored by nodal-diameter order. A
3+// cymbal's spectrum is far denser than the sibling triangle's: below a
4+// kilohertz the partials are still resolvable, above it the modal overlap
5+// exceeds one and the spectrum fills in, which is what a wash is.
6+
7+import { magnitudeSpectrum } from '../../lib/fft';
8+import { useCanvas, drawFrame, PLOT_COLORS, modeColor } from './canvasUtils';
9+
10+interface Props {
11+ y: Float32Array;
12+ sampleRate: number;
13+ modeF: Float32Array;
14+ modeN: Float32Array;
15+}
16+
17+const F_LO = 60;
18+const F_HI = 24000;
19+const DB_RANGE = 80;
20+
21+export default function SpectrumPlot({ y, sampleRate, modeF, modeN }: Props) {
22+ let maxN = 1;
23+ for (let i = 0; i < modeN.length; i++) maxN = Math.max(maxN, modeN[i]);
24+ const ref = useCanvas(
25+ (ctx, w, h) => {
26+ drawFrame(ctx, w, h);
27+ const padL = 34;
28+ const padB = 18;
29+ const pw = w - padL - 6;
30+ const ph = h - padB - 8;
31+ const fx = (f: number) =>
32+ padL + ((Math.log10(f) - Math.log10(F_LO)) / (Math.log10(F_HI) - Math.log10(F_LO))) * pw;
33+
34+ ctx.font = '10px system-ui';
35+ ctx.textAlign = 'center';
36+ for (const f of [100, 200, 500, 1000, 2000, 5000, 10000, 20000]) {
37+ const x = fx(f);
38+ ctx.strokeStyle = PLOT_COLORS.grid;
39+ ctx.beginPath();
40+ ctx.moveTo(x, 8);
41+ ctx.lineTo(x, 8 + ph);
42+ ctx.stroke();
43+ ctx.fillStyle = PLOT_COLORS.axis;
44+ ctx.fillText(f >= 1000 ? `${f / 1000}k` : `${f}`, x, h - 5);
45+ }
46+ ctx.textAlign = 'right';
47+ for (let db = 0; db >= -DB_RANGE; db -= 20) {
48+ const yy = 8 + (-db / DB_RANGE) * ph;
49+ ctx.strokeStyle = PLOT_COLORS.grid;
50+ ctx.beginPath();
51+ ctx.moveTo(padL, yy);
52+ ctx.lineTo(padL + pw, yy);
53+ ctx.stroke();
54+ ctx.fillStyle = PLOT_COLORS.axis;
55+ ctx.fillText(`${db}`, padL - 4, yy + 3);
56+ }
57+
58+ const { mag, size } = magnitudeSpectrum(y, 32768);
59+ let peak = 0;
60+ for (let i = 1; i < mag.length; i++) peak = Math.max(peak, mag[i]);
61+ if (peak <= 0) return;
62+
63+ // mode lines first so the trace draws on top
64+ for (let i = 0; i < modeF.length; i++) {
65+ const f = modeF[i];
66+ if (f < F_LO || f > F_HI) continue;
67+ ctx.strokeStyle = modeColor(modeN[i], maxN, 0.35);
68+ ctx.beginPath();
69+ ctx.moveTo(fx(f), 8 + ph * 0.45);
70+ ctx.lineTo(fx(f), 8 + ph);
71+ ctx.stroke();
72+ }
73+
74+ ctx.strokeStyle = PLOT_COLORS.trace;
75+ ctx.beginPath();
76+ let started = false;
77+ for (let px = 0; px <= pw; px++) {
78+ const f = 10 ** (Math.log10(F_LO) + (px / pw) * (Math.log10(F_HI) - Math.log10(F_LO)));
79+ const bin = Math.min(mag.length - 1, Math.round((f / sampleRate) * size));
80+ const fNext =
81+ 10 ** (Math.log10(F_LO) + ((px + 1) / pw) * (Math.log10(F_HI) - Math.log10(F_LO)));
82+ const binNext = Math.min(
83+ mag.length - 1,
84+ Math.max(bin + 1, Math.round((fNext / sampleRate) * size)),
85+ );
86+ let m = 0;
87+ for (let b = bin; b < binNext; b++) m = Math.max(m, mag[b]);
88+ const db = 20 * Math.log10(m / peak + 1e-12);
89+ const yy = 8 + Math.min(ph, (-db / DB_RANGE) * ph);
90+ if (!started) {
91+ ctx.moveTo(padL + px, yy);
92+ started = true;
93+ } else {
94+ ctx.lineTo(padL + px, yy);
95+ }
96+ }
97+ ctx.stroke();
98+
99+ ctx.textAlign = 'left';
100+ ctx.fillStyle = PLOT_COLORS.axis;
101+ ctx.fillText('mode lines colored by nodal diameters, 0', padL + 4, 18);
102+ const x0 = padL + 200;
103+ for (let n = 0; n <= maxN; n++) {
104+ ctx.fillStyle = modeColor(n, maxN, 0.95);
105+ ctx.fillRect(x0 + (n / maxN) * 70, 11, 70 / (maxN + 1) + 1, 8);
106+ }
107+ ctx.fillStyle = PLOT_COLORS.axis;
108+ ctx.fillText(`${maxN}`, x0 + 74, 18);
109+ },
110+ [y, sampleRate, modeF, modeN, maxN],
111+ 200,
112+ );
113+ return <canvas ref={ref} />;
114+}
webapp/src/components/plots/WaveformPlot.tsxadded+66−0View file
@@ -0,0 +1,66 @@
1+import { useCanvas, drawFrame, PLOT_COLORS } from './canvasUtils';
2+
3+interface Props {
4+ y: Float32Array;
5+ sampleRate: number;
6+}
7+
8+export default function WaveformPlot({ y, sampleRate }: Props) {
9+ const ref = useCanvas(
10+ (ctx, w, h) => {
11+ drawFrame(ctx, w, h);
12+ const padL = 34;
13+ const padB = 16;
14+ const pw = w - padL - 6;
15+ const ph = h - padB - 6;
16+ const midY = 6 + ph / 2;
17+
18+ ctx.strokeStyle = PLOT_COLORS.grid;
19+ ctx.beginPath();
20+ ctx.moveTo(padL, midY);
21+ ctx.lineTo(padL + pw, midY);
22+ ctx.stroke();
23+
24+ // time gridlines, step chosen so about ten fit
25+ const durS = y.length / sampleRate;
26+ const steps = [0.05, 0.1, 0.2, 0.5, 1, 2];
27+ const tStep = steps.find((s) => durS / s <= 12) ?? 5;
28+ ctx.fillStyle = PLOT_COLORS.axis;
29+ ctx.font = '10px system-ui';
30+ ctx.textAlign = 'center';
31+ for (let t = 0; t <= durS + 1e-9; t += tStep) {
32+ const x = padL + (t / durS) * pw;
33+ ctx.strokeStyle = PLOT_COLORS.grid;
34+ ctx.beginPath();
35+ ctx.moveTo(x, 6);
36+ ctx.lineTo(x, 6 + ph);
37+ ctx.stroke();
38+ ctx.fillText(`${t.toFixed(1)}s`, x, h - 4);
39+ }
40+ ctx.textAlign = 'right';
41+ ctx.fillText('+1', padL - 4, 12);
42+ ctx.fillText('-1', padL - 4, 6 + ph);
43+
44+ // min/max envelope per pixel column
45+ ctx.strokeStyle = PLOT_COLORS.trace;
46+ ctx.beginPath();
47+ const perPx = y.length / pw;
48+ for (let px = 0; px < pw; px++) {
49+ const i0 = Math.floor(px * perPx);
50+ const i1 = Math.min(y.length, Math.max(i0 + 1, Math.floor((px + 1) * perPx)));
51+ let mn = Infinity, mx = -Infinity;
52+ for (let i = i0; i < i1; i++) {
53+ if (y[i] < mn) mn = y[i];
54+ if (y[i] > mx) mx = y[i];
55+ }
56+ if (!Number.isFinite(mn)) continue;
57+ ctx.moveTo(padL + px + 0.5, midY - mx * (ph / 2));
58+ ctx.lineTo(padL + px + 0.5, midY - mn * (ph / 2) + 0.5);
59+ }
60+ ctx.stroke();
61+ },
62+ [y, sampleRate],
63+ 170,
64+ );
65+ return <canvas ref={ref} />;
66+}
webapp/src/components/plots/canvasUtils.tsadded+77−0View file
@@ -0,0 +1,77 @@
1+import { useEffect, useRef } from 'react';
2+
3+// Canvas hook: resizes to element width at devicePixelRatio and calls draw.
4+export function useCanvas(
5+ draw: (ctx: CanvasRenderingContext2D, w: number, h: number) => void,
6+ deps: unknown[],
7+ height: number,
8+) {
9+ const ref = useRef<HTMLCanvasElement | null>(null);
10+
11+ useEffect(() => {
12+ const canvas = ref.current;
13+ if (!canvas) return;
14+ const render = () => {
15+ const dpr = window.devicePixelRatio || 1;
16+ const w = canvas.clientWidth;
17+ canvas.width = Math.max(1, Math.round(w * dpr));
18+ canvas.height = Math.max(1, Math.round(height * dpr));
19+ canvas.style.height = `${height}px`;
20+ const ctx = canvas.getContext('2d')!;
21+ ctx.setTransform(dpr, 0, 0, dpr, 0, 0);
22+ ctx.clearRect(0, 0, w, height);
23+ draw(ctx, w, height);
24+ };
25+ render();
26+ const ro = new ResizeObserver(render);
27+ ro.observe(canvas);
28+ return () => ro.disconnect();
29+ // eslint-disable-next-line react-hooks/exhaustive-deps
30+ }, deps);
31+
32+ return ref;
33+}
34+
35+export const PLOT_COLORS = {
36+ bg: '#181b20',
37+ grid: '#2a2f38',
38+ axis: '#8b93a1',
39+ trace: '#4dd0e1',
40+ accent: '#e8a33d',
41+ good: '#6fbf73',
42+ warn: '#e57373',
43+};
44+
45+export function drawFrame(ctx: CanvasRenderingContext2D, w: number, h: number) {
46+ ctx.fillStyle = PLOT_COLORS.bg;
47+ ctx.fillRect(0, 0, w, h);
48+}
49+
50+// Viridis-like colormap for the spectrogram, t in [0,1].
51+export function colormap(t: number): [number, number, number] {
52+ const stops: [number, number, number][] = [
53+ [13, 8, 60], [84, 39, 143], [140, 41, 129], [196, 66, 90],
54+ [237, 105, 37], [251, 155, 6], [247, 209, 61], [252, 255, 164],
55+ ];
56+ const x = Math.min(0.9999, Math.max(0, t)) * (stops.length - 1);
57+ const i = Math.floor(x);
58+ const f = x - i;
59+ const a = stops[i];
60+ const b = stops[i + 1];
61+ return [
62+ Math.round(a[0] + f * (b[0] - a[0])),
63+ Math.round(a[1] + f * (b[1] - a[1])),
64+ Math.round(a[2] + f * (b[2] - a[2])),
65+ ];
66+}
67+
68+// Color a mode by its number of nodal diameters: cyan at n = 0 through amber
69+// at the highest order retained. Used by every plot that shows modes, so the
70+// orders read the same everywhere.
71+export function modeColor(n: number, maxN: number, alpha = 1): string {
72+ const t = Math.min(1, Math.max(0, n / Math.max(1, maxN)));
73+ const r = Math.round(77 + t * (232 - 77));
74+ const g = Math.round(208 + t * (163 - 208));
75+ const b = Math.round(225 + t * (61 - 225));
76+ return `rgba(${r}, ${g}, ${b}, ${alpha})`;
77+}
webapp/src/hooks/useRenderer.tsadded+113−0View file
@@ -0,0 +1,113 @@
1+// Hook owning the render worker. Unlike the sibling cymbal app there is no
2+// debounced auto-preview: a cymbal2 render is a full time-domain integration
3+// and takes from tens of seconds (preview) to minutes (full quality) in the
4+// browser, so renders are explicit. A request in flight can be cancelled,
5+// which tears the worker down and builds a fresh one.
6+
7+import { useCallback, useEffect, useRef, useState } from 'react';
8+import type { Cymbal2Params } from '../params';
9+import type { MFile } from '../mfiles';
10+import type { ErrorMessage, ProgressMessage, RenderResponse } from '../worker/renderWorker';
11+
12+export interface RenderError {
13+ message: string;
14+ file?: string;
15+ line?: number;
16+}
17+
18+export interface RenderState {
19+ result: RenderResponse | null;
20+ error: RenderError | null;
21+ busy: boolean;
22+ stage: string;
23+ elapsedS: number;
24+}
25+
26+export function useRenderer(params: Cymbal2Params, mfiles: MFile[]) {
27+ const workerRef = useRef<Worker | null>(null);
28+ const nextId = useRef(1);
29+ const inFlight = useRef<number | null>(null);
30+ const startedAt = useRef(0);
31+ const [state, setState] = useState<RenderState>({
32+ result: null,
33+ error: null,
34+ busy: false,
35+ stage: '',
36+ elapsedS: 0,
37+ });
38+
39+ const attachHandlers = useCallback((worker: Worker) => {
40+ worker.onmessage = (ev: MessageEvent<RenderResponse | ProgressMessage | ErrorMessage>) => {
41+ const msg = ev.data;
42+ if (msg.type === 'progress') {
43+ setState((s) => ({ ...s, stage: msg.stage }));
44+ return;
45+ }
46+ inFlight.current = null;
47+ if (msg.type === 'error') {
48+ setState((s) => ({
49+ ...s,
50+ error: { message: msg.message, file: msg.file, line: msg.line },
51+ busy: false,
52+ stage: '',
53+ }));
54+ return;
55+ }
56+ setState((s) => ({ ...s, result: msg, error: null, busy: false, stage: '' }));
57+ };
58+ }, []);
59+
60+ const makeWorker = useCallback(() => {
61+ const worker = new Worker(new URL('../worker/renderWorker.ts', import.meta.url), {
62+ type: 'module',
63+ });
64+ attachHandlers(worker);
65+ workerRef.current = worker;
66+ return worker;
67+ }, [attachHandlers]);
68+
69+ useEffect(() => {
70+ const worker = makeWorker();
71+ return () => {
72+ worker.terminate();
73+ workerRef.current = null;
74+ inFlight.current = null;
75+ };
76+ }, [makeWorker]);
77+
78+ // elapsed-time ticker while a render is in flight
79+ useEffect(() => {
80+ if (!state.busy) return;
81+ const t = setInterval(
82+ () => setState((s) => ({ ...s, elapsedS: (performance.now() - startedAt.current) / 1000 })),
83+ 250,
84+ );
85+ return () => clearInterval(t);
86+ }, [state.busy]);
87+
88+ const render = useCallback(
89+ (quality: 'preview' | 'full') => {
90+ const worker = workerRef.current;
91+ if (!worker || inFlight.current !== null) return;
92+ const id = nextId.current++;
93+ inFlight.current = id;
94+ startedAt.current = performance.now();
95+ setState((s) => ({ ...s, busy: true, stage: 'starting', elapsedS: 0 }));
96+ worker.postMessage({ type: 'render', id, params, mfiles, quality });
97+ },
98+ [params, mfiles],
99+ );
100+
101+ const cancel = useCallback(() => {
102+ const worker = workerRef.current;
103+ if (worker) worker.terminate();
104+ inFlight.current = null;
105+ makeWorker();
106+ setState((s) => ({ ...s, busy: false, stage: '', elapsedS: 0 }));
107+ }, [makeWorker]);
108+
109+ const renderPreview = useCallback(() => render('preview'), [render]);
110+ const renderFull = useCallback(() => render('full'), [render]);
111+
112+ return { ...state, renderPreview, renderFull, cancel };
113+}
webapp/src/lib/fft.tsadded+91−0View file
@@ -0,0 +1,91 @@
1+// Minimal iterative radix-2 complex FFT, in-place on separate re/im arrays.
2+
3+export function fft(re: Float64Array, im: Float64Array): void {
4+ const n = re.length;
5+ if ((n & (n - 1)) !== 0) throw new Error('FFT length must be a power of 2');
6+ // bit-reversal permutation
7+ for (let i = 1, j = 0; i < n; i++) {
8+ let bit = n >> 1;
9+ for (; j & bit; bit >>= 1) j ^= bit;
10+ j ^= bit;
11+ if (i < j) {
12+ const tr = re[i]; re[i] = re[j]; re[j] = tr;
13+ const ti = im[i]; im[i] = im[j]; im[j] = ti;
14+ }
15+ }
16+ for (let len = 2; len <= n; len <<= 1) {
17+ const ang = (-2 * Math.PI) / len;
18+ const wr = Math.cos(ang);
19+ const wi = Math.sin(ang);
20+ for (let i = 0; i < n; i += len) {
21+ let cwr = 1, cwi = 0;
22+ for (let k = 0; k < len / 2; k++) {
23+ const ur = re[i + k], ui = im[i + k];
24+ const vr = re[i + k + len / 2] * cwr - im[i + k + len / 2] * cwi;
25+ const vi = re[i + k + len / 2] * cwi + im[i + k + len / 2] * cwr;
26+ re[i + k] = ur + vr;
27+ im[i + k] = ui + vi;
28+ re[i + k + len / 2] = ur - vr;
29+ im[i + k + len / 2] = ui - vi;
30+ const nwr = cwr * wr - cwi * wi;
31+ cwi = cwr * wi + cwi * wr;
32+ cwr = nwr;
33+ }
34+ }
35+ }
36+}
37+
38+// Magnitude spectrum of a real signal with a Hann window, zero-padded to the
39+// next power of two (at least minSize). Returns bins 0..N/2.
40+export function magnitudeSpectrum(x: Float32Array | Float64Array, minSize = 16384): {
41+ mag: Float64Array;
42+ size: number;
43+} {
44+ let size = minSize;
45+ while (size < x.length) size <<= 1;
46+ const re = new Float64Array(size);
47+ const im = new Float64Array(size);
48+ const n = x.length;
49+ for (let i = 0; i < n; i++) {
50+ const w = 0.5 * (1 - Math.cos((2 * Math.PI * i) / (n - 1)));
51+ re[i] = x[i] * w;
52+ }
53+ fft(re, im);
54+ const half = size / 2;
55+ const mag = new Float64Array(half + 1);
56+ for (let i = 0; i <= half; i++) mag[i] = Math.hypot(re[i], im[i]);
57+ return { mag, size };
58+}
59+
60+// Simple STFT magnitude in dB for the spectrogram, frames x bins.
61+export function stft(
62+ x: Float32Array,
63+ win = 512,
64+ hop = 128,
65+): { frames: number; bins: number; db: Float32Array; floorDb: number } {
66+ const nFrames = Math.max(1, Math.floor((x.length - win) / hop) + 1);
67+ const bins = win / 2 + 1;
68+ const db = new Float32Array(nFrames * bins);
69+ const re = new Float64Array(win);
70+ const im = new Float64Array(win);
71+ const hann = new Float64Array(win);
72+ for (let i = 0; i < win; i++) hann[i] = 0.5 * (1 - Math.cos((2 * Math.PI * i) / (win - 1)));
73+ let maxDb = -Infinity;
74+ for (let f = 0; f < nFrames; f++) {
75+ const off = f * hop;
76+ for (let i = 0; i < win; i++) {
77+ re[i] = (off + i < x.length ? x[off + i] : 0) * hann[i];
78+ im[i] = 0;
79+ }
80+ fft(re, im);
81+ for (let b = 0; b < bins; b++) {
82+ const m = Math.hypot(re[b], im[b]);
83+ const v = 20 * Math.log10(m + 1e-12);
84+ db[f * bins + b] = v;
85+ if (v > maxDb) maxDb = v;
86+ }
87+ }
88+ // normalize so 0 dB is the frame-bin maximum
89+ for (let i = 0; i < db.length; i++) db[i] -= maxDb;
90+ return { frames: nFrames, bins, db, floorDb: -90 };
91+}
webapp/src/lib/wav.tsadded+40−0View file
@@ -0,0 +1,40 @@
1+// Mono PCM WAV encoding (16- or 24-bit).
2+
3+export function encodeWav(y: Float32Array, sampleRate: number, bitDepth: 16 | 24): Blob {
4+ const bytesPerSample = bitDepth / 8;
5+ const dataSize = y.length * bytesPerSample;
6+ const buf = new ArrayBuffer(44 + dataSize);
7+ const dv = new DataView(buf);
8+ const writeStr = (off: number, s: string) => {
9+ for (let i = 0; i < s.length; i++) dv.setUint8(off + i, s.charCodeAt(i));
10+ };
11+ writeStr(0, 'RIFF');
12+ dv.setUint32(4, 36 + dataSize, true);
13+ writeStr(8, 'WAVE');
14+ writeStr(12, 'fmt ');
15+ dv.setUint32(16, 16, true);
16+ dv.setUint16(20, 1, true); // PCM
17+ dv.setUint16(22, 1, true); // mono
18+ dv.setUint32(24, sampleRate, true);
19+ dv.setUint32(28, sampleRate * bytesPerSample, true);
20+ dv.setUint16(32, bytesPerSample, true);
21+ dv.setUint16(34, bitDepth, true);
22+ writeStr(36, 'data');
23+ dv.setUint32(40, dataSize, true);
24+
25+ let off = 44;
26+ for (let i = 0; i < y.length; i++) {
27+ const v = Math.max(-1, Math.min(1, y[i]));
28+ if (bitDepth === 16) {
29+ dv.setInt16(off, Math.round(v * 32767), true);
30+ off += 2;
31+ } else {
32+ const iv = Math.round(v * 8388607);
33+ dv.setUint8(off, iv & 0xff);
34+ dv.setUint8(off + 1, (iv >> 8) & 0xff);
35+ dv.setUint8(off + 2, (iv >> 16) & 0xff);
36+ off += 3;
37+ }
38+ }
39+ return new Blob([buf], { type: 'audio/wav' });
40+}
webapp/src/main.tsxadded+10−0View file
@@ -0,0 +1,10 @@
1+import React from 'react';
2+import { createRoot } from 'react-dom/client';
3+import App from './App';
4+import './style.css';
5+
6+createRoot(document.getElementById('root')!).render(
7+ <React.StrictMode>
8+ <App />
9+ </React.StrictMode>,
10+);
webapp/src/mfiles.tsadded+44−0View file
@@ -0,0 +1,44 @@
1+// The MATLAB implementation, loaded from ../matlab at build time. The app
2+// keeps an editable in-memory copy of these sources (see App.tsx); edits feed
3+// straight into the render worker and are deliberately not persisted, so a
4+// reload restores the shipped code.
5+
6+export interface MFile {
7+ name: string;
8+ source: string;
9+}
10+
11+const modules = import.meta.glob('../../matlab/*.m', {
12+ eager: true,
13+ query: '?raw',
14+ import: 'default',
15+}) as Record<string, string>;
16+
17+export const DEFAULT_MFILES: MFile[] = Object.entries(modules)
18+ .map(([path, source]) => ({ name: path.split('/').pop()!, source }))
19+ .sort((a, b) => a.name.localeCompare(b.name));
20+
21+// Files the app itself executes. The others are entry points for running the
22+// model from the numbl CLI, and editing them changes nothing in the browser.
23+export const EXECUTED_BY_APP = new Set([
24+ 'cymbal2_grid.m',
25+ 'cymbal2_perm.m',
26+ 'cymbal2_system.m',
27+ 'cymbal2_mallet.m',
28+ 'cymbal2_sim.m',
29+ 'cymbal2_output.m',
30+ 'cymbal2_modes.m',
31+]);
32+
33+// Cheap content hash (FNV-1a) so the worker can key its caches on the code.
34+export function hashMFiles(files: MFile[]): string {
35+ let h = 0x811c9dc5;
36+ for (const f of files) {
37+ const s = f.name + ' ' + f.source + '';
38+ for (let i = 0; i < s.length; i++) {
39+ h ^= s.charCodeAt(i);
40+ h = Math.imul(h, 0x01000193);
41+ }
42+ }
43+ return (h >>> 0).toString(36);
44+}
webapp/src/params.tsadded+274−0View file
@@ -0,0 +1,274 @@
1+// Parameter definitions for the cymbal2 model, mirroring cymbal2_defaults.m.
2+// Unlike the sibling cymbal app there are no recomputation tiers: the
3+// expensive stage is the time integration itself, and everything feeds it.
4+// The render duration is not a parameter: it is
5+// (strike_count - 1) / strike_rate_hz + tail_s, plus the first stroke's
6+// travel time.
7+
8+export interface Cymbal2Params {
9+ // Shell geometry
10+ radius_m: number;
11+ rise_m: number;
12+ thickness_m: number;
13+ // Material
14+ youngs_modulus: number;
15+ poisson_ratio: number;
16+ density: number;
17+ // Damping
18+ sigma0: number;
19+ sigma1: number;
20+ // Support
21+ washer_radius_m: number;
22+ washer_stiffness: number;
23+ washer_mass: number;
24+ washer_damping: number;
25+ support_clamped: boolean;
26+ // Mallet head
27+ head_radius_m: number;
28+ head_mass_kg: number;
29+ head_stiffness: number;
30+ head_alpha: number;
31+ // Mallet stick
32+ stick_length_m: number;
33+ stick_area_m2: number;
34+ stick_youngs: number;
35+ stick_density: number;
36+ stick_inertia_m4: number;
37+ stick_damping_rate: number;
38+ n_beam: number;
39+ // Grip
40+ grip_clamped: boolean;
41+ grip_spring: number;
42+ // Strike gesture
43+ strike_velocity: number;
44+ strike_accel: number;
45+ retract_accel: number;
46+ strike_hold_ms: number;
47+ strike_gap_m: number;
48+ strike_r: number;
49+ strike_theta_deg: number;
50+ strike_count: number;
51+ strike_rate_hz: number;
52+ velocity_spread: number;
53+ seed: number;
54+ // Pickup
55+ pickup_r: number;
56+ pickup_theta_deg: number;
57+ // Discretization
58+ n_radial: number;
59+ n_theta: number;
60+ g_param: number;
61+ phi_param: number;
62+ // Rendering
63+ sample_rate: number;
64+ oversample: number;
65+ tail_s: number;
66+ pre_delay_ms: number;
67+ normalize: boolean;
68+ peak_dbfs: number;
69+ full_scale_accel: number;
70+ bit_depth: number;
71+}
72+
73+export const DEFAULTS: Cymbal2Params = {
74+ radius_m: 0.23,
75+ rise_m: 0.04,
76+ thickness_m: 0.001,
77+ youngs_modulus: 1.1e11,
78+ poisson_ratio: 0.33,
79+ density: 8860.0,
80+ sigma0: 1.34,
81+ sigma1: 1.2e-3,
82+ washer_radius_m: 0.029,
83+ washer_stiffness: 1.65e5,
84+ washer_mass: 0.0,
85+ washer_damping: 0.0,
86+ support_clamped: false,
87+ head_radius_m: 0.025,
88+ head_mass_kg: 0.028,
89+ head_stiffness: 1.6e8,
90+ head_alpha: 2.54,
91+ stick_length_m: 0.39,
92+ stick_area_m2: 1.43139e-4,
93+ stick_youngs: 8.8e9,
94+ stick_density: 400.0,
95+ stick_inertia_m4: 1.63044e-9,
96+ stick_damping_rate: 10.0,
97+ n_beam: 20,
98+ grip_clamped: false,
99+ grip_spring: 140.0,
100+ strike_velocity: 4.0,
101+ strike_accel: 250.0,
102+ retract_accel: 2500.0,
103+ strike_hold_ms: 0.5,
104+ strike_gap_m: -1,
105+ strike_r: 0.6,
106+ strike_theta_deg: 180.0,
107+ strike_count: 1,
108+ strike_rate_hz: 4.0,
109+ velocity_spread: 0.0,
110+ seed: 1,
111+ pickup_r: 0.6,
112+ pickup_theta_deg: 0.0,
113+ n_radial: 24,
114+ n_theta: 60,
115+ g_param: 0.2498,
116+ phi_param: 0.7,
117+ sample_rate: 48000,
118+ oversample: 2,
119+ tail_s: 6.0,
120+ pre_delay_ms: 3.0,
121+ normalize: true,
122+ peak_dbfs: -3.0,
123+ full_scale_accel: 4.0e4,
124+ bit_depth: 24,
125+};
126+
127+// Keys that determine the linear mode table (cymbal2_modes); the cheap modes
128+// stage is cached on these.
129+export const MODES_KEYS: (keyof Cymbal2Params)[] = [
130+ 'radius_m', 'rise_m', 'thickness_m',
131+ 'youngs_modulus', 'poisson_ratio', 'density',
132+ 'washer_radius_m', 'washer_stiffness', 'support_clamped',
133+ 'n_radial', 'n_theta',
134+];
135+
136+export interface ParamSpec {
137+ key: keyof Cymbal2Params;
138+ label: string;
139+ unit?: string;
140+ min: number;
141+ max: number;
142+ step?: number;
143+ log?: boolean;
144+ integer?: boolean;
145+ help?: string;
146+}
147+
148+export interface ParamGroup {
149+ name: string;
150+ specs: ParamSpec[];
151+}
152+
153+export const PARAM_GROUPS: ParamGroup[] = [
154+ {
155+ name: 'Stroke',
156+ specs: [
157+ { key: 'strike_velocity', label: 'Hand speed', unit: 'm/s', min: 0.25, max: 7.0, step: 0.05, help: 'Peak speed of the hand. The compliant stick whips, so the tip arrives about 1.3x faster. 0.5 taps, 2 begins to wash, 4 and up crash. Above about 7 the explicit nonlinear integration diverges, which is also where the shallow-shell physics stops being valid; the renderer reports it rather than emitting garbage.' },
158+ { key: 'strike_r', label: 'Contact position', min: 0.05, max: 0.99, step: 0.005, help: 'Radial fraction from the washer (0) to the rim (1); draggable on the cross-section above. Near the rim drives the high nodal-diameter modes; near the washer rings sparser.' },
159+ { key: 'strike_accel', label: 'Downstroke accel', unit: 'm/s²', min: 50, max: 4000, log: true, help: 'How hard the arm accelerates. Faster strokes give the compliant stick less time to follow, so the tip whips more.' },
160+ { key: 'retract_accel', label: 'Retract accel', unit: 'm/s²', min: 200, max: 20000, log: true, help: 'The wrist snap pulling the stick back after contact. Too slow and the stick shoves the wobbling shell and re-hits it.' },
161+ { key: 'strike_hold_ms', label: 'Hold', unit: 'ms', min: 0, max: 20, step: 0.5, help: 'Dwell between first contact and the retract. Long holds choke the shell against the head.' },
162+ ],
163+ },
164+ {
165+ name: 'Repeats',
166+ specs: [
167+ { key: 'strike_count', label: 'Strokes', min: 1, max: 16, step: 1, integer: true, help: 'The hand returns to rest between strokes, which the paper\'s constant-acceleration hand cannot do.' },
168+ { key: 'strike_rate_hz', label: 'Rate', unit: 'Hz', min: 0.5, max: 12, step: 0.25, help: 'Strokes per second. A stroke whose predecessor has not finished retracting is skipped.' },
169+ { key: 'velocity_spread', label: 'Velocity spread', min: 0, max: 0.7, step: 0.01, help: 'Relative scatter of the per-stroke hand speeds.' },
170+ { key: 'seed', label: 'Random seed', min: 1, max: 9999, step: 1, integer: true, help: 'Redraws the velocity scatter; the same seed always reproduces it.' },
171+ ],
172+ },
173+ {
174+ name: 'Support (washer)',
175+ specs: [
176+ { key: 'washer_stiffness', label: 'Washer stiffness', unit: 'N/m/m', min: 1e4, max: 1e7, log: true, help: 'Spring per unit length of the mounting circle; 1.65e5 is the paper\'s measured felt washer. It sets the bounce (about 13 Hz) and tilt (about 4 Hz) of the whole shell on its mount, and through the boundary condition the low audible modes too.' },
177+ { key: 'washer_damping', label: 'Washer damping', unit: 'N s/m/m', min: 0, max: 500, step: 1, help: 'Dashpot per unit length alongside the spring. The paper sets 0; a real felt sleeve has some.' },
178+ { key: 'washer_mass', label: 'Washer mass', unit: 'kg/m', min: 0, max: 1, step: 0.005, help: 'Mass per unit length moving with the mounting circle. The paper sets 0.' },
179+ { key: 'washer_radius_m', label: 'Washer radius', unit: 'm', min: 0.015, max: 0.06, step: 0.001, help: 'Radius of the mounting circle. The grid needs it above roughly 3 radial spacings; the builder refuses configurations it cannot resolve.' },
180+ ],
181+ },
182+ {
183+ name: 'Mallet',
184+ specs: [
185+ { key: 'head_mass_kg', label: 'Head mass', unit: 'kg', min: 0.002, max: 0.2, log: true, help: 'The 28 g head of the paper\'s measured mallet. Heavier darkens and lengthens the contact.' },
186+ { key: 'head_stiffness', label: 'Head stiffness', unit: 'N m⁻ᵅ', min: 1e6, max: 1e11, log: true, help: 'Coefficient of the Hertz-type contact f = K delta^alpha. Softer rolls the attack off; a wood tip is orders of magnitude stiffer than felt.' },
187+ { key: 'head_alpha', label: 'Contact exponent', min: 1.2, max: 3.5, step: 0.02, help: 'Hertz theory gives 1.5 for a sphere on a plate; the paper\'s measured mallet is 2.54, felt being strongly hardening.' },
188+ { key: 'head_radius_m', label: 'Head radius', unit: 'm', min: 0.005, max: 0.05, step: 0.001 },
189+ ],
190+ },
191+ {
192+ name: 'Stick and grip',
193+ specs: [
194+ { key: 'grip_spring', label: 'Grip stiffness', unit: 'N m/rad', min: 10, max: 5000, log: true, help: 'Rotational spring of the loose grip. Softer lets the stick pivot freely and whip; the audible effect is on the attack, as the paper\'s grip-tightness experiment found.' },
195+ { key: 'stick_youngs', label: 'Stick modulus', unit: 'Pa', min: 1e9, max: 1e11, log: true, help: 'Wood is about 9e9. A stiffer stick follows the hand faster and doubles the contact peaks differently.' },
196+ { key: 'stick_length_m', label: 'Stick length', unit: 'm', min: 0.2, max: 0.6, step: 0.005 },
197+ { key: 'stick_density', label: 'Stick density', unit: 'kg/m³', min: 200, max: 1200, step: 10 },
198+ { key: 'stick_damping_rate', label: 'Stick damping', unit: 's⁻¹', min: 0, max: 100, step: 1, help: 'Mass-proportional damping of the stick; the paper uses 10.' },
199+ ],
200+ },
201+ {
202+ name: 'Shell geometry',
203+ specs: [
204+ { key: 'radius_m', label: 'Radius', unit: 'm', min: 0.12, max: 0.35, step: 0.001, help: '0.23 is the paper\'s 18-inch crash. Mode frequencies scale roughly as thickness / radius².' },
205+ { key: 'rise_m', label: 'Rise', unit: 'm', min: 0.005, max: 0.09, step: 0.001, help: 'Height of the cap; sets the radius of curvature Rs = a²/(2 rise). More curvature stiffens the low modes through the membrane coupling and raises the ring frequency.' },
206+ { key: 'thickness_m', label: 'Thickness', unit: 'm', min: 0.0004, max: 0.003, step: 0.00005, help: 'Uniform; the shallow-shell formulation cannot carry a taper. Thin crashes, thick rides. The nonlinearity strengthens as displacement over thickness.' },
207+ ],
208+ },
209+ {
210+ name: 'Material and damping',
211+ specs: [
212+ { key: 'density', label: 'Density', unit: 'kg/m³', min: 2000, max: 12000, log: true, help: 'The paper\'s measured bronze is 8860.' },
213+ { key: 'youngs_modulus', label: "Young's modulus", unit: 'Pa', min: 3e10, max: 3e11, log: true },
214+ { key: 'poisson_ratio', label: 'Poisson ratio', min: 0.0, max: 0.45, step: 0.01 },
215+ { key: 'sigma0', label: 'Damping σ₀', unit: 's⁻¹', min: 0.05, max: 10, log: true, help: 'Frequency-independent decay rate: T60 = 6.91 / (2 σ₀) ... more precisely 3 ln10 / σ₀ at the bottom of the band. The paper fits 1.34.' },
216+ { key: 'sigma1', label: 'Damping σ₁', min: 0, max: 0.02, step: 0.0001, help: 'High-frequency damping, growing as the squared wavenumber. The paper fits 1.2e-3.' },
217+ ],
218+ },
219+ {
220+ name: 'Pickup',
221+ specs: [
222+ { key: 'pickup_r', label: 'Pickup position', min: 0.05, max: 0.99, step: 0.005, help: 'Radial fraction of the pickup point, draggable on the cross-section. The output is the shell\'s normal acceleration there: a contact pickup, not a microphone.' },
223+ { key: 'pickup_theta_deg', label: 'Pickup azimuth', unit: '°', min: 0, max: 359, step: 1, help: 'The strike is at 180°. On top of it you hear the contact hardest.' },
224+ ],
225+ },
226+ {
227+ name: 'Output',
228+ specs: [
229+ { key: 'tail_s', label: 'Tail', unit: 's', min: 0.5, max: 12.0, step: 0.1, help: 'Ring time rendered after the last stroke. The low modes carry T60 near 5 s at the default damping.' },
230+ { key: 'pre_delay_ms', label: 'Pre-delay', unit: 'ms', min: 0, max: 50, step: 0.5 },
231+ { key: 'peak_dbfs', label: 'Peak level', unit: 'dBFS', min: -24, max: 0, step: 0.5, help: 'Used when normalize is on.' },
232+ { key: 'full_scale_accel', label: 'Full scale', unit: 'm/s²', min: 1e2, max: 1e6, log: true, help: 'Used when normalize is off: this acceleration maps to 0 dBFS.' },
233+ ],
234+ },
235+ {
236+ name: 'Advanced / quality',
237+ specs: [
238+ { key: 'n_radial', label: 'Radial intervals', min: 20, max: 48, step: 1, integer: true, help: 'The paper\'s grid uses 24. The washer circle needs about 3 radial spacings inside it; too-coarse grids are refused. Cost is roughly linear.' },
239+ { key: 'n_theta', label: 'Azimuthal points', min: 24, max: 120, step: 2, integer: true, help: 'The paper uses 60. Sets the highest nodal-diameter order carried; the crash\'s top octaves live there, and hard strokes need the full ladder or the cascade piles up at the cutoff and diverges. Cost is roughly linear.' },
240+ { key: 'oversample', label: 'Oversample', min: 1, max: 4, step: 1, integer: true, help: 'Simulation rate multiple of the output rate. 2 matches the paper\'s 88.2 kHz; hard strokes need it for stability of the explicit nonlinearity.' },
241+ { key: 'g_param', label: 'Scheme g', min: 0.24, max: 0.3, step: 0.0002, help: 'Free parameter of the implicit shell scheme; the paper uses 0.2498. At 0.25 and above the linear part is unconditionally stable.' },
242+ { key: 'phi_param', label: 'Scheme φ', min: 0.5, max: 1.0, step: 0.01, help: 'Implicit weight of the washer and beam schemes.' },
243+ { key: 'n_beam', label: 'Beam elements', min: 8, max: 40, step: 1, integer: true },
244+ { key: 'strike_gap_m', label: 'Rest clearance', unit: 'm', min: -1, max: 0.5, step: 0.005, help: 'Head-to-shell clearance at rest. Nonpositive sizes it to the stroke so contact arrives just after full speed, which keeps the tip whip consistent.' },
245+ ],
246+ },
247+];
248+
249+const STORAGE_KEY = 'cymbal2-webapp-params-v1';
250+
251+export function loadParams(): Cymbal2Params {
252+ try {
253+ const raw = localStorage.getItem(STORAGE_KEY);
254+ if (!raw) return { ...DEFAULTS };
255+ const parsed = JSON.parse(raw);
256+ const out = { ...DEFAULTS };
257+ for (const k of Object.keys(DEFAULTS) as (keyof Cymbal2Params)[]) {
258+ if (typeof parsed[k] === typeof DEFAULTS[k]) {
259+ (out as Record<string, unknown>)[k] = parsed[k];
260+ }
261+ }
262+ return out;
263+ } catch {
264+ return { ...DEFAULTS };
265+ }
266+}
267+
268+export function saveParams(p: Cymbal2Params): void {
269+ try {
270+ localStorage.setItem(STORAGE_KEY, JSON.stringify(p));
271+ } catch {
272+ // storage unavailable; persistence is a convenience only
273+ }
274+}
webapp/src/selftest.tsadded+26−0View file
@@ -0,0 +1,26 @@
1+// Runs the MATLAB validation script (../matlab/cymbal2_validate.m) through
2+// numbl in Node: closed-form limits, the paper's linear-FEM reference,
3+// stability and structural identities, exact reproducibility, and the pinned
4+// end-to-end diagnostics. The two full renders inside make this take a few
5+// minutes. Run with: npm run selftest
6+
7+import { readFileSync, readdirSync } from 'fs';
8+import { join, dirname } from 'path';
9+import { fileURLToPath } from 'url';
10+import { executeCode } from 'numbl';
11+
12+const matlabDir = join(dirname(fileURLToPath(import.meta.url)), '..', '..', 'matlab');
13+const workspaceFiles = readdirSync(matlabDir)
14+ .filter((f) => f.endsWith('.m'))
15+ .map((f) => ({ name: f, source: readFileSync(join(matlabDir, f), 'utf8') }));
16+
17+const main = workspaceFiles.find((f) => f.name === 'cymbal2_validate.m')!;
18+const t0 = performance.now();
19+const result = executeCode(main.source, { displayResults: false }, workspaceFiles, 'cymbal2_validate.m');
20+const text = result.output.join('');
21+process.stdout.write(text);
22+console.log(`total time : ${((performance.now() - t0) / 1000).toFixed(2)} s`);
23+
24+if (text.includes('MISMATCH') || text.includes('FAILED')) {
25+ process.exit(1);
26+}
webapp/src/style.cssadded+368−0View file
@@ -0,0 +1,368 @@
1+:root {
2+ --bg: #14161a;
3+ --panel: #1d2026;
4+ --panel2: #23272f;
5+ --border: #31363f;
6+ --text: #d6dae2;
7+ --muted: #8b93a1;
8+ --accent: #e8a33d;
9+ --accent2: #4dd0e1;
10+ --good: #6fbf73;
11+ --warn: #e57373;
12+ font-family: system-ui, -apple-system, 'Segoe UI', Roboto, sans-serif;
13+}
14+
15+* { box-sizing: border-box; }
16+
17+html, body, #root {
18+ margin: 0;
19+ height: 100%;
20+ background: var(--bg);
21+ color: var(--text);
22+}
23+
24+.app {
25+ display: flex;
26+ flex-direction: column;
27+ height: 100vh;
28+}
29+
30+.topbar {
31+ display: flex;
32+ align-items: center;
33+ gap: 10px;
34+ padding: 8px 14px;
35+ background: var(--panel);
36+ border-bottom: 1px solid var(--border);
37+ flex-wrap: wrap;
38+}
39+
40+.topbar h1 {
41+ font-size: 16px;
42+ margin: 0 12px 0 0;
43+ font-weight: 600;
44+ color: var(--accent);
45+ white-space: nowrap;
46+}
47+
48+button {
49+ background: var(--panel2);
50+ color: var(--text);
51+ border: 1px solid var(--border);
52+ border-radius: 6px;
53+ padding: 6px 14px;
54+ font-size: 13px;
55+ cursor: pointer;
56+}
57+button:hover { border-color: var(--muted); }
58+button:disabled { opacity: 0.45; cursor: default; }
59+button.primary {
60+ background: var(--accent);
61+ color: #1a1204;
62+ border-color: var(--accent);
63+ font-weight: 600;
64+}
65+button.play {
66+ background: var(--good);
67+ color: #06210a;
68+ border-color: var(--good);
69+ font-weight: 600;
70+}
71+button.play.preview {
72+ background: transparent;
73+ color: var(--accent);
74+ border-color: var(--accent);
75+}
76+button.loop.active {
77+ background: var(--accent2);
78+ color: #062a2f;
79+ border-color: var(--accent2);
80+ font-weight: 600;
81+}
82+button.code { font-family: ui-monospace, SFMono-Regular, Menlo, monospace; }
83+button.code.active { border-color: var(--accent2); color: var(--accent2); }
84+
85+.status {
86+ margin-left: auto;
87+ font-size: 12px;
88+ color: var(--muted);
89+ display: flex;
90+ align-items: center;
91+ gap: 8px;
92+ min-width: 180px;
93+ justify-content: flex-end;
94+}
95+
96+.badge {
97+ font-size: 11px;
98+ border-radius: 10px;
99+ padding: 2px 8px;
100+ border: 1px solid var(--border);
101+}
102+.badge.preview { color: var(--accent); border-color: var(--accent); }
103+.badge.full { color: var(--good); border-color: var(--good); }
104+.badge.err { color: var(--warn); border-color: var(--warn); }
105+
106+.progressbar {
107+ width: 90px;
108+ height: 6px;
109+ background: var(--panel2);
110+ border-radius: 3px;
111+ overflow: hidden;
112+}
113+.progressbar > div {
114+ height: 100%;
115+ background: var(--accent2);
116+ transition: width 0.1s linear;
117+}
118+
119+.main {
120+ display: flex;
121+ flex: 1;
122+ min-height: 0;
123+}
124+
125+.sidebar {
126+ width: 340px;
127+ min-width: 340px;
128+ overflow-y: auto;
129+ background: var(--panel);
130+ border-right: 1px solid var(--border);
131+ padding: 8px;
132+}
133+
134+.content {
135+ flex: 1;
136+ overflow-y: auto;
137+ padding: 10px;
138+ display: grid;
139+ grid-template-columns: 1fr 1fr;
140+ gap: 10px;
141+ align-content: start;
142+}
143+
144+.card {
145+ background: var(--panel);
146+ border: 1px solid var(--border);
147+ border-radius: 8px;
148+ padding: 8px 10px 10px;
149+ min-width: 0;
150+}
151+.card.wide { grid-column: 1 / -1; }
152+.card h2 {
153+ font-size: 12px;
154+ text-transform: uppercase;
155+ letter-spacing: 0.06em;
156+ color: var(--muted);
157+ margin: 0 0 6px;
158+ font-weight: 600;
159+}
160+.card canvas { display: block; width: 100%; }
161+
162+.scene3d { height: 380px; position: relative; }
163+.scene3d canvas { height: 100%; border-radius: 6px; }
164+.scene-hint {
165+ position: absolute;
166+ bottom: 10px;
167+ right: 14px;
168+ font-size: 11px;
169+ color: var(--muted);
170+ pointer-events: none;
171+}
172+
173+/* Parameter panel */
174+details.group {
175+ border: 1px solid var(--border);
176+ border-radius: 6px;
177+ margin-bottom: 6px;
178+ background: var(--panel2);
179+}
180+details.group summary {
181+ cursor: pointer;
182+ padding: 6px 10px;
183+ font-size: 13px;
184+ font-weight: 600;
185+ user-select: none;
186+}
187+details.group[open] summary { border-bottom: 1px solid var(--border); }
188+.group-body { padding: 6px 10px 8px; }
189+
190+.param-row { margin-bottom: 8px; }
191+.param-label {
192+ display: flex;
193+ justify-content: space-between;
194+ font-size: 12px;
195+ margin-bottom: 2px;
196+ align-items: baseline;
197+}
198+.param-label .name { color: var(--text); }
199+.param-controls { display: flex; gap: 8px; align-items: center; }
200+.param-controls input[type='range'] { flex: 1; accent-color: var(--accent2); }
201+.param-controls input[type='number'] {
202+ width: 88px;
203+ background: var(--bg);
204+ color: var(--text);
205+ border: 1px solid var(--border);
206+ border-radius: 4px;
207+ padding: 3px 5px;
208+ font-size: 12px;
209+}
210+.param-controls select, .param-controls input[type='checkbox'] {
211+ background: var(--bg);
212+ color: var(--text);
213+ border: 1px solid var(--border);
214+ border-radius: 4px;
215+ font-size: 12px;
216+}
217+.param-help { font-size: 11px; color: var(--muted); margin-top: 1px; }
218+
219+/* Diagnostics */
220+.diag table { border-collapse: collapse; font-size: 12px; width: 100%; }
221+.diag td, .diag th {
222+ padding: 2px 8px 2px 0;
223+ text-align: left;
224+ color: var(--text);
225+ font-weight: 400;
226+}
227+.diag th { color: var(--muted); white-space: nowrap; }
228+.diag .partials { font-variant-numeric: tabular-nums; }
229+.warning { color: var(--warn); font-size: 12px; margin-top: 4px; }
230+.cachenote { color: var(--muted); font-size: 11px; margin-top: 6px; }
231+
232+.mode-list { display: flex; gap: 6px; flex-wrap: wrap; margin-bottom: 6px; }
233+.mode-list button { padding: 3px 10px; font-size: 12px; }
234+.mode-list button.active { border-color: var(--accent2); color: var(--accent2); }
235+
236+/* Code panel */
237+.modal-backdrop {
238+ position: fixed;
239+ inset: 0;
240+ background: rgba(8, 10, 13, 0.72);
241+ display: flex;
242+ align-items: center;
243+ justify-content: center;
244+ z-index: 50;
245+}
246+
247+.code-panel {
248+ width: min(1180px, 94vw);
249+ height: min(860px, 92vh);
250+ background: var(--panel);
251+ border: 1px solid var(--border);
252+ border-radius: 10px;
253+ display: flex;
254+ flex-direction: column;
255+ overflow: hidden;
256+}
257+
258+.code-header {
259+ display: flex;
260+ align-items: center;
261+ gap: 10px;
262+ padding: 9px 12px;
263+ border-bottom: 1px solid var(--border);
264+}
265+.code-header h2 {
266+ margin: 0;
267+ font-size: 13px;
268+ text-transform: uppercase;
269+ letter-spacing: 0.06em;
270+ color: var(--accent);
271+}
272+.code-header button:last-child { margin-left: 0; }
273+.code-note {
274+ font-size: 11.5px;
275+ color: var(--muted);
276+ margin-right: auto;
277+}
278+
279+.code-body { display: flex; flex: 1; min-height: 0; }
280+
281+.code-files {
282+ width: 236px;
283+ min-width: 236px;
284+ border-right: 1px solid var(--border);
285+ overflow-y: auto;
286+ padding: 6px;
287+}
288+.code-files button {
289+ display: flex;
290+ align-items: center;
291+ justify-content: space-between;
292+ gap: 6px;
293+ width: 100%;
294+ text-align: left;
295+ background: transparent;
296+ border: 1px solid transparent;
297+ border-radius: 5px;
298+ padding: 5px 8px;
299+ margin-bottom: 2px;
300+ font-family: ui-monospace, SFMono-Regular, Menlo, monospace;
301+ font-size: 12px;
302+}
303+.code-files button:hover { background: var(--panel2); }
304+.code-files button.active {
305+ background: var(--panel2);
306+ border-color: var(--accent2);
307+ color: var(--accent2);
308+}
309+.code-files .fname { display: flex; align-items: center; gap: 5px; }
310+.code-files .dot {
311+ width: 6px;
312+ height: 6px;
313+ border-radius: 50%;
314+ background: var(--accent);
315+}
316+.code-files .ftag {
317+ font-size: 9.5px;
318+ color: var(--muted);
319+ text-transform: uppercase;
320+ letter-spacing: 0.05em;
321+ font-family: system-ui, sans-serif;
322+}
323+
324+.code-main { flex: 1; min-width: 0; display: flex; flex-direction: column; }
325+.code-toolbar {
326+ display: flex;
327+ align-items: center;
328+ gap: 10px;
329+ padding: 6px 10px;
330+ border-bottom: 1px solid var(--border);
331+ font-family: ui-monospace, SFMono-Regular, Menlo, monospace;
332+ font-size: 12px;
333+ color: var(--text);
334+}
335+.code-hint {
336+ font-size: 11px;
337+ color: var(--muted);
338+ font-family: system-ui, sans-serif;
339+}
340+.code-editor { flex: 1; min-height: 0; overflow: hidden; }
341+
342+.code-error {
343+ border-top: 1px solid var(--warn);
344+ background: #2a1a1a;
345+ color: #ffb4a9;
346+ padding: 7px 12px;
347+ font-size: 12px;
348+ font-family: ui-monospace, SFMono-Regular, Menlo, monospace;
349+}
350+
351+.error-card { border-color: var(--warn); }
352+.error-card h2 { color: var(--warn); }
353+.error-card div { font-size: 12px; }
354+.error-card strong { font-family: ui-monospace, SFMono-Regular, Menlo, monospace; }
355+
356+@media (max-width: 1100px) {
357+ .content { grid-template-columns: 1fr; }
358+}
359+
360+/* Strike picker */
361+.strike-picker { position: relative; }
362+.strike-picker canvas { display: block; width: 100%; }
363+.picker-hint {
364+ font-size: 11px;
365+ color: var(--muted);
366+ text-align: right;
367+ margin-top: 2px;
368+}
webapp/src/vite-env.d.tsadded+1−0View file
@@ -0,0 +1 @@
1+/// <reference types="vite/client" />
webapp/src/worker/renderWorker.tsadded+308−0View file
@@ -0,0 +1,308 @@
1+// Web Worker that renders the cymbal2 stroke by executing the MATLAB
2+// implementation in ../matlab through numbl. Unlike the sibling cymbal app
3+// there is no cheap stage to cache around: the cost IS the time integration,
4+// one long executeCode call, so requests carry a quality tier
5+// preview -> oversample 1, tail clamped to 1.5 s (the azimuthal grid is
6+// kept: hard strokes need the full wavenumber band to carry
7+// the cascade, and coarsening it destabilizes the integration)
8+// full -> the parameters as given
9+// and the UI shows elapsed time rather than a meaningful progress fraction.
10+// A cheap linear-modes stage (cymbal2_modes per wavenumber) is cached on the
11+// geometry so the spectrum plot can overlay the mode table immediately.
12+
13+import { executeCode, type RuntimeValue, type RuntimeTensor } from 'numbl';
14+import { hashMFiles, type MFile } from '../mfiles';
15+import { MODES_KEYS, type Cymbal2Params } from '../params';
16+
17+export interface RenderRequest {
18+ type: 'render';
19+ id: number;
20+ params: Cymbal2Params;
21+ quality: 'preview' | 'full';
22+ mfiles: MFile[];
23+}
24+
25+export interface RenderResponse {
26+ type: 'result';
27+ id: number;
28+ quality: 'preview' | 'full';
29+ y: Float32Array; // final signal, [-1, 1]
30+ sampleRate: number;
31+ forceT: Float32Array; // contact force history (max-pooled), s
32+ forceN: Float32Array; // and its values, N
33+ modeF: Float32Array; // linear mode frequencies, Hz
34+ modeM: Float32Array; // nodal diameters of each
35+ diag: {
36+ strokesPlayed: number;
37+ meanV: number;
38+ contactVel: number;
39+ peakForceN: number;
40+ contactMs: number;
41+ contactEvents: number;
42+ peakXiMm: number;
43+ peakDispMm: number;
44+ peakDispThickness: number;
45+ peakAccel: number;
46+ strikeRadiusM: number;
47+ pickupRadiusM: number;
48+ RsM: number;
49+ nsteps: number;
50+ simRate: number;
51+ nDropped: number;
52+ durationS: number;
53+ renderMs: number;
54+ oversample: number;
55+ nTheta: number;
56+ warnings: string[];
57+ };
58+}
59+
60+export interface ProgressMessage {
61+ type: 'progress';
62+ id: number;
63+ stage: 'modes' | 'integrate' | 'output';
64+ frac: number;
65+}
66+
67+export interface ErrorMessage {
68+ type: 'error';
69+ id: number;
70+ message: string;
71+ file?: string;
72+ line?: number;
73+}
74+
75+interface ModesEntry {
76+ modeF: Float32Array;
77+ modeM: Float32Array;
78+}
79+
80+const modesCache = new Map<string, ModesEntry>();
81+const MODES_CACHE_MAX = 4;
82+
83+function lruGet<T>(cache: Map<string, T>, key: string): T | undefined {
84+ const v = cache.get(key);
85+ if (v !== undefined) {
86+ cache.delete(key);
87+ cache.set(key, v);
88+ }
89+ return v;
90+}
91+
92+function lruSet<T>(cache: Map<string, T>, key: string, value: T, max: number): void {
93+ cache.delete(key);
94+ cache.set(key, value);
95+ while (cache.size > max) {
96+ const oldest = cache.keys().next().value as string;
97+ cache.delete(oldest);
98+ }
99+}
100+
101+function keyFor(params: Cymbal2Params, keys: readonly (keyof Cymbal2Params)[]): string {
102+ return keys.map((k) => String(params[k])).join('|');
103+}
104+
105+// MATLAB source lines assigning every parameter onto a struct P.
106+function pLiteral(params: Cymbal2Params): string {
107+ const lines: string[] = [];
108+ for (const [k, v] of Object.entries(params)) {
109+ if (k === 'bit_depth') continue;
110+ if (typeof v === 'number') lines.push(`P.${k} = ${v};`);
111+ else if (typeof v === 'boolean') lines.push(`P.${k} = ${v ? 1 : 0};`);
112+ }
113+ return lines.join('\n');
114+}
115+
116+function runStage(mfiles: MFile[], source: string, inputs?: Record<string, RuntimeValue>) {
117+ return executeCode(
118+ source,
119+ { initialVariableValues: inputs, displayResults: false },
120+ mfiles,
121+ 'stage.m',
122+ );
123+}
124+
125+function asScalar(v: RuntimeValue): number {
126+ if (typeof v === 'number') return v;
127+ const t = v as RuntimeTensor;
128+ if (t && t.kind === 'tensor' && t.data.length === 1) return t.data[0];
129+ if (typeof v === 'boolean') return v ? 1 : 0;
130+ throw new Error('expected scalar result from numbl stage');
131+}
132+
133+function toF32(v: RuntimeValue): Float32Array {
134+ const t = v as RuntimeTensor;
135+ if (!t || t.kind !== 'tensor') throw new Error('expected tensor result from numbl stage');
136+ const out = new Float32Array(t.data.length);
137+ for (let i = 0; i < t.data.length; i++) out[i] = t.data[i];
138+ return out;
139+}
140+
141+const SKIP_WARNING = /skipping (\S+\.m) \(syntax error at line (\d+)\)/;
142+
143+self.onmessage = (ev: MessageEvent<RenderRequest>) => {
144+ const msg = ev.data;
145+ if (msg.type !== 'render') return;
146+
147+ const warnings: string[] = [];
148+ const origWarn = console.warn;
149+ console.warn = (...args: unknown[]) => {
150+ warnings.push(args.map(String).join(' '));
151+ };
152+
153+ try {
154+ render(msg);
155+ } catch (err) {
156+ const e = err as { message?: string; file?: string; line?: number };
157+ let m: ErrorMessage = {
158+ type: 'error',
159+ id: msg.id,
160+ message: e?.message ? String(e.message) : String(err),
161+ file: typeof e?.file === 'string' ? e.file : undefined,
162+ line: typeof e?.line === 'number' ? e.line : undefined,
163+ };
164+ for (const w of warnings) {
165+ const hit = SKIP_WARNING.exec(w);
166+ if (hit) {
167+ m = { type: 'error', id: msg.id, message: 'syntax error', file: hit[1], line: Number(hit[2]) };
168+ break;
169+ }
170+ }
171+ (self as unknown as Worker).postMessage(m);
172+ } finally {
173+ console.warn = origWarn;
174+ }
175+};
176+
177+function render(msg: RenderRequest): void {
178+ const { id, mfiles, quality } = msg;
179+ const params = { ...msg.params };
180+ if (quality === 'preview') {
181+ params.oversample = 1;
182+ params.tail_s = Math.min(params.tail_s, 1.5);
183+ }
184+
185+ const post = (m: ProgressMessage) => (self as unknown as Worker).postMessage(m);
186+ const codeKey = hashMFiles(mfiles);
187+
188+ // ---- linear modes, cheap, cached on geometry ----
189+ const modesKey = codeKey + '|' + keyFor(params, MODES_KEYS);
190+ let modesEntry = lruGet(modesCache, modesKey);
191+ if (!modesEntry) {
192+ post({ type: 'progress', id, stage: 'modes', frac: 0.05 });
193+ const r = runStage(mfiles, `
194+${pLiteral(params)}
195+mF = [];
196+mM = [];
197+for mw = 0:min(10, floor(P.n_theta / 2))
198+ f = cymbal2_modes(P, mw);
199+ f = f(f > 0.5 & f < 4000);
200+ nk = min(4, length(f));
201+ mF = [mF; f(1:nk)];
202+ mM = [mM; mw * ones(nk, 1)];
203+end
204+`);
205+ modesEntry = {
206+ modeF: toF32(r.variableValues.mF),
207+ modeM: toF32(r.variableValues.mM),
208+ };
209+ lruSet(modesCache, modesKey, modesEntry, MODES_CACHE_MAX);
210+ }
211+
212+ // ---- the time integration, one long call ----
213+ post({ type: 'progress', id, stage: 'integrate', frac: 0.2 });
214+ const t0 = performance.now();
215+ const r = runStage(mfiles, `
216+${pLiteral(params)}
217+S = cymbal2_system(P);
218+MB = cymbal2_mallet(P);
219+[acc, dinfo, fhist] = cymbal2_sim(P, S, MB);
220+[y, peakAccel] = cymbal2_output(P, acc);
221+% max-pool the contact force history to about 2000 points for the plot
222+nsim = length(fhist);
223+np = min(2000, nsim);
224+blk = ceil(nsim / np);
225+np = ceil(nsim / blk);
226+fpool = max(reshape([fhist, zeros(1, np * blk - nsim)], blk, np), [], 1);
227+tpool = ((0:np - 1) * blk + blk / 2) / dinfo.simRate;
228+d_strokes = dinfo.strokesPlayed;
229+d_meanV = dinfo.meanV;
230+d_cvel = dinfo.contactVel(1);
231+d_peakF = dinfo.peakForceN;
232+d_contactMs = dinfo.contactMs;
233+d_events = dinfo.contactEvents;
234+d_xi = dinfo.peakXiMm;
235+d_disp = dinfo.peakDispM;
236+d_dispH = dinfo.peakDispThickness;
237+d_nsteps = dinfo.nsteps;
238+d_simRate = dinfo.simRate;
239+d_dropped = S.nDropped;
240+d_strikeR = S.strikeRadiusM;
241+d_pickupR = S.pickupRadiusM;
242+d_Rs = S.Rs;
243+`);
244+ const renderMs = performance.now() - t0;
245+
246+ post({ type: 'progress', id, stage: 'output', frac: 0.95 });
247+ const y32 = toF32(r.variableValues.y);
248+ const peakAccel = asScalar(r.variableValues.peakAccel);
249+ const nDropped = asScalar(r.variableValues.d_dropped);
250+
251+ const warnings: string[] = [];
252+ if (nDropped > 0) {
253+ warnings.push(
254+ `${nDropped} azimuthal columns were truncated: the radial grid cannot resolve their inner boundary layer; raise n_radial`,
255+ );
256+ }
257+ if (!params.normalize) {
258+ let peak = 0;
259+ for (let i = 0; i < y32.length; i++) peak = Math.max(peak, Math.abs(y32[i]));
260+ if (peak > 1.0) {
261+ warnings.push(
262+ `clipping: peak acceleration ${peakAccel.toExponential(2)} m/s² exceeds full_scale_accel`,
263+ );
264+ }
265+ }
266+
267+ const resp: RenderResponse = {
268+ type: 'result',
269+ id,
270+ quality,
271+ y: y32,
272+ sampleRate: params.sample_rate,
273+ forceT: toF32(r.variableValues.tpool),
274+ forceN: toF32(r.variableValues.fpool),
275+ modeF: modesEntry.modeF.slice(),
276+ modeM: modesEntry.modeM.slice(),
277+ diag: {
278+ strokesPlayed: asScalar(r.variableValues.d_strokes),
279+ meanV: asScalar(r.variableValues.d_meanV),
280+ contactVel: asScalar(r.variableValues.d_cvel),
281+ peakForceN: asScalar(r.variableValues.d_peakF),
282+ contactMs: asScalar(r.variableValues.d_contactMs),
283+ contactEvents: asScalar(r.variableValues.d_events),
284+ peakXiMm: asScalar(r.variableValues.d_xi),
285+ peakDispMm: asScalar(r.variableValues.d_disp) * 1000,
286+ peakDispThickness: asScalar(r.variableValues.d_dispH),
287+ peakAccel,
288+ strikeRadiusM: asScalar(r.variableValues.d_strikeR),
289+ pickupRadiusM: asScalar(r.variableValues.d_pickupR),
290+ RsM: asScalar(r.variableValues.d_Rs),
291+ nsteps: asScalar(r.variableValues.d_nsteps),
292+ simRate: asScalar(r.variableValues.d_simRate),
293+ nDropped,
294+ durationS: y32.length / params.sample_rate,
295+ renderMs,
296+ oversample: params.oversample,
297+ nTheta: params.n_theta,
298+ warnings,
299+ },
300+ };
301+ (self as unknown as Worker).postMessage(resp, [
302+ resp.y.buffer,
303+ resp.forceT.buffer,
304+ resp.forceN.buffer,
305+ resp.modeF.buffer,
306+ resp.modeM.buffer,
307+ ]);
308+}
webapp/tsconfig.jsonadded+20−0View file
@@ -0,0 +1,20 @@
1+{
2+ "compilerOptions": {
3+ "target": "ES2022",
4+ "useDefineForClassFields": true,
5+ "lib": ["ES2022", "DOM", "DOM.Iterable", "WebWorker"],
6+ "module": "ESNext",
7+ "skipLibCheck": true,
8+ "moduleResolution": "bundler",
9+ "resolveJsonModule": true,
10+ "isolatedModules": true,
11+ "noEmit": true,
12+ "jsx": "react-jsx",
13+ "strict": true,
14+ "noUnusedLocals": true,
15+ "noUnusedParameters": true,
16+ "noFallthroughCasesInSwitch": true
17+ },
18+ "include": ["src"],
19+ "exclude": ["src/selftest.ts"]
20+}
webapp/vite.config.tsadded+9−0View file
@@ -0,0 +1,9 @@
1+import { defineConfig } from 'vite';
2+import react from '@vitejs/plugin-react';
3+
4+export default defineConfig({
5+ plugins: [react()],
6+ base: './',
7+ worker: { format: 'es' },
8+ server: { fs: { allow: ['..'] } }, // the .m sources live in ../matlab
9+});