Remove the stand
The stand (the felt pads as parameters, the lift-off rule, the viscoelastic felt) is taken out of the parameters, both solvers, the GPU kernel, the renderer, the mode browser, the tutorial and the tests, to be redesigned from scratch. Every plate is free; the kernel's inertia relief is unconditional, and the tutorial has eighteen pages.
30 changed files+167−1084
README.mdmodified+15−18View file
@@ -2,7 +2,7 @@
22
33 A browser application for designing a plate or a cymbal and simulating a drumstick strike on it. The design is a short script in MATLAB syntax (material with its loss constants, bounding box, outline, thickness field, an optional stress-free profile), run in the page by [numbl](https://numbl.org), that the user edits in the page or that an AI writes for them. The sound is computed in one of two ways. The **linear model** computes the modes of the plate once per design and synthesizes the strike from them in about a second, on the CPU, in a pool of Web Workers. The **full model** solves the geometrically nonlinear (von Karman) plate equations with in-plane inertia and a spatially varying thickness, stepped explicitly on the GPU through WebGPU; it takes tens of seconds per second of sound and captures the pitch glide and the wash of a hard hit. Both use the same grid, so the linear sound is the small-amplitude limit of the full one. The model and the numerical method are described in [THEORY.md](THEORY.md).
44
5-The focus is on free plates, sheets, discs and ellipses, flat by default, and on a realistic cymbal stand: a rod through the center hole with a felt pad of given radius below the cymbal and optionally another above it, which hold the cymbal level and let it bounce, rock, and (in the full model) lift off. The stand is a setting of the setup rather than of the design, so any plate can be put on it or taken off it.
5+The focus is on free plates, sheets, discs and ellipses, flat by default; every plate hangs freely.
66
77 ## Running
88
@@ -15,11 +15,11 @@ Requires Node 20 or later. The linear model runs in any modern browser; the full
1515
1616 ## Using the app
1717
18-The sidebar holds the run controls at its top, then the design script, the strike, listener, stand and numerics parameters, and the mode browser. The main area shows the plate (a lit 3D height field with WebGPU, or a 2D color map; drag to orbit, wheel to zoom), the movie and playback toolbar, and the spectrogram of the listener signal.
18+The sidebar holds the run controls at its top, then the design script, the strike, listener and numerics parameters, and the mode browser. The main area shows the plate (a lit 3D height field with WebGPU, or a 2D color map; drag to orbit, wheel to zoom), the movie and playback toolbar, and the spectrogram of the listener signal.
1919
2020 **Run mode.** The select at the top chooses the linear or the full model. In linear mode the modes of the current design are computed as soon as the script compiles (a progress bar shows the slices), and Synthesize produces the sound; after the first synthesis, a change of the strike or the listener parameters re-synthesizes automatically. In full mode Start begins the GPU run; Pause halts it, Stop ends it and keeps the audio and movie, and Continue extends a finished run after the duration has been raised.
2121
22-**Design script.** Pick a preset from the dropdown (a flat steel rectangle, an aluminum square, a bowed rectangle, a flat bronze disc, a tapered 16 inch crash, a 20 inch ride, the same crash with a bow and a bell, a steel ellipse, a thunder sheet) and edit it freely; the script is compiled as you type (or with Ctrl+Enter) and the view shows the resulting stress-free shape. Errors are reported under the script, and a summary lists the bounding box, the material, the thickness range, the mass, the grid, the profile's peak-to-peak height and largest slope (with a warning above the shallow-shell range), whether it sits on the stand, the loss constants with the coincidence frequency, and the estimated decay. Revert restores the selected preset.
22+**Design script.** Pick a preset from the dropdown (a flat steel rectangle, an aluminum square, a bowed rectangle, a flat bronze disc, a tapered 16 inch crash, a 20 inch ride, the same crash with a bow and a bell, a steel ellipse, a thunder sheet) and edit it freely; the script is compiled as you type (or with Ctrl+Enter) and the view shows the resulting stress-free shape. Errors are reported under the script, and a summary lists the bounding box, the material, the thickness range, the mass, the grid, the profile's peak-to-peak height and largest slope (with a warning above the shallow-shell range), the loss constants with the coincidence frequency, and the estimated decay. Revert restores the selected preset.
2323
2424 The script fills in a struct named `design`. All lengths are meters. The outline, the thickness, and the profile are numbers or function handles `@(x, y, r, theta)` of the point relative to the center of the bounding box, its distance from the center, and its angle from the x axis; the handles are applied to whole grid matrices, so they use elementwise operators. A script is self-contained: nothing is in scope but numbl's own functions, so a script defines the shapes it needs as ordinary MATLAB functions at the bottom of the file, and the unit constants at the top.
2525
@@ -45,39 +45,36 @@ v = min(max((x - a) ./ (b - a), 0), 1);
4545 end
4646 ```
4747
48-The decay of the sound is not prescribed but derived: from the loss factor of the metal's internal friction (`lossFactor`, 0 when omitted), the thermoelastic damping that follows from the three thermal constants (all three or none; without them there is no thermoelastic loss) and `design.temperature` (293 K by default), the radiation of each mode into the air, which the linear model computes from the mode's shape, and the felt of the stand. A script that still sets `design.loss` or `design.stand` is refused with a message saying so.
48+The decay of the sound is not prescribed but derived: from the loss factor of the metal's internal friction (`lossFactor`, 0 when omitted), the thermoelastic damping that follows from the three thermal constants (all three or none; without them there is no thermoelastic loss) and `design.temperature` (293 K by default), the radiation of each mode into the air, which the linear model computes from the mode's shape. A script that still sets `design.loss` is refused with a message saying so.
4949
50-**Modes.** In linear mode the mode browser lists every mode up to the cutoff frequency with its decay time and, beside it, the decay time each mechanism alone would give: the air (radiation), the metal (internal friction), thermal (thermoelastic) and the felt of a stand; click a row to draw the mode's shape on the plate (scaled to 1 mm), and filter the list by frequency. On a stand the first modes are the rock and the bounce of the cymbal on the felt. The modal frequencies are also drawn as ticks on the left edge of the spectrogram, in either run mode, so that the partials of a full run can be compared with the modes.
51-
52-**Stand.** The Stand group puts the plate on the stand or takes it off (a free plate), and sets the felt annulus (inner and outer radius, thickness, Young's modulus in compression, loss factor), whether there is a top pad under the wing nut, its preload, and the position of the rod relative to the center of the bounding box. The pads act on the plate nodes within the annulus whether or not the outline has a hole there; the Notes under the group give the total stiffness, the bounce frequency of the rigid plate, the number of pad nodes, the preloads, and the fraction of the stiffness the felt keeps at rest. The felt is viscoelastic, a chain of relaxation times that gives it the same loss factor at every frequency, so the modulus is the value at the bounce frequency and rises slowly above it; both models use the same chain (THEORY.md, Section 3). The stand is on by default, so a preset without a hole, such as the flat rectangle, sits on it until the box is unticked.
50+**Modes.** In linear mode the mode browser lists every mode up to the cutoff frequency with its decay time and, beside it, the decay time each mechanism alone would give: the air (radiation), the metal (internal friction) and thermal (thermoelastic); click a row to draw the mode's shape on the plate (scaled to 1 mm), and filter the list by frequency. The modal frequencies are also drawn as ticks on the left edge of the spectrogram, in either run mode, so that the partials of a full run can be compared with the modes.
5351
5452 **Strike and listener.** The strike point is given as fractions of the bounding box (the yellow ring in the view), with the stick speed, its effective mass, the Hertz contact stiffness and exponent, and the tip radius. The listener distance, its angle from the plate's normal, and its azimuth can be changed at any time; the microphone in the view moves with them, drawn nearer than the listener is so that both fit in the frame.
5553
56-**Numerics.** The grid spacing (8 mm by default), the cutoff frequency of the modes, the lift-off of the pads in the full model, the sample rate, and the frame rate. The Notes under the group report the grid, the frequencies at which grid dispersion lowers the partials by 3 and 10 percent, the count and time of the modes, the time step and step count of the full model, and the decay times of the fitted two-parameter loss of the full model next to the estimate it was fitted to.
54+**Numerics.** The grid spacing (8 mm by default), the cutoff frequency of the modes, the sample rate, and the frame rate. The Notes under the group report the grid, the frequencies at which grid dispersion lowers the partials by 3 and 10 percent, the count and time of the modes, the time step and step count of the full model, and the decay times of the fitted two-parameter loss of the full model next to the estimate it was fitted to.
5755
5856 **Playback.** After a run, drag the time slider to explore the movie; Play starts the audio from the slider position and moves the movie and the spectrogram cursor with it; Design returns to the stress-free shape and, in full mode, Live to the running view. The audio can be played or downloaded as a 16-bit WAV at any time; Absolute level maps 1 Pa to full scale shifted by the gain slider (the peak level in dBFS is shown with a clipping warning), and the two normalized modes scale each run to its loudness or its peak. When the signal to be played is loud enough to hurt in headphones (a peak within 10 dB of full scale, or a short-term level of -20 dBFS RMS over some 50 ms window, which every loudness-normalized hit reaches), Play first shows a warning to turn the volume down, with a Play anyway button. Once accepted, the same signal replays without a prompt; a new recording, another level mode or a change of the gain brings the warning back.
5957
60-**Tutorial.** The Tutorial button at the top of the sidebar opens a panel beside the view with a guided tour of the physics and the mathematics of the simulator in nineteen pages: the bending stiffness and the plate equation, dispersion, modes and their scaling laws, the density of modes, the grid and its dispersion error, the eigenvalue computation, decay times, the strike and the driving-point impedance, the modal force and the choice of the strike point, the modal synthesis, radiation and directivity, the spectrogram, the stand, the taper, curvature, the nonlinear model, and the limitations. Each page reads the design in the editor and shows its numbers (the bending stiffness, the wavelength at 1 kHz, the count of modes against Weyl's estimate, the gains of the first modes at the strike point and at the listener, the bounce and the rock of a stand against their formulas), and the buttons under Try it load presets, change parameters, compute the modes, strike and play, so that every formula can be checked against the model; a few pages remember mode lists for a comparison, such as the plate twice as thick against the reference. The panel keeps its page between visits, and `?tutorial=1` or `?tutorial=<page id>` opens it from the address.
58+**Tutorial.** The Tutorial button at the top of the sidebar opens a panel beside the view with a guided tour of the physics and the mathematics of the simulator in eighteen pages: the bending stiffness and the plate equation, dispersion, modes and their scaling laws, the density of modes, the grid and its dispersion error, the eigenvalue computation, decay times, the strike and the driving-point impedance, the modal force and the choice of the strike point, the modal synthesis, radiation and directivity, the spectrogram, the taper, curvature, the nonlinear model, and the limitations. Each page reads the design in the editor and shows its numbers (the bending stiffness, the wavelength at 1 kHz, the count of modes against Weyl's estimate, the gains of the first modes at the strike point and at the listener), and the buttons under Try it load presets, change parameters, compute the modes, strike and play, so that every formula can be checked against the model; a few pages remember mode lists for a comparison, such as the plate twice as thick against the reference. The panel keeps its page between visits, and `?tutorial=1` or `?tutorial=<page id>` opens it from the address.
6159
62-**URL parameters.** `?preset=crash-16` selects a preset, `mode=full` the run mode, and any strike, listener, stand, or numerics parameter can be preset in SI units (booleans as 0 or 1), for example `?preset=ride-20&spacing=0.006&duration=3&autostart=1` or `?preset=flat-rectangle&standOn=0`; `nogpu=1` disables WebGPU.
60+**URL parameters.** `?preset=crash-16` selects a preset, `mode=full` the run mode, and any strike, listener, or numerics parameter can be preset in SI units, for example `?preset=ride-20&spacing=0.006&duration=3&autostart=1`; `nogpu=1` disables WebGPU.
6361
6462 ## Tests
6563
6664 `npm test` runs eight groups of tests with Node's built-in runner; `npm run test:cpu` runs the first seven without a browser.
6765
68-- `tests/design.test.ts`: the design script (compilation, the checks with their messages, hoisted functions, the mask of a disc with a hole, the orientation of the sampled matrices, the removal of the mean plane, the stand from the parameters with its rates and preloads, the loss constants of the material), and that every preset evaluates on the default grid.
69-- `tests/reference.test.ts`: the float64 reference of the discretization (rigid modes of the stiffness, forces as the gradient of the energy, dissipativity of the losses, the work per cycle of the felt's chain, thickness scaling, the stand's stiffness against translation and tilt, the lift-off rule).
70-- `tests/assemble.test.ts`: the assembled stiffness equals the linearized reference force on a flat plate and on a bowed, tapered, holed plate on a stand.
71-- `tests/modal.test.ts`: the banded factorization and its inertia count, the tridiagonal eigensolver, all modes of a small square and a small disc against a dense solve (degenerate pairs included), the free square and disc frequencies against Gorman's and Leissa's values, and a plate on a stand.
66+- `tests/design.test.ts`: the design script (compilation, the checks with their messages, hoisted functions, the mask of a disc with a hole, the orientation of the sampled matrices, the removal of the mean plane, the loss constants of the material), and that every preset evaluates on the default grid.
67+- `tests/reference.test.ts`: the float64 reference of the discretization (rigid modes of the stiffness, forces as the gradient of the energy, dissipativity of the losses, thickness scaling).
68+- `tests/assemble.test.ts`: the assembled stiffness equals the linearized reference force on a flat plate and on a bowed, tapered, holed plate.
69+- `tests/modal.test.ts`: the banded factorization and its inertia count, the tridiagonal eigensolver, all modes of a small square and a small disc against a dense solve (degenerate pairs included), and the free square and disc frequencies against Gorman's and Leissa's values.
7270 - `tests/synth.test.ts`: a single mode under the stick follows the closed-form damped sinusoid, the radiation gain of a piston equals the analytic Rayleigh sum, the attack and the tail join, the parts computed by the workers sum to the whole.
7371 - `tests/loss.test.ts`: Zener's thermoelastic loss factor peaks where it should, the outline perimeter of a rectangle and a disc, the limits of the average radiation efficiency, the two-parameter fit, the radiation rate of a piston against the closed form and the FFT convolution against the direct double sum for every mode, and the thermoelastic rate of a uniform plate.
74-- `tests/felt.test.ts`: the felt's chain (flat loss factor across the band, storage of one at the reference, the exact step of an arm against a ramp and a sinusoid).
7572 - `tests/tutorial.test.ts`: every page of the tutorial renders with no design, with a plate, and with its modes, without undefined values and with balanced math delimiters; the actions run against an in-process host; the scaling pages' comparisons give the predicted ratios (2 for twice the thickness, 4 for half the size on a grid of the same cell count) to rounding.
76-- `tests/gpu.test.ts`: serves the app with Vite and drives a headless Chromium through Playwright: the in-page self-test of the WGSL kernel against the float64 reference on five presets (`?selftest=1` shows the same numbers in the browser), the linear model's modes and synthesis, the tutorial panel with its typeset mathematics and an action that computes the modes, the loudness warning before playback, an error in a script, short full runs on a free plate and on a stand, the absence of a high-frequency noise floor in a curved design's full run, and a light hit simulated both ways, whose spectra must agree. It needs Playwright's Chromium (`npx playwright install chromium`) and a machine whose Chromium has WebGPU.
73+- `tests/gpu.test.ts`: serves the app with Vite and drives a headless Chromium through Playwright: the in-page self-test of the WGSL kernel against the float64 reference on five presets (`?selftest=1` shows the same numbers in the browser), the linear model's modes and synthesis, the tutorial panel with its typeset mathematics and an action that computes the modes, the loudness warning before playback, an error in a script, short full runs on a flat and on a tapered plate, the absence of a high-frequency noise floor in a curved design's full run, and a light hit simulated both ways, whose spectra must agree. It needs Playwright's Chromium (`npx playwright install chromium`) and a machine whose Chromium has WebGPU.
7774
7875 ## Layout
7976
80-- `src/params.ts`: simulation parameters (strike, listener, stand, numerics), grid and mask, the plate as evaluated on the grid, material functions, the stand from the parameters and its rates, patches, listener geometry, grid dispersion.
77+- `src/params.ts`: simulation parameters (strike, listener, numerics), grid and mask, the plate as evaluated on the grid, material functions, patches, listener geometry, grid dispersion.
8178 - `src/design.ts`, `src/mruntime.ts`, `src/presets.ts`: running the script with numbl, the checks on the design struct, evaluation on the grid, and the preset scripts.
8279 - `src/loss.ts`: the losses from physical constants (internal friction, Zener's thermoelastic damping, the average radiation efficiency, the outline perimeter) and the two-parameter fit of the full model.
8380 - `src/model/reference.ts`: float64 reference of the discretization (forces, losses, energy, the time-step estimate).
@@ -86,7 +83,7 @@ The decay of the sound is not prescribed but derived: from the loss factor of th
8683 - `src/modal/damping.ts`: the per-mode radiation rate from the mode shape by FFT convolution of the radiation resistance, and the per-mode thermoelastic rate weighted by the local bending energy.
8784 - `src/modal/synth.ts`, `worker.ts`, `pool.ts`: the modal synthesis (contact, radiation gains, attack table, tail, movie), the worker, and the pool that runs the eigenvalue jobs, the per-mode losses, and the synthesis split between the workers.
8885 - `src/sim/shaders/`, `src/sim/simulation.ts`, `src/sim/selftest.ts`: the WGSL kernels (step, listener sum, extraction), the WebGPU driver, and the in-browser comparison with the reference.
89-- `src/render/`: 2D and 3D rendering of the plate with the markers of the strike, the stand, and the listener, with a 2D canvas fallback.
86+- `src/render/`: 2D and 3D rendering of the plate with the markers of the strike and the listener, with a 2D canvas fallback.
9087 - `src/audio/`: FFT, spectrogram, WAV encoder, playback, level mapping.
9188 - `src/ui/`: the script editor, the parameter form, the mode browser; `src/main.ts` wires it all together.
9289 - `src/tutorial/steps.ts`, `tutorial.ts`: the pages of the tutorial (text with TeX mathematics, live numbers from the application's state, and actions on a small host interface that `main.ts` implements), and the panel that shows them with KaTeX.
THEORY.mdmodified+16−36View file
@@ -1,8 +1,8 @@
11 # Model and numerical method
22
3-This note describes the model behind cymbalsim: a thin metal plate of arbitrary outline and spatially varying thickness, flat or with a stress-free curved profile, either hanging freely or sitting on a cymbal stand, struck once with a drumstick. The plate is described by a script; the stand, the strike and the listener are parameters of the setup. The sound is computed in one of two ways. The *linear model* takes the equations of the plate linearized about its rest shape, computes their modes once per design, and synthesizes the strike from the modes; this takes a few seconds and needs no GPU. The *full model* steps the geometrically nonlinear plate equations of Foppl and von Karman explicitly on the GPU; this takes tens of seconds per second of sound and captures what the linear model cannot, the pitch glide and the transfer of energy toward high frequencies that turn the initial tone of a hard hit into a wash. Both use the same spatial discretization, so the linear model is exactly the small-amplitude limit of the full one, and a light hit sounds the same either way.
3+This note describes the model behind cymbalsim: a thin metal plate of arbitrary outline and spatially varying thickness, flat or with a stress-free curved profile, hanging freely, struck once with a drumstick. The plate is described by a script; the strike and the listener are parameters of the setup. The sound is computed in one of two ways. The *linear model* takes the equations of the plate linearized about its rest shape, computes their modes once per design, and synthesizes the strike from the modes; this takes a few seconds and needs no GPU. The *full model* steps the geometrically nonlinear plate equations of Foppl and von Karman explicitly on the GPU; this takes tens of seconds per second of sound and captures what the linear model cannot, the pitch glide and the transfer of energy toward high frequencies that turn the initial tone of a hard hit into a wash. Both use the same spatial discretization, so the linear model is exactly the small-amplitude limit of the full one, and a light hit sounds the same either way.
44
5-Section 1 gives the plate model, Section 2 its discretization on the grid, Section 3 the stand, Section 4 the linear model and the modal synthesis, Section 5 the full model, Section 6 the checks that were made, and Section 7 the limitations. The discretization and the GPU kernel follow those of steelsim41, the earlier simulator of a scripted cymbal, simplified where the modal approach allows it; what carries over is repeated here so that this note stands alone.
5+Section 1 gives the plate model, Section 2 its discretization on the grid, Section 3 the linear model and the modal synthesis, Section 3 the full model, Section 4 the checks that were made, and Section 5 the limitations. The discretization and the GPU kernel follow those of steelsim41, the earlier simulator of a scripted cymbal, simplified where the modal approach allows it; what carries over is repeated here so that this note stands alone.
66
77 ## 1. Plate model
88
@@ -39,9 +39,9 @@ which reduces to $D\nabla^4 w$ where $D$ is constant. Bending acts on the deviat
3939
4040 **Thickness field.** A real cymbal is lathed to a taper, thick at the bell and thin at the edge. The thickness enters in three places: the mass per unit area $\rho h$, the membrane stiffnesses $C$ and $G$, proportional to $h$, and the bending stiffness $D$, proportional to $h^3$. A plate of twice the thickness has twice the bending frequencies and the same in-plane frequencies. The thickness must be positive over the plate and should vary slowly on the scale of the grid.
4141
42-**Profile.** The stress-free shape is given by the script as a function of position. Only its slopes enter the strains, so its mean is immaterial, but its tilt is not: the von Karman strains are not invariant under rotation, and a tilted stress-free shape would add a spurious membrane coupling linear in the tilt. The mass-weighted least-squares plane is therefore subtracted from every profile. The treatment is that of a shallow shell, valid while the slopes are small compared with one; the interface reports the largest slope and warns above 0.1. The bell of a real crash rises 25 mm over a radius of 50 mm, a slope of about 0.5, so the curved preset uses a flatter bell, and this is the main point of tension between the model and a real cymbal (Section 7).
42+**Profile.** The stress-free shape is given by the script as a function of position. Only its slopes enter the strains, so its mean is immaterial, but its tilt is not: the von Karman strains are not invariant under rotation, and a tilted stress-free shape would add a spurious membrane coupling linear in the tilt. The mass-weighted least-squares plane is therefore subtracted from every profile. The treatment is that of a shallow shell, valid while the slopes are small compared with one; the interface reports the largest slope and warns above 0.1. The bell of a real crash rises 25 mm over a radius of 50 mm, a slope of about 0.5, so the curved preset uses a flatter bell, and this is the main point of tension between the model and a real cymbal (Section 6).
4343
44-**Losses.** The decay of the sound is not prescribed. A partial of angular frequency $\omega$ decays as $e^{-\sigma t}$, and its decay time $T_{60} = \ln(1000)/\sigma$ is the time in which it falls by 60 dB; the rate $\sigma$ is the sum of four mechanisms, each with its own physics and its own constants in the script. *Internal friction* of the metal is a hysteretic loss: each cycle dissipates the fraction $2\pi\eta$ of the stored energy, so a mode decays at $\sigma = \eta\omega/2$, and the loss factor $\eta$ is a material constant on the same footing as $E$ and $\rho$, measured rather than derived, of order $10^{-4}$ for bronze, steel and aluminum. *Thermoelastic damping* comes from the heat that flows across the thickness between the compressed and the stretched face of a bent plate; Zener's analysis gives the loss factor $\eta_{te}(\omega, h) = \Delta\,\omega\tau/(1 + \omega^2\tau^2)$ with the relaxation strength $\Delta = [E\alpha^2 T/(\rho c_p)]\,(1+\nu)/(1-\nu)$, the biaxial form of the beam value, and the relaxation time $\tau = h^2\rho c_p/(\pi^2 k)$, the time heat takes to diffuse across the thickness, from the thermal expansion $\alpha$, the specific heat $c_p$, the conductivity $k$ and the temperature $T$ (293 K unless the script says otherwise). The loss peaks at $\omega\tau = 1$, below 100 Hz for a millimeter of metal, and is small for the presets; it is what damps thin aluminum. *Radiation* into the air is the energy the plate loses as sound, which the linear model computes exactly for each mode from its shape (Section 4) and the full model takes from an analytic average (Section 5). It is smaller than one might expect. A plate with free edges radiates poorly below its coincidence frequency, the frequency at which the bending wavelength equals the wavelength of sound (about 17 kHz for 1 mm bronze): its modes are orthogonal to the rigid translation and rotations, so they have no net volume velocity and no dipole moment, and at a free edge the bending moment and the shear vanish, which removes the edge radiation that a simply supported plate has. The efficiency rises steeply toward coincidence, so on the presets the radiation overtakes the internal friction only in the top octaves, above about 5 kHz, and it is there that it makes the high partials die first; below that the decay of a free plate is set by the metal. The *felt* of a stand is the fourth mechanism (Section 3). The linear model adds the four rates for each mode, and the mode browser shows the decay time each would give alone. The full model needs a loss operator rather than a rate per mode and uses the standard two-parameter model of plate synthesis, a rate $\sigma_0 + \sigma_1 k^2$ for a bending wave of wavenumber $k$, fitted to the estimated curve (Section 5). The loss factors and thermal constants of the presets are typical published values, not measurements of a particular cymbal.
44+**Losses.** The decay of the sound is not prescribed. A partial of angular frequency $\omega$ decays as $e^{-\sigma t}$, and its decay time $T_{60} = \ln(1000)/\sigma$ is the time in which it falls by 60 dB; the rate $\sigma$ is the sum of three mechanisms, each with its own physics and its own constants in the script. *Internal friction* of the metal is a hysteretic loss: each cycle dissipates the fraction $2\pi\eta$ of the stored energy, so a mode decays at $\sigma = \eta\omega/2$, and the loss factor $\eta$ is a material constant on the same footing as $E$ and $\rho$, measured rather than derived, of order $10^{-4}$ for bronze, steel and aluminum. *Thermoelastic damping* comes from the heat that flows across the thickness between the compressed and the stretched face of a bent plate; Zener's analysis gives the loss factor $\eta_{te}(\omega, h) = \Delta\,\omega\tau/(1 + \omega^2\tau^2)$ with the relaxation strength $\Delta = [E\alpha^2 T/(\rho c_p)]\,(1+\nu)/(1-\nu)$, the biaxial form of the beam value, and the relaxation time $\tau = h^2\rho c_p/(\pi^2 k)$, the time heat takes to diffuse across the thickness, from the thermal expansion $\alpha$, the specific heat $c_p$, the conductivity $k$ and the temperature $T$ (293 K unless the script says otherwise). The loss peaks at $\omega\tau = 1$, below 100 Hz for a millimeter of metal, and is small for the presets; it is what damps thin aluminum. *Radiation* into the air is the energy the plate loses as sound, which the linear model computes exactly for each mode from its shape (Section 3) and the full model takes from an analytic average (Section 4). It is smaller than one might expect. A plate with free edges radiates poorly below its coincidence frequency, the frequency at which the bending wavelength equals the wavelength of sound (about 17 kHz for 1 mm bronze): its modes are orthogonal to the rigid translation and rotations, so they have no net volume velocity and no dipole moment, and at a free edge the bending moment and the shear vanish, which removes the edge radiation that a simply supported plate has. The efficiency rises steeply toward coincidence, so on the presets the radiation overtakes the internal friction only in the top octaves, above about 5 kHz, and it is there that it makes the high partials die first; below that the decay of a free plate is set by the metal. The linear model adds the three rates for each mode, and the mode browser shows the decay time each would give alone. The full model needs a loss operator rather than a rate per mode and uses the standard two-parameter model of plate synthesis, a rate $\sigma_0 + \sigma_1 k^2$ for a bending wave of wavenumber $k$, fitted to the estimated curve (Section 4). The loss factors and thermal constants of the presets are typical published values, not measurements of a particular cymbal.
4545
4646 ## 2. Spatial discretization
4747
@@ -63,29 +63,9 @@ with $a_k$, $b_k$ the second differences at node $k$ (absent where undefined) an
6363
6464 **Accuracy.** Two errors limit the accuracy of the spatial scheme. The first is the treatment of the edge. For the rectangle the scheme converges at second order: the first seven frequencies of a completely free square plate are within 8 percent of Gorman's reference values at 10 cells per side and within 2 percent at 20. The staircase edge of a curved outline is only first-order accurate; the first nine frequencies of a free circular plate are within 3 percent of Leissa's values at 40 cells across the diameter, and the degenerate pairs, which the staircase splits, differ by about 1 percent. A 16 inch cymbal at 8 mm has about 50 cells across, so its low modes are placed to within a percent or two, and the splitting of the pairs is of the same order as the splitting that the hammering of a real cymbal produces. The second error is grid dispersion, which affects every partial regardless of the outline. A plane bending wave of wavenumber $k$ along a grid axis has the discrete frequency ratio $(2 - 2\cos kh)/(kh)^2$, so a partial comes out low by 3 percent at $kh \approx 0.60$ and by 10 percent at $kh \approx 1.12$. Since the bending frequency at a given $k$ is proportional to the thickness, these points move up with the thickness: for 1 mm bronze on the 8 mm grid the 3 percent point is near 1.1 kHz and the 10 percent point near 3.6 kHz. The interface reports the two frequencies at the mean thickness of the design. The high partials are therefore audibly flattened at 8 mm, in both models alike, and a finer grid is a matter of patience rather than of the method.
6565
66-## 3. The stand
66+## 3. The linear model
6767
68-A cymbal sits on a stand: a rod through the center hole, with a felt pad below the cymbal and usually another above it under a wing nut. The stand is a setting of the setup rather than of the design, given in the interface and applied to whatever plate is loaded. We model each pad as an annulus of felt of inner radius $r_1$, outer radius $r_2$ and thickness $t$, with Young's modulus $E_p$ and loss factor $\eta$. Its stiffness per unit area in compression is $E_p/t$, and its stiffness per unit area in shear is $G_p/t$ with $G_p = E_p/(2(1+\nu_p))$ and $\nu_p = 0.3$; the shear is what holds the cymbal in place sideways, and it stands in for the plastic sleeve of the rod, which is not modeled (Section 7). The stand energy is
69-
70-$$
71-U_s = \frac12 \sum_{k \in \text{pads}} \left[ k_k\, w_k^2 + g_k (u_k^2 + v_k^2) \right],
72-$$
73-
74-where the sums run over the plate nodes within the annulus, $k_k = n_p (E_p/t)\, \Delta A_k$ and $g_k = n_p (G_p/t)\, \Delta A_k$, $n_p$ is 1 or 2 for one or two pads, and $\Delta A_k$ is the node's share of the area, scaled so that the shares add up to the exact area of the annulus (a coarse grid covers the annulus with a dozen nodes). A uniform translation of the plate then feels the total stiffness $K = n_p (E_p/t) \pi (r_2^2 - r_1^2)$, and a tilt feels a restoring moment set by the second moment of the annulus, so the finite radius of the felt is what holds the cymbal level. With the defaults of the presets (a felt of 12 mm at 0.5 MPa between 7 and 20 mm, two pads) a 16 inch bronze cymbal of about 1 kg has a total stiffness of about 90 kN/m, and treated as a rigid body it would bounce on the stand at about 45 Hz and rock at about 5 Hz. The rock does come out near that value. The bounce comes out well below it, near 15 Hz for the flat disc preset and 19 Hz for the thicker ride, because a thin plate is not rigid: the pads push on its center, whose compliance adds, loosely speaking, in series with that of the felt, so the bounce mode is a rigid translation blended with the umbrella-shaped plate mode (one nodal circle, near 38 Hz for the uniform disc), and the higher plate modes that move the center take a share of the felt's damping. The felt modulus is an estimate; felt is a nonlinear material whose modulus rises with compression, and the value can be changed in the Stand parameters.
75-
76-**Preload and lift-off.** The wing nut compresses the top pad by a preload $\delta_t$ at rest (1 mm by default), and the bottom pad carries that preload plus the weight of the plate: $\delta_b = \delta_t + Mg/(k_\infty K_1)$ with $K_1$ the stiffness of one pad and $k_\infty$ the relaxed fraction of the felt's stiffness defined below. Without a top pad $\delta_t = 0$ and the bottom pad carries the weight alone. Gravity enters only through this preload; the static sag of the plate under its own weight, a fraction of a millimeter, is neglected, and $w = 0$ is the rest state. In the full model each pad pushes only while compressed: the bottom pad by $\delta_b - w$, the top pad by $\delta_t + w$, and the constant part of the two forces, which balances the weight, is subtracted, so that within the preload the pads are the linear spring above and beyond it the cymbal lifts off one pad and loads the other. A hard hit near the edge lifts the cymbal off the bottom pad and starts a rock of a centimeter or two at the edge, as on a real stand. The linear model uses the linear spring throughout; it cannot represent the lift-off, which is a nonlinear effect of the stand rather than of the plate.
77-
78-**Damping of the felt.** A loss factor $\eta$ means a stiffness with a phase lag, $K^*(\omega) = K'(\omega)(1 + i\eta)$, so that a mode of angular frequency $\omega$ whose shape puts the energy $\tfrac12 \phi^T K_s \phi$ into the pads decays at the additional rate $\sigma = \eta\, \phi^T K_s \phi / (2\omega)$. A loss factor that is the same at every frequency cannot, however, come from a spring and a dashpot: a dashpot of rate $c$ gives $\eta = c\omega/K$, which grows with frequency. An earlier version of the full model used a dashpot matched to the bounce frequency $\omega_b$, and it over-damped the modes at $\omega$ by the factor $\omega/\omega_b$, which for a light plate on the default stand, whose bounce is near 100 Hz and many of whose modes move the center, shortened the decay of the kilohertz modes from seconds to a few hundredths of a second. We therefore model the felt as a viscoelastic solid, a spring $k_\infty$ in parallel with a chain of Maxwell arms (a spring $k_i$ in series with a dashpot, relaxation time $\tau_i$), whose complex stiffness per unit of the pad stiffness is
79-
80-$$
81-\hat K(\omega) = k_\infty + \sum_i k_i \frac{i\omega\tau_i}{1 + i\omega\tau_i}.
82-$$
83-
84-With the arm frequencies $1/(2\pi\tau_i)$ spaced two per decade from 0.5 Hz to 160 kHz and the stiffnesses growing geometrically, $k_{i+1}/k_i = r^{2\delta/\pi}$ with $r = \sqrt{10}$ the frequency ratio and $\delta = \arctan\eta$ the loss angle, the chain approximates the power law $\hat K \propto (i\omega)^{2\delta/\pi}$ of a material with a constant loss angle, and its loss factor $\mathrm{Im}\,\hat K/\mathrm{Re}\,\hat K$ stays within 6 percent of $\eta$ from 5 Hz to 10 kHz for $\eta \le 0.3$ (within 10 percent for $\eta = 0.6$). The chain is normalized so that its storage stiffness $\mathrm{Re}\,\hat K$ is 1 at the bounce frequency, where the felt's stiffness matters most: the modulus $E_p$ of the parameters is the modulus at the bounce, and the felt is softer below it and stiffer above, by the factor $10^{2\delta/\pi}$ per decade, about 1.5 for $\eta = 0.3$, as the Kramers-Kronig relations require of any material with a loss. At rest it holds the fraction $k_\infty$ of that stiffness, about 0.4 for $\eta = 0.3$, which is the stiffness the static preload and the weight see. The linear model uses the storage stiffness at the bounce in $K_s$ and gives each mode the rate $\sigma = \mathrm{Im}\,\hat K(\omega)\, \phi^T K_s \phi / (2\omega)$, which is the hysteretic rate above with the felt's actual loss stiffness at the mode's frequency. The full model carries the chain in time: each arm has at every pad node and for each of $u$, $v$, $w$ one state, its extension $p_i$, which obeys $\dot p_i = \dot w - p_i/\tau_i$, and the felt's force on the node is $-k_k[k_\infty w + \sum_i k_i p_i]$, with the compression-only rule of the previous paragraph applied to the relaxed part and the arm force scaled by the number of pads in contact. The update over a step in which the node moves by $\Delta w$ is exact for a linear motion within the step, $p_i \leftarrow p_i e^{-\Delta t/\tau_i} + \Delta w\,(\tau_i/\Delta t)(1 - e^{-\Delta t/\tau_i})$, so that arms far faster than the step are handled as well as slow ones; a node that has left every pad relaxes its transverse arms at once. The two solvers now describe the same felt, with one difference that remains. The rate of the linear model is first order in the felt's coupling, and a strongly damped element at one place also reshapes the modes that load it; for a mode with a large share of its stiffness in the felt the first-order rate is an over-estimate. We checked this on the aluminum square preset (12 mm grid) by a dense eigensolution of the linearized time-domain system, plate, pads and arms together: the mode at 117 Hz, which holds about a sixth of its stiffness in the felt, decays in 0.41 s against the 0.37 s of the first-order rate, one at 359 Hz in 0.21 s against 0.15 s, and one at 745 Hz in 0.59 s against 0.36 s, while modes with a small share agree within a few percent; the rock, whose stiffness is all felt, lies at 9.6 Hz against 11.3 Hz in the linear model, because the felt is softer at the rock's frequency than at the bounce where the modulus is defined. The full model, which steps the chain, reproduces the exact values. The test of Section 6 confirms that what the stand adds to the decay of a light hit agrees between the two.
85-
86-## 4. The linear model
87-
88-**Modes.** Collecting the linear part of Section 1 on the grid of Section 2 gives the matrix equation $M \ddot q + K q = f$, with $q$ the vector of unknowns (the transverse displacement of every plate node for a flat plate, and $u$, $v$, $w$ interleaved for a plate with a profile), $M$ the diagonal lumped mass, and $K$ the symmetric stiffness, the Hessian of $U_b + U_m + U_s$ at the rest state. The tests check that $K q$ equals the linearized force of the full model to rounding. The modes are the solutions of $K\phi = \omega^2 M \phi$, normalized so that $\phi^T M \phi = 1$; a free plate has three (flat) or six (with a profile) rigid modes at $\omega = 0$, which are computed and dropped, since a hanging plate's rigid motion is taken up by the thread and does not radiate in the model of Section 4 below. A plate on a stand has none.
68+**Modes.** Collecting the linear part of Section 1 on the grid of Section 2 gives the matrix equation $M \ddot q + K q = f$, with $q$ the vector of unknowns (the transverse displacement of every plate node for a flat plate, and $u$, $v$, $w$ interleaved for a plate with a profile), $M$ the diagonal lumped mass, and $K$ the symmetric stiffness, the Hessian of $U_b + U_m + U_s$ at the rest state. The tests check that $K q$ equals the linearized force of the full model to rounding. The modes are the solutions of $K\phi = \omega^2 M \phi$, normalized so that $\phi^T M \phi = 1$; a free plate has three (flat) or six (with a profile) rigid modes at $\omega = 0$, which are computed and dropped, since a hanging plate's rigid motion is taken up by the thread and does not radiate in the model of Section 3 below. A plate on a stand has none.
8969
9070 The stiffness is banded, with a half-bandwidth of about twice the number of nodes across the shorter side of the grid (times three with a profile), and its $LDL^T$ factorization at a shift, $K - \sigma M = LDL^T$, costs little. We use it in two ways. The number of negative pivots of $D$ is the number of modes with $\omega^2 < \sigma$ (Sylvester's law of inertia), which tells exactly how many modes lie in any interval. And the factorization applies $(K - \sigma M)^{-1}$, which turns the modes near $\sigma$ into the dominant eigenvalues of a symmetric operator that the Lanczos method finds in a few dozen steps. All modes up to a cutoff frequency (16 kHz by default) are computed by *spectrum slicing*: the interval is cut into slices of about 30 modes (100 with a profile) at boundaries equally spaced in frequency, refined where a slice is too full, the inertia at each boundary fixes the count in each slice, and each slice is solved by shift-invert Lanczos with full reorthogonalization from a random start at the slice's center, until the count of converged modes in the slice equals the count from the inertias; a slice that comes up short (a degenerate pair of a disc can be missed from an unlucky start) is restarted with a new random vector orthogonal to the modes found. The slices are independent and run in a pool of Web Workers. A 16 inch cymbal at 8 mm has about 2150 nodes and about 1900 bending modes below 16 kHz; they take about 2.5 s on a laptop with 12 threads. With a profile the same cymbal has 6500 unknowns and a band three times as wide, and the factorizations and solves become bound by memory traffic; the modes then take about 30 s on the same laptop.
9171
@@ -95,7 +75,7 @@ $$
9575 \ddot q_m + 2\sigma_m \dot q_m + \omega_m^2 q_m = \phi_m^T f(t),
9676 $$
9777
98-with $\sigma_m$ the sum of the four loss rates of the mode. The internal friction gives $\eta U_m/\omega_m$ with $U_m$ the strain energy of the plate in the mode, $\omega_m^2/2$ for a mass-normalized mode less the energy in the pads; the pads give the term of Section 3; the thermoelastic loss gives $\sum_k \eta_{te}(\omega_m, h_k)\,U_{b,k}/\omega_m$, the local loss factor weighted by the bending energy $U_{b,k}$ of the mode at each node (the same discrete bending energy as Section 2, with the cell terms shared among their corners), which reduces to $\eta_{te}\omega_m/2$ on a flat uniform plate. The radiation rate is the power the mode radiates over twice its energy. A baffled plate whose nodes move with the velocities $v_k$ radiates the power $v^T R v$ from one face, with the radiation resistance matrix $R_{jk} = (\rho_{\mathrm{air}}\omega^2/4\pi c)\,\Delta A_j \Delta A_k\, \mathrm{sinc}(k_a r_{jk})$ of Elliott and Johnson, $k_a = \omega/c$; both faces radiate, the velocity of a mass-normalized mode of amplitude $A$ is $\omega A\phi$ and its energy $\omega^2 A^2/2$, so</p>
78+with $\sigma_m$ the sum of the three loss rates of the mode. The internal friction gives $\eta U_m/\omega_m$ with $U_m$ the strain energy of the mode, $\omega_m^2/2$ for a mass-normalized mode, that is $\eta\omega_m/2$; the thermoelastic loss gives $\sum_k \eta_{te}(\omega_m, h_k)\,U_{b,k}/\omega_m$, the local loss factor weighted by the bending energy $U_{b,k}$ of the mode at each node (the same discrete bending energy as Section 2, with the cell terms shared among their corners), which reduces to $\eta_{te}\omega_m/2$ on a flat uniform plate. The radiation rate is the power the mode radiates over twice its energy. A baffled plate whose nodes move with the velocities $v_k$ radiates the power $v^T R v$ from one face, with the radiation resistance matrix $R_{jk} = (\rho_{\mathrm{air}}\omega^2/4\pi c)\,\Delta A_j \Delta A_k\, \mathrm{sinc}(k_a r_{jk})$ of Elliott and Johnson, $k_a = \omega/c$; both faces radiate, the velocity of a mass-normalized mode of amplitude $A$ is $\omega A\phi$ and its energy $\omega^2 A^2/2$, so</p>
9979
10080 $$
10181 \sigma_{\mathrm{rad},m} = \frac{\rho_{\mathrm{air}}\,\omega_m^2}{2\pi c} \sum_{j,k} \Delta A_j\,\Delta A_k\,\phi_m(j)\,\phi_m(k)\,\frac{\sin(k_a r_{jk})}{k_a r_{jk}},
@@ -122,28 +102,28 @@ which is exact for the free decay and includes the directivity of every mode. Th
122102
123103 **Movie.** The transverse displacement at the frame times is the modal sum $\sum_m \phi_m(k) q_m(t)$, with $q_m$ from the contact tables before $t_1$ and from the free decay after it. Its cost, 1900 modes times 2150 nodes times 480 frames, is the largest part of the synthesis and is split among the workers, after the audio, so that playback can start first.
124104
125-## 5. The full model
105+## 4. The full model
126106
127107 The full model keeps the nonlinear strains of Section 1 and steps the semi-discrete system with the symplectic Euler scheme, $v^{n+1} = v^n + \Delta t\, a(x^n, v^n)$ and $x^{n+1} = x^n + \Delta t\, v^{n+1}$, in single precision on the GPU. With the velocities interpreted at half steps this is the leapfrog scheme, second order and symplectic for the conservative part; the losses are evaluated explicitly at the old velocity, and the felt's arms at their state from the previous step, which is then advanced with the new displacement. The stability limit is $\Delta t < 2/\omega_{\max}$, where $\omega_{\max}$ is the largest frequency of the linearized discrete system, estimated at start-up by power iteration on a small uniform grid with the same spacing and the largest thickness of the design; 90 percent of the limit is used. The in-plane waves, with speed $c_p = \sqrt{E/\rho(1-\nu^2)}$, about 3800 m/s for bronze, set this limit, and $c_p$ does not depend on the thickness: at 8 mm spacing $\Delta t \approx 1.9\ \mu$s for bronze, and 4 s of sound take about 2 million steps. The step is rounded down so that an integer number of steps makes one audio sample. The in-plane inertia terms are usually dropped in favor of an Airy stress function, since in-plane waves are much faster than bending waves; we keep them because they make the whole system explicit and local, which suits the GPU, at the cost of the time step just described.
128108
129-**Losses.** The two-parameter model gives a bending mode of angular frequency $\omega$ and wavenumber $k$ the rate $\sigma_0 + \sigma_1 k^2$; for the flat plate $k^2 = \omega/\kappa$ with $\kappa = \sqrt{D/\rho h}$, so $\sigma_0$ sets the low-frequency decay and $\sigma_1$ makes the high frequencies decay faster in proportion to frequency. Note that the internal friction, $\eta\omega/2$, is exactly of the second form with $\sigma_1 = \eta\kappa/2$. The two parameters are fitted by non-negative least squares to an estimate of the total rate on log-spaced frequencies from 50 Hz to 10 kHz, with $\kappa$ at the mean thickness. The estimate has no modes to work with, so it takes the internal friction as is, the thermoelastic loss factor averaged over the plate area with the local thickness, and the radiation from the average radiation efficiency of a baffled plate of the design's area and perimeter (the perimeter measured along the marching-squares contour of the cell mask, so that a disc is not counted as a staircase). The average of Leppington, Broadbent and Heron for a simply supported plate grows about linearly with frequency below coincidence, peaks there at a value set by the plate's size in coincidence wavelengths, and tends to one above; the modal computation of Section 4 reproduces Wallace's classical efficiencies of simply supported mode shapes to 2 percent, and for the free-edged modes of the presets it gives values below that average by the factor $(f/f_c)^2$, within a factor of two from 100 Hz to 8 kHz, which is what the boundary terms of the far-field integral predict for edges at which the moment and the shear vanish. The estimate applies that factor. The interface reports the fitted decay times at 100 Hz, 1 kHz and 5 kHz next to the estimate, and the modal computation, when it has run, shows how far the exact per-mode rates scatter about it. The $\sigma_1$ term is implemented as $2\sigma_1$ times the graph Laplacian of the transverse velocity, whose eigenvalue on this grid is exactly $\omega/\kappa$ for a grid mode, so the intended rate holds for the discrete modes too; it does not depend on the thickness, and where the plate is thicker a given frequency corresponds to a smaller $k^2$, so the frequency-dependent loss is weaker there. The in-plane motion neither radiates nor has a thermoelastic loss in this model and carries only the internal friction, as a frequency-independent rate $\eta\omega_p/2$ at the lowest in-plane resonance $\omega_p = \pi c_p/L$ of the plate's longer side, since the stretching modes lie at several kilohertz and couple linearly to the transverse motion once the plate is curved.
109+**Losses.** The two-parameter model gives a bending mode of angular frequency $\omega$ and wavenumber $k$ the rate $\sigma_0 + \sigma_1 k^2$; for the flat plate $k^2 = \omega/\kappa$ with $\kappa = \sqrt{D/\rho h}$, so $\sigma_0$ sets the low-frequency decay and $\sigma_1$ makes the high frequencies decay faster in proportion to frequency. Note that the internal friction, $\eta\omega/2$, is exactly of the second form with $\sigma_1 = \eta\kappa/2$. The two parameters are fitted by non-negative least squares to an estimate of the total rate on log-spaced frequencies from 50 Hz to 10 kHz, with $\kappa$ at the mean thickness. The estimate has no modes to work with, so it takes the internal friction as is, the thermoelastic loss factor averaged over the plate area with the local thickness, and the radiation from the average radiation efficiency of a baffled plate of the design's area and perimeter (the perimeter measured along the marching-squares contour of the cell mask, so that a disc is not counted as a staircase). The average of Leppington, Broadbent and Heron for a simply supported plate grows about linearly with frequency below coincidence, peaks there at a value set by the plate's size in coincidence wavelengths, and tends to one above; the modal computation of Section 3 reproduces Wallace's classical efficiencies of simply supported mode shapes to 2 percent, and for the free-edged modes of the presets it gives values below that average by the factor $(f/f_c)^2$, within a factor of two from 100 Hz to 8 kHz, which is what the boundary terms of the far-field integral predict for edges at which the moment and the shear vanish. The estimate applies that factor. The interface reports the fitted decay times at 100 Hz, 1 kHz and 5 kHz next to the estimate, and the modal computation, when it has run, shows how far the exact per-mode rates scatter about it. The $\sigma_1$ term is implemented as $2\sigma_1$ times the graph Laplacian of the transverse velocity, whose eigenvalue on this grid is exactly $\omega/\kappa$ for a grid mode, so the intended rate holds for the discrete modes too; it does not depend on the thickness, and where the plate is thicker a given frequency corresponds to a smaller $k^2$, so the frequency-dependent loss is weaker there. The in-plane motion neither radiates nor has a thermoelastic loss in this model and carries only the internal friction, as a frequency-independent rate $\eta\omega_p/2$ at the lowest in-plane resonance $\omega_p = \pi c_p/L$ of the plate's longer side, since the stretching modes lie at several kilohertz and couple linearly to the transverse motion once the plate is curved.
130110
131-**Free plate.** A sheet hanging from a thread is treated as free, and the reaction to the stick is distributed over the plate as the inertial load of a rigid body, which is the inertia-relief equilibrium of a free body: the force of the previous step is applied as the rigid acceleration field $a_k = F\,(1/M + y_k y_c/I_x + x_k x_c/I_y)$ that balances the patch load and its moments about the center, so the plate gains no net momentum or angular momentum. The one-step lag has zero net impulse. A plate on a stand takes the reaction through the pads. This is the counterpart of dropping the rigid modes in the linear model.
111+**Free plate.** A sheet hanging from a thread is treated as free, and the reaction to the stick is distributed over the plate as the inertial load of a rigid body, which is the inertia-relief equilibrium of a free body: the force of the previous step is applied as the rigid acceleration field $a_k = F\,(1/M + y_k y_c/I_x + x_k x_c/I_y)$ that balances the patch load and its moments about the center, so the plate gains no net momentum or angular momentum. The one-step lag has zero net impulse. This is the counterpart of dropping the rigid modes in the linear model.
132112
133113 **Single precision.** Three places in the kernel would otherwise turn rounding into audible noise. The first is the membrane strain of a curved design. Written directly, $\varepsilon_{xx} = u_x + \tfrac12[(w_0+w)_x^2 - w_{0,x}^2]$ subtracts two nearly equal numbers once $w$ is small compared with the profile, and the rounding error of the two squares is proportional to $w_{0,x}^2$, a constant: it does not fall as the vibration decays, and it acts as a broadband forcing that leaves a noise floor which never dies away. The strains are therefore written so that every term carries the deviation, $\varepsilon_{xx} = u_x + w_x(w_{0,x} + \tfrac12 w_x)$ and $\gamma_{xy} = u_y + v_x + w_{0,x}w_y + w_{0,y}w_x + w_x w_y$ with $w$ now the deviation alone, which is algebraically the same and has nothing large to cancel. On the 16 inch crash with a bell, whose largest profile slope is about 0.23, the direct form left the top octave at a constant level for the whole run, about 25 dB below the peak of a light hit; the rearranged form puts it more than 60 dB down and decaying. The second is the update of $w$. With $w$ of a few millimeters and $\Delta t \approx 2\ \mu$s, a velocity below about $10^{-3}$ m/s produces an increment smaller than the resolution of $w$, and the rounding would again act as a broadband noise source. The update of $w$ is therefore compensated: the rounding error of each addition is carried in a separate array and added back at the next step. The third is the listener sum, which forms the acceleration as the difference of successive delayed velocities at the audio rate, the average acceleration over one sample interval, rather than sampling the instantaneous acceleration, whose per-step rounding error of order $1\ \mathrm{m/s^2}$ per node would go straight into the audio. What remains after all three is a noise floor proportional to the signal, roughly 75 dB below the peak in the top octave for a strongly curved design and 100 dB for a flat one.
134114
135115 **The kernel.** The state is stored twice and each step reads one copy and writes the other. A step is one compute dispatch in which each workgroup of 256 threads advances a 12 by 12 tile of nodes: the threads load the tile with a two-node halo into shared memory, compute the membrane forces of the cells and the bending moments at the nodes once each, and then evaluate the update of the tile's nodes from shared memory. On grids of a few thousand nodes the cost of a step is the latency of one such workgroup rather than the arithmetic. A per-step uniform, addressed by a dynamic offset, carries the losses and the audio bookkeeping, so that a whole movie frame of steps is encoded in one pass; the retarded-time listener sum runs once per frame for all of the frame's samples from a ring buffer of the transverse velocities, and the following frame is encoded while the GPU works on the current one. A 16 inch cymbal at 8 mm runs at about 10 s per second of sound on the integrated GPU of a laptop.
136116
137-## 6. Validation
117+## 5. Validation
138118
139-The tests check, without a GPU: that the linear stiffness of every configuration (uniform and tapered rectangles, a disc with a hole) has exactly the rigid modes in its null space; that the forces of the full model are the exact gradient of the discrete energy at a strongly nonlinear state, flat and bowed, uniform and tapered, free and on a stand; that the loss terms only remove energy; that the felt's chain has the loss factor it was built for within 6 percent across the band and a storage stiffness of 1 at the reference, that the exact step of an arm reproduces its response to a ramp for steps from a fiftieth to twice the relaxation time, and that the chain stepped in time on a plate on a stand dissipates the work $\pi\,\mathrm{Im}\,\hat K\, K A^2$ per cycle of a rigid translation of amplitude $A$; that the assembled $K$ equals the linearized force of the full model to $10^{-6}$ on a bowed, tapered, holed plate on a stand; that the banded factorization solves and that its inertia counts the eigenvalues; that every mode of a small square and of a small disc, degenerate pairs included, matches a dense eigensolution to $10^{-7}$ with no mode missed; that the free square and disc frequencies converge to Gorman's and Leissa's values as stated in Section 2; that the stand resists a translation with its total stiffness and a tilt with the pads' second moment, and that the lift-off rule unloads the right pad; that Zener's loss factor peaks at $\omega\tau = 1$ with the value $\Delta/2$, that the outline perimeter of a disc is within a few percent of $2\pi R$, that the average radiation efficiency has its textbook limits, that the radiation rate of a uniform piston at low frequency equals the closed form $\rho_{\mathrm{air}}\omega^2 S^2/(2\pi c M)$ and that the FFT convolution equals the direct double sum to $10^{-9}$ for every mode, and that the thermoelastic rate of a flat uniform plate is $\eta_{te}\omega/2$; that a single mode struck by the stick follows the closed-form damped sinusoid, that the radiation gain of a uniform piston equals the analytic Rayleigh sum, that the attack and the tail of the synthesized signal join, and that the parts computed by the workers sum to the whole.
119+The tests check, without a GPU: that the linear stiffness of every configuration (uniform and tapered rectangles, a disc with a hole) has exactly the rigid modes in its null space; that the forces of the full model are the exact gradient of the discrete energy at a strongly nonlinear state, flat and bowed, uniform and tapered; that the loss terms only remove energy; that the assembled $K$ equals the linearized force of the full model to $10^{-6}$ on a bowed, tapered, holed plate; that the banded factorization solves and that its inertia counts the eigenvalues; that every mode of a small square and of a small disc, degenerate pairs included, matches a dense eigensolution to $10^{-7}$ with no mode missed; that the free square and disc frequencies converge to Gorman's and Leissa's values as stated in Section 2; that Zener's loss factor peaks at $\omega\tau = 1$ with the value $\Delta/2$, that the outline perimeter of a disc is within a few percent of $2\pi R$, that the average radiation efficiency has its textbook limits, that the radiation rate of a uniform piston at low frequency equals the closed form $\rho_{\mathrm{air}}\omega^2 S^2/(2\pi c M)$ and that the FFT convolution equals the direct double sum to $10^{-9}$ for every mode, and that the thermoelastic rate of a flat uniform plate is $\eta_{te}\omega/2$; that a single mode struck by the stick follows the closed-form damped sinusoid, that the radiation gain of a uniform piston equals the analytic Rayleigh sum, that the attack and the tail of the synthesized signal join, and that the parts computed by the workers sum to the whole.
140120
141-With a GPU, a Playwright test runs one step of the WGSL kernel from a smooth random state of about 2.5 thicknesses amplitude, with the losses and the stand with lift-off and the felt's arms at random extensions, on the flat, bowed, tapered, curved and elliptical presets, and compares the accelerations with the float64 reference: they agree to between $1.5\times10^{-7}$ and $6\times10^{-7}$ relative, which is single-precision rounding, and the updated arm extensions agree to better than $10^{-6}$. A short strike runs to completion on a free plate and on a stand. A light hit on the curved crash is run to the end and the top octave of its tail is required to lie more than 50 dB below the peak, which catches a cancellation in the strains of the kind described in Section 5. And a light hit (0.05 m/s) on the flat disc is simulated both ways: over the interval from 50 to 250 ms the magnitude spectra of the two listener signals below 6 kHz have a correlation of 0.985 and their levels agree within 13 percent. The remaining difference comes from the fitted two-parameter loss against the per-mode decay times and the modal truncation at 16 kHz. Finally, a light hit on the aluminum square, the case that exposed the dashpot, is simulated both ways with the stand on and off, and what the stand adds to the decay of each band between 0.25 and 0.55 s (the drop with the stand on minus the drop without it, which cancels the fitted plate loss against the per-mode one) must agree between the models within 4 dB in the four bands from 50 Hz to 6 kHz; the differences are 1.4, 0.4, 0.7 and 2.8 dB. With the dashpot the same measure differed by tens of decibels above 600 Hz.
121+With a GPU, a Playwright test runs one step of the WGSL kernel from a smooth random state of about 2.5 thicknesses amplitude, with the losses, on the flat, bowed, tapered, curved and elliptical presets, and compares the accelerations with the float64 reference: they agree to between $1.5\times10^{-7}$ and $6\times10^{-7}$ relative, which is single-precision rounding. A short strike runs to completion on a flat and on a tapered plate. A light hit on the curved crash is run to the end and the top octave of its tail is required to lie more than 50 dB below the peak, which catches a cancellation in the strains of the kind described in Section 4. And a light hit (0.05 m/s) on the flat disc is simulated both ways: over the interval from 50 to 250 ms the magnitude spectra of the two listener signals below 6 kHz have a correlation of 0.985 and their levels agree within 13 percent. The remaining difference comes from the fitted two-parameter loss against the per-mode decay times and the modal truncation at 16 kHz.
142122
143-## 7. Limitations
123+## 6. Limitations
144124
145-The linear model is linear: it has no pitch glide, no wash, no lift-off, and a hard hit is only a loud one. It is the right tool for the tone, the modes and the decay of a design, and for a first listen, and the full model is the right tool for the sound of a hard hit; the two agree for a light one. The modal truncation at the cutoff frequency removes the highest partials, which the grid places poorly anyway; the linear interpolation of the acceleration table during the first two milliseconds slightly blurs the attack.
125+The linear model is linear: it has no pitch glide, no wash, and a hard hit is only a loud one. It is the right tool for the tone, the modes and the decay of a design, and for a first listen, and the full model is the right tool for the sound of a hard hit; the two agree for a light one. The modal truncation at the cutoff frequency removes the highest partials, which the grid places poorly anyway; the linear interpolation of the acceleration table during the first two milliseconds slightly blurs the attack.
146126
147127 The von Karman model assumes moderate rotations, so it is not valid for displacements that are large compared with the plate size; a strike that moves a 400 mm cymbal by a few millimeters is well within its range. The Marguerre treatment of the profile is a shallow-shell approximation, adequate for a bow of a few millimeters on a sheet and for a gentle cymbal bow, but strained by the bell of a real cymbal, whose slope is several times the range in which the approximation holds; the curved preset uses a flatter bell than a real cymbal, and the interface reports the largest slope of a design. A deep bell also behaves as a stiff cup that shifts the bell modes up in ways the shallow shell does not capture. The thickness field is treated as a slowly varying coefficient in a thin-plate model, which is appropriate for a lathed taper but not for a step or a groove on the scale of the grid. The hammering of a real cymbal, which leaves residual stresses and a slightly non-axisymmetric shape, is absent unless the script adds it. The losses are derived from physical constants, but two of those constants are themselves empirical: the loss factor of the metal and, to a lesser degree, its thermal constants are typical published values rather than measurements of a particular alloy after hammering, and the loss factor is the one knob a user has on the decay. The radiation loss is that of a plate in a rigid wall, consistent with the listener model but an overstatement at low frequency for a plate in free air, where the two faces partly cancel; the losses to the air's viscosity and to the support other than the felt are absent; and the full model carries the losses only through the two-parameter fit to an average, so its decay agrees with the modal one on average and not mode by mode. Grid dispersion lowers the high partials by the amounts given in Section 2, and the staircase edge places the low frequencies within a few percent and splits modes that would be degenerate on a true disc.
148128
149-The stand is a linear viscoelastic felt with a loss factor that is the same at every frequency; real felt stiffens with compression, its loss factor is not exactly flat, and the value of its modulus is an estimate. The rod and its sleeve, which stop the cymbal from sliding sideways once the felt's shear gives way, are not modeled, and neither is the friction of the felt. Gravity enters only through the preload. Air loading, which lowers the frequencies of a thin plate slightly and adds radiation damping, is not included. The plate is unbaffled in reality, so the sound from its two faces partly cancels at low frequencies, and the Rayleigh integral overestimates the low-frequency level; an accurate treatment would require solving the exterior acoustic problem, which is outside the scope of this simulator. Finally, single precision limits the dynamic range of the tail of the full model's sound to roughly 100 dB below the peak, and the linear model, in double precision, has no such limit but no more physics.
129+The plate is free: nothing supports it, and neither the thread a sheet hangs from nor the stand a cymbal sits on is modeled, so the losses and the low modes a support adds are absent. Air loading, which lowers the frequencies of a thin plate slightly and adds radiation damping, is not included. The plate is unbaffled in reality, so the sound from its two faces partly cancels at low frequencies, and the Rayleigh integral overestimates the low-frequency level; an accurate treatment would require solving the exterior acoustic problem, which is outside the scope of this simulator. Finally, single precision limits the dynamic range of the tail of the full model's sound to roughly 100 dB below the peak, and the linear model, in double precision, has no such limit but no more physics.
index.htmlmodified+1−1View file
@@ -21,7 +21,7 @@
2121 <summary>About</summary>
2222 <p>
2323 A plate or cymbal is described by a short MATLAB-syntax script (material, outline, thickness, an optional
24- profile, the stand, and the decay times), and a drumstick strike on it is simulated in one of two ways.
24+ profile, and the loss constants), and a drumstick strike on it is simulated in one of two ways.
2525 The <b>linear</b> model computes the modes of the plate once per design and synthesizes the sound from them
2626 in about a second, on the CPU. The <b>full</b> model steps the geometrically nonlinear (von Karman) plate
2727 equations on the GPU through WebGPU, which takes about a minute per second of sound but captures the pitch
package.jsonmodified+1−1View file
@@ -9,7 +9,7 @@
99 "build": "tsc --noEmit && vite build",
1010 "preview": "vite preview",
1111 "test": "node --test tests/*.test.ts",
12- "test:cpu": "node --test tests/design.test.ts tests/reference.test.ts tests/assemble.test.ts tests/modal.test.ts tests/synth.test.ts tests/loss.test.ts tests/felt.test.ts tests/tutorial.test.ts",
12+ "test:cpu": "node --test tests/design.test.ts tests/reference.test.ts tests/assemble.test.ts tests/modal.test.ts tests/synth.test.ts tests/loss.test.ts tests/tutorial.test.ts",
1313 "test:gpu": "node --test tests/gpu.test.ts"
1414 },
1515 "devDependencies": {
src/design.tsmodified+5−10View file
@@ -1,8 +1,7 @@
11 // The design: a MATLAB-syntax script, run by numbl, that fills in a struct
22 // named `design` with the material, the bounding box, the outline, the
33 // thickness field, the stress-free profile, and the loss constants of the
4-// material (from which the decay is derived, loss.ts). The stand is not part
5-// of the design: it is a parameter of the setup (params.ts).
4+// material (from which the decay is derived, loss.ts).
65 // The outline, the thickness, and the profile are function handles
76 // @(x, y, r, theta), vectorized over grid matrices (or plain numbers).
87 // compileDesign runs the script and checks the struct; evaluatePlate runs it
@@ -20,7 +19,6 @@ import {
2019 type Grid,
2120 type Material,
2221 type Plate,
23- type Stand,
2422 type Thermal,
2523 ROOM_TEMPERATURE,
2624 insideNode,
@@ -150,7 +148,7 @@ export function compileDesign(source: string): Design {
150148 };
151149 if (material.nu >= 0.5) throw new DesignError("design.material.nu must be below 0.5");
152150 if (material.eta >= 0.1) throw new DesignError("design.material.lossFactor must be small (metals have loss factors of 1e-4 to 1e-3)");
153- if (field(d, "loss") !== undefined) throw new DesignError("design.loss is no longer used: the decay of each mode follows from design.material.lossFactor, the thermal constants alpha, cp and k, the radiation of the plate, and the felt of the stand");
151+ if (field(d, "loss") !== undefined) throw new DesignError("design.loss is no longer used: the decay of each mode follows from design.material.lossFactor, the thermal constants alpha, cp and k, and the radiation of the plate");
154152 const temperature = has(d, "temperature") ? num(field(d, "temperature"), "design.temperature (K)", 1) : ROOM_TEMPERATURE;
155153
156154 const b = field(d, "bounds");
@@ -183,8 +181,6 @@ export function compileDesign(source: string): Design {
183181 }
184182 }
185183
186- if (field(d, "stand") !== undefined) throw new DesignError("design.stand is no longer part of the script: the stand (the felt pads, the preload, the position of the rod) is set in the Stand parameters beside the script and applies to any design");
187-
188184 return { name, material, Lx, Ly, hasInside, thicknessConst, profileConst, temperature, source };
189185 }
190186
@@ -245,8 +241,7 @@ function matrix(vars: Record<string, RuntimeValue>, name: string, rows: number,
245241 }
246242
247243 /**
248- * Samples the design on a grid of about the requested spacing, with the given
249- * stand (from the parameters, not the script): the mask from
244+ * Samples the design on a grid of about the requested spacing: the mask from
250245 * the outline at the cell centers, the thickness at the nodes (which must be
251246 * positive on the plate), the profile at the nodes with its mean plane
252247 * removed. Only the slopes of
@@ -254,7 +249,7 @@ function matrix(vars: Record<string, RuntimeValue>, name: string, rows: number,
254249 * not: the von Karman strains are not invariant under rotation, and a tilted
255250 * stress-free shape would add a spurious membrane coupling.
256251 */
257-export function evaluatePlate(d: Design, spacing: number, stand: Stand | null = null): Plate {
252+export function evaluatePlate(d: Design, spacing: number): Plate {
258253 const full = makeGrid(d.Lx, d.Ly, spacing);
259254 const { nx, ny } = full;
260255 const vars = runDesignScript(d.source + trailer(d, full), d.source.split("\n").length);
@@ -285,7 +280,7 @@ export function evaluatePlate(d: Design, spacing: number, stand: Stand | null =
285280 if (!any) throw new DesignError("the outline contains no grid cells; enlarge it or reduce the grid spacing");
286281 removePlane(w0, g);
287282 for (let j = 0; j < ny; j++) for (let i = 0; i < nx; i++) if (!insideNode(g, i, j)) w0[j * nx + i] = 0;
288- return { name: d.name, Lx: d.Lx, Ly: d.Ly, grid: g, mat: d.material, h, w0, stand, temperature: d.temperature };
283+ return { name: d.name, Lx: d.Lx, Ly: d.Ly, grid: g, mat: d.material, h, w0, temperature: d.temperature };
289284 }
290285
291286 /** Whether the plate is flat (no profile), in which case bending and in-plane motion decouple in the linear model. */
src/felt.tsdeleted+0−85View file
@@ -1,85 +0,0 @@
1-// The felt of the stand as a viscoelastic material: a spring in parallel
2-// with a chain of Maxwell arms (a spring and a dashpot in series), whose
3-// relaxation times are spaced two per decade from 0.5 Hz to above the audio
4-// band and whose weights grow geometrically. Such a chain has a loss factor
5-// that is nearly flat over the band, which is what a single loss factor
6-// eta means, and unlike a dashpot it is causal and can be stepped in time:
7-// each arm carries one state per node, its extension p, which relaxes
8-// toward zero at its own rate while the node moves.
9-//
10-// Every quantity here is normalized to the stiffness at the reference
11-// frequency (the bounce of the plate on the stand), where the storage
12-// stiffness of the chain is 1, so that the felt modulus given by the user is
13-// the modulus at the frequency where it matters most. Below the reference
14-// the felt is softer and above it stiffer, as the Kramers-Kronig relations
15-// require of any material with a loss; at rest its stiffness is kInf.
16-
17-export interface FeltChain {
18- /** Relaxation time of each arm [s]. */
19- tau: Float64Array;
20- /** Stiffness of each arm, per unit of the stiffness at the reference frequency. */
21- k: Float64Array;
22- /** Relaxed (static) stiffness, per unit of the stiffness at the reference frequency. */
23- kInf: number;
24- /** The loss factor the chain was built for. */
25- eta: number;
26-}
27-
28-/** Lowest arm frequency [Hz]; the highest is FELT_ARMS_PER_DECADE arms per decade above it. */
29-export const FELT_LOWEST_HZ = 0.5;
30-export const FELT_ARMS_PER_DECADE = 2;
31-/** Number of arms: 0.5 Hz to 158 kHz. */
32-export const FELT_ARMS = 12;
33-
34-/**
35- * The chain for a loss factor eta, normalized at the frequency fRef [Hz]. A
36- * chain of arms with frequencies in the ratio r and stiffnesses in the ratio
37- * 1 + c approximates the power law K(i omega) ~ (i omega)^(2 delta / pi) of a
38- * material with the constant loss angle delta = atan(eta), whose modulus
39- * rises by the factor r^(2 delta / pi) per arm; c is that factor minus one.
40- */
41-export function feltChain(eta: number, fRef: number): FeltChain {
42- const r = Math.pow(10, 1 / FELT_ARMS_PER_DECADE);
43- const c = eta > 0 ? Math.pow(r, (2 * Math.atan(eta)) / Math.PI) - 1 : 0;
44- const nArms = c > 0 ? FELT_ARMS : 0;
45- const tau = new Float64Array(nArms), k = new Float64Array(nArms);
46- let S = 1;
47- for (let i = 0; i < nArms; i++) {
48- tau[i] = 1 / (2 * Math.PI * FELT_LOWEST_HZ * Math.pow(r, i));
49- k[i] = c * S;
50- S *= 1 + c;
51- }
52- const raw: FeltChain = { tau, k, kInf: 1, eta };
53- const scale = 1 / feltModulus(raw, 2 * Math.PI * Math.max(fRef, 1e-6)).storage;
54- for (let i = 0; i < nArms; i++) k[i] *= scale;
55- return { tau, k, kInf: scale, eta };
56-}
57-
58-/** Storage and loss stiffness of the chain at the angular frequency omega, per unit of the reference stiffness. */
59-export function feltModulus(ch: FeltChain, omega: number): { storage: number; loss: number } {
60- let storage = ch.kInf, loss = 0;
61- for (let i = 0; i < ch.tau.length; i++) {
62- const x = omega * ch.tau[i], d = 1 + x * x;
63- storage += (ch.k[i] * x * x) / d;
64- loss += (ch.k[i] * x) / d;
65- }
66- return { storage, loss };
67-}
68-
69-/**
70- * The constants of the exact time step of the arms for a step dt: with the
71- * node moving linearly by dw over the step, the extension of arm i goes from
72- * p to p decay_i + dw growth_i, where decay_i = exp(-dt / tau_i) and
73- * growth_i = (tau_i / dt)(1 - decay_i). The force of the arms on the node is
74- * -k_node sum_i k_i p_i.
75- */
76-export function feltStepConstants(ch: FeltChain, dt: number): { decay: Float64Array; growth: Float64Array } {
77- const n = ch.tau.length;
78- const decay = new Float64Array(n), growth = new Float64Array(n);
79- for (let i = 0; i < n; i++) {
80- const x = dt / ch.tau[i];
81- decay[i] = Math.exp(-x);
82- growth[i] = -Math.expm1(-x) / x;
83- }
84- return { decay, growth };
85-}
src/loss.tsmodified+0−2View file
@@ -10,8 +10,6 @@
1010 // the plate loses as sound; the linear model computes it exactly for each
1111 // mode from the mode's shape (modal/damping.ts), and this module gives the
1212 // analytic average of Leppington for a plate of given area and perimeter,
13-// which the full model's loss operator is fitted to. The felt of a stand is
14-// a fourth mechanism, handled with the stand (params.ts, synth.ts).
1513
1614 import { type Grid, type Material, type Plate, RHO_AIR, SPEED_OF_SOUND, bendingConstant, insideCell, massProperties, meanThickness, nodeWeight } from "./params.ts";
1715
src/main.tsmodified+6−40View file
@@ -12,8 +12,6 @@ import {
1212 massProperties,
1313 meanThickness,
1414 thicknessRange,
15- standRates,
16- standFromParams,
1715 strikePoint,
1816 strikePatch,
1917 listenerPosition,
@@ -46,8 +44,6 @@ type ViewState = "design" | "mode" | "live" | "movie";
4644
4745 const $ = <T extends HTMLElement>(id: string): T => document.getElementById(id) as T;
4846
49-/** The parameters that describe the stand; a change of any of them rebuilds the plate. */
50-const STAND_KEYS = new Set<keyof SimParams>(["standOn", "standX", "standY", "padInnerRadius", "padOuterRadius", "padThickness", "padE", "padLossFactor", "topPad", "preload"]);
5147
5248
5349 /** A recording: the listener signal and the movie of one run, from either model. */
@@ -285,29 +281,12 @@ class App {
285281 { key: "listenerAzimuth", label: "Azimuth", unit: "deg", step: 15, min: -180, max: 180, live: true },
286282 ],
287283 },
288- {
289- title: "Stand",
290- notes: true,
291- fields: [
292- { key: "standOn", label: "On a stand (off: free plate)", kind: "checkbox" },
293- { key: "padInnerRadius", label: "Felt inner radius", unit: "mm", scale: 1e-3, step: 1, min: 0 },
294- { key: "padOuterRadius", label: "Felt outer radius", unit: "mm", scale: 1e-3, step: 1, min: 1 },
295- { key: "padThickness", label: "Felt thickness", unit: "mm", scale: 1e-3, step: 1, min: 0.5 },
296- { key: "padE", label: "Felt modulus E", unit: "MPa", scale: 1e6, step: 0.1, min: 0.01 },
297- { key: "padLossFactor", label: "Felt loss factor", step: 0.05, min: 0, max: 2 },
298- { key: "topPad", label: "Top pad under the wing nut", kind: "checkbox" },
299- { key: "preload", label: "Preload of the top pad", unit: "mm", scale: 1e-3, step: 0.25, min: 0 },
300- { key: "standX", label: "Rod position x", unit: "mm", scale: 1e-3, step: 5 },
301- { key: "standY", label: "Rod position y", unit: "mm", scale: 1e-3, step: 5 },
302- ],
303- },
304284 {
305285 title: "Grid and numerics",
306286 notes: true,
307287 fields: [
308288 { key: "spacing", label: "Grid spacing", unit: "mm", scale: 1e-3, step: 1, min: 2 },
309289 { key: "modeMaxFreq", label: "Modes up to", unit: "Hz", step: 1000, min: 200 },
310- { key: "padLiftoff", label: "Pads lift off (full model)", kind: "checkbox" },
311290 { key: "fs", label: "Sample rate", kind: "select", options: [{ value: 44100, label: "44.1 kHz" }, { value: 48000, label: "48 kHz" }, { value: 96000, label: "96 kHz" }] },
312291 { key: "fps", label: "Movie frame rate", kind: "select", options: [{ value: 60, label: "60" }, { value: 120, label: "120" }, { value: 240, label: "240" }] },
313292 ],
@@ -320,13 +299,11 @@ class App {
320299
321300 private onParam(key: keyof SimParams, live: boolean): void {
322301 const p = this.form.get();
323- if (key === "spacing" || key === "modeMaxFreq" || key === "fs" || STAND_KEYS.has(key)) {
324- // The stand is part of the plate the solvers see, so it rebuilds the plate like the grid does.
302+ if (key === "spacing" || key === "modeMaxFreq" || key === "fs") {
325303 if (this.design) this.rebuildPlate(this.design);
326304 return;
327305 }
328306 if (key === "fps") { this.updateNotes(); return; }
329- if (key === "padLiftoff") { this.updateNotes(); return; }
330307 if (live && this.sim) this.sim.setListener(p);
331308 if (key.startsWith("strike") || key.startsWith("stick") || key.startsWith("contact") || key.startsWith("listener")) {
332309 this.updateMarkers();
@@ -464,7 +441,7 @@ class App {
464441 const p = this.form.get();
465442 let plate: Plate;
466443 try {
467- plate = evaluatePlate(d, p.spacing, standFromParams(p));
444+ plate = evaluatePlate(d, p.spacing);
468445 } catch (e) {
469446 this.editor.setStatus(e instanceof DesignError ? e.message : String(e), "error");
470447 return;
@@ -494,8 +471,7 @@ class App {
494471 private modesKey(): string | null {
495472 if (!this.design || !this.plate) return null;
496473 const p = this.form.get();
497- const stand = Array.from(STAND_KEYS).map((k) => String(p[k])).join(",");
498- return `${this.plateSpacing}|${Math.min(p.modeMaxFreq, p.fs / 2)}|${stand}|${this.design.source}`;
474+ return `${this.plateSpacing}|${Math.min(p.modeMaxFreq, p.fs / 2)}|${this.design.source}`;
499475 }
500476
501477 /** Installs cached modes for the current design if there are any; otherwise leaves the browser empty. */
@@ -515,7 +491,7 @@ class App {
515491 this.sys = sys; this.modes = modes; this.model = model;
516492 const sigma = modalRates(model);
517493 const t60 = (s: number) => (s > 0 ? LN1000 / s : Infinity);
518- const rows = Array.from(modes.freqs).map((f, m) => ({ freq: f, t60: t60(sigma[m]), radiation: t60(model.sigmaRad[m]), material: t60(model.sigmaMat[m]), thermoelastic: t60(model.sigmaTe[m]), felt: t60(model.padSigma[m]) }));
494+ const rows = Array.from(modes.freqs).map((f, m) => ({ freq: f, t60: t60(sigma[m]), radiation: t60(model.sigmaRad[m]), material: t60(model.sigmaMat[m]), thermoelastic: t60(model.sigmaTe[m]) }));
519495 const kind = sys.fields === 1 ? "bending modes" : "modes (bending coupled to stretching)";
520496 const when = secs === null ? "cached" : `in ${secs.toFixed(1)} s`;
521497 this.modeBrowser.setModes(rows, `${modes.count} ${kind} up to ${(fMax / 1e3).toFixed(1)} kHz, ${when} (${sys.n} unknowns${modes.missing ? `, ${modes.missing} missing` : ""}). Click a mode to see its shape.`);
@@ -533,7 +509,6 @@ class App {
533509 const { mass } = massProperties(pl);
534510 const hr = thicknessRange(pl);
535511 const ps = profileStats(pl);
536- const sr = standRates(pl);
537512 const lines: string[] = [];
538513 lines.push(`${pl.name}: ${(pl.Lx * 1e3).toFixed(0)} x ${(pl.Ly * 1e3).toFixed(0)} mm, E = ${(pl.mat.E / 1e9).toFixed(0)} GPa, nu = ${pl.mat.nu}, rho = ${pl.mat.rho} kg/m^3`);
539514 lines.push(`Thickness ${(hr.min * 1e3).toFixed(2)} to ${(hr.max * 1e3).toFixed(2)} mm, mass ${(mass * 1e3).toFixed(0)} g, grid ${g.nx} x ${g.ny} at ${(g.hx * 1e3).toFixed(2)} mm (${countInside(g)} nodes)`);
@@ -541,14 +516,6 @@ class App {
541516 const warn = ps.maxSlope > 0.1 ? ` <span class="warn">(above the shallow-shell range of about 0.1)</span>` : "";
542517 lines.push(`Profile ${(ps.peakToPeak * 1e3).toFixed(1)} mm peak to peak, largest slope ${ps.maxSlope.toFixed(3)}${warn}`);
543518 } else lines.push("Flat plate (no profile): the linear model is bending only");
544- if (pl.stand && sr) {
545- const pad = pl.stand.pad;
546- this.form.setNotes("Stand", `${pl.stand.topPad ? "Two pads" : "Bottom pad only"}, ${(pad.innerRadius * 1e3).toFixed(0)} to ${(pad.outerRadius * 1e3).toFixed(0)} mm, ${(pad.thickness * 1e3).toFixed(0)} mm felt at ${(pad.E / 1e6).toFixed(2)} MPa: ${(sr.totalKw / 1e3).toFixed(0)} kN/m in all, bounce of the rigid plate at ${sr.bounceHz.toFixed(0)} Hz, ${sr.nodes.length} pad nodes, preload ${(sr.preloadTop * 1e3).toFixed(2)} mm on top and ${(sr.preloadBottom * 1e3).toFixed(2)} mm below. The felt is viscoelastic: its modulus is the value at the bounce, and at rest it keeps ${(100 * sr.chain.kInf).toFixed(0)}% of it. The pads act on the plate nodes within the annulus whether or not the outline has a hole there.`);
547- lines.push(`On the stand (see the Stand parameters): bounce ${sr.bounceHz.toFixed(0)} Hz`);
548- } else {
549- this.form.setNotes("Stand", pl.stand ? "The felt annulus covers no plate node: move the rod or enlarge the pads." : "Free plate: the stand is off.");
550- lines.push("Free plate (the stand is off)");
551- }
552519 const geom = plateGeometry(pl);
553520 const est = (f: number) => estimatedRate(pl, geom, f);
554521 const t60 = rateText;
@@ -596,7 +563,6 @@ class App {
596563 const L = listenerPosition(p);
597564 this.renderer.setMarkers({
598565 strike: { x: sp.x, y: sp.y, radius: p.contactRadius },
599- stand: pl.stand ? { x: pl.stand.x, y: pl.stand.y, innerRadius: pl.stand.pad.innerRadius, outerRadius: pl.stand.pad.outerRadius } : null,
600566 listener: L,
601567 });
602568 }
@@ -1034,9 +1000,9 @@ class App {
10341000 for (const id of ["flat-rectangle", "bowed-rectangle", "crash-16", "crash-bell-16", "ellipse"]) {
10351001 const pr = presetById(id)!;
10361002 const d = compileDesign(pr.source);
1037- const pl = evaluatePlate(d, 0.012, standFromParams(this.form.get()));
1003+ const pl = evaluatePlate(d, 0.012);
10381004 const r = await runSelfTest(this.device, pl, this.form.get());
1039- lines.push(`${id}: max relative error ${r.maxRel.toExponential(2)}, felt state error ${r.padRel.toExponential(2)} (${r.n} nodes)`);
1005+ lines.push(`${id}: max relative error ${r.maxRel.toExponential(2)} (${r.n} nodes)`);
10401006 }
10411007 return lines.join("\n");
10421008 }
src/modal/synth.tsmodified+12−34View file
@@ -1,7 +1,7 @@
11 // Modal synthesis of a strike. Each mode m of the linearized plate is a
22 // damped oscillator q'' + 2 sigma q' + omega^2 q = f, with sigma the sum of
33 // the metal's internal friction, the thermoelastic loss, the radiation of
4-// the mode (damping.ts) and the pads' loss factor on a stand. The stick is
4+// the mode (damping.ts). The stick is
55 // the Hertz contact of the full simulation, stepped against the modal
66 // system with an exact update of every oscillator under a piecewise
77 // constant force. Once the stick has left, every mode is a damped sinusoid
@@ -15,8 +15,7 @@
1515
1616 import { type LinearSystem } from "../model/assemble.ts";
1717 import { type ModeSet } from "./modes.ts";
18-import { type Plate, type SimParams, STICK_RELEASE, listenerTerms, standRates, strikePatch, nodeMass } from "../params.ts";
19-import { feltModulus } from "../felt.ts";
18+import { type Plate, type SimParams, STICK_RELEASE, listenerTerms, strikePatch, nodeMass } from "../params.ts";
2019
2120 /** What the synthesis needs of a plate and its modes; plain arrays so that it can be posted to a worker. */
2221 export interface ModalModel {
@@ -37,8 +36,6 @@ export interface ModalModel {
3736 sigmaTe: Float64Array;
3837 /** Decay rate of each mode from its radiation into the air [1/s]; filled by damping.ts, zero until then. */
3938 sigmaRad: Float64Array;
40- /** Decay rate added by the pads' loss factor to each mode [1/s]; zero on a free plate. */
41- padSigma: Float64Array;
4239 /** Node masses on the grid [kg]. */
4340 masses: Float64Array;
4441 }
@@ -56,36 +53,17 @@ export function modalModel(pl: Plate, sys: LinearSystem, modes: ModeSet): ModalM
5653 }
5754 const masses = new Float64Array(nGrid);
5855 for (let j = 0; j < pl.grid.ny; j++) for (let i = 0; i < pl.grid.nx; i++) masses[j * pl.grid.nx + i] = nodeMass(pl, i, j);
59- // Hysteretic damping: a loss factor eta on a store of energy U gives a mode
60- // the rate sigma = eta U / omega. The pads hold the energy phi' K_pad phi / 2
61- // and the plate the rest of the mode's strain energy, omega^2 / 2 for a
62- // mass-normalized mode; the metal has its loss factor, and the felt the
63- // loss stiffness of its chain at the mode's frequency, which is eta times
64- // the storage stiffness there. The thermoelastic and the radiation rates
65- // need the shapes on the grid and are filled in by damping.ts.
66- const padSigma = new Float64Array(M), sigmaMat = new Float64Array(M);
67- const sr = standRates(pl);
56+ // Hysteretic damping: the metal's loss factor eta on the mode's strain
57+ // energy, omega^2 / 2 for a mass-normalized mode, gives the rate
58+ // sigma = eta omega / 2. The thermoelastic and the radiation rates need
59+ // the shapes on the grid and are filled in by damping.ts.
60+ const sigmaMat = new Float64Array(M);
6861 for (let m = 0; m < M; m++) {
69- let e = 0;
70- if (sr) {
71- for (const k of sr.nodes) {
72- const d = sys.nodeDof[k];
73- if (d < 0) continue;
74- const w = modes.shapes[m * sys.n + d + (sys.fields - 1)];
75- e += sr.kw[k] * w * w;
76- if (sys.fields === 3) {
77- const u = modes.shapes[m * sys.n + d], v = modes.shapes[m * sys.n + d + 1];
78- e += sr.ku[k] * (u * u + v * v);
79- }
80- }
81- }
8262 const omega = 2 * Math.PI * modes.freqs[m];
8363 if (!(omega > 0)) continue;
84- // The felt's chain has the loss stiffness L(omega) per unit of its stiffness at the bounce (felt.ts).
85- padSigma[m] = sr ? (feltModulus(sr.chain, omega).loss * e) / (2 * omega) : 0;
86- sigmaMat[m] = Math.max(0, pl.mat.eta * (omega / 2 - e / (2 * omega)));
64+ sigmaMat[m] = (pl.mat.eta * omega) / 2;
8765 }
88- return { nGrid, nRows, count: M, freqs: modes.freqs, wShapes, wRow, rowNode: Int32Array.from(rows), sigmaMat, sigmaTe: new Float64Array(M), sigmaRad: new Float64Array(M), padSigma, masses };
66+ return { nGrid, nRows, count: M, freqs: modes.freqs, wShapes, wRow, rowNode: Int32Array.from(rows), sigmaMat, sigmaTe: new Float64Array(M), sigmaRad: new Float64Array(M), masses };
8967 }
9068
9169 /** The model restricted to the modes [from, to). */
@@ -93,7 +71,7 @@ export function sliceModel(model: ModalModel, from: number, to: number): ModalMo
9371 const M = to - from;
9472 const wShapes = new Float32Array(model.nRows * M);
9573 for (let r = 0; r < model.nRows; r++) wShapes.set(model.wShapes.subarray(r * model.count + from, r * model.count + to), r * M);
96- return { ...model, count: M, freqs: model.freqs.slice(from, to), wShapes, sigmaMat: model.sigmaMat.slice(from, to), sigmaTe: model.sigmaTe.slice(from, to), sigmaRad: model.sigmaRad.slice(from, to), padSigma: model.padSigma.slice(from, to) };
74+ return { ...model, count: M, freqs: model.freqs.slice(from, to), wShapes, sigmaMat: model.sigmaMat.slice(from, to), sigmaTe: model.sigmaTe.slice(from, to), sigmaRad: model.sigmaRad.slice(from, to) };
9775 }
9876
9977 /** The transverse shape of mode m on the grid (zero outside the plate). */
@@ -103,12 +81,12 @@ export function modeShapeOnGrid(model: ModalModel, m: number): Float32Array {
10381 return w;
10482 }
10583
106-/** Per-mode decay rates: internal friction, thermoelastic, radiation and the pads, kept below the mode's frequency (underdamped). */
84+/** Per-mode decay rates: internal friction, thermoelastic and radiation, kept below the mode's frequency (underdamped). */
10785 export function modalRates(model: ModalModel): Float64Array {
10886 const sigma = new Float64Array(model.count);
10987 for (let m = 0; m < model.count; m++) {
11088 const w = 2 * Math.PI * model.freqs[m];
111- sigma[m] = Math.min(model.sigmaMat[m] + model.sigmaTe[m] + model.sigmaRad[m] + model.padSigma[m], 0.99 * w);
89+ sigma[m] = Math.min(model.sigmaMat[m] + model.sigmaTe[m] + model.sigmaRad[m], 0.99 * w);
11290 }
11391 return sigma;
11492 }
src/modal/worker.tsmodified+1−1View file
@@ -43,7 +43,7 @@ self.onmessage = (ev: MessageEvent<WorkerRequest>) => {
4343 const msg = ev.data;
4444 try {
4545 if (msg.type === "system") {
46- sys = { n: msg.n, bw: msg.bw, band: msg.band, mass: msg.mass, fields: 1, dofNode: new Int32Array(0), nodeDof: new Int32Array(0), onStand: false };
46+ sys = { n: msg.n, bw: msg.bw, band: msg.band, mass: msg.mass, fields: 1, dofNode: new Int32Array(0), nodeDof: new Int32Array(0) };
4747 post({ type: "ready" });
4848 } else if (msg.type === "job") {
4949 if (!sys) throw new Error("no system");
src/model/assemble.tsmodified+2−15View file
@@ -8,7 +8,6 @@
88
99 import {
1010 type Plate,
11- type StandRates,
1211 bendingStiffness,
1312 cellThickness,
1413 insideCell,
@@ -17,7 +16,6 @@ import {
1716 nodeMass,
1817 nodeWeight,
1918 shearStiffness,
20- standRates,
2119 } from "../params.ts";
2220 import { shapeDerivatives } from "./reference.ts";
2321
@@ -36,8 +34,6 @@ export interface LinearSystem {
3634 dofNode: Int32Array;
3735 /** First unknown of a grid node, or -1 outside the plate. */
3836 nodeDof: Int32Array;
39- /** Whether the plate is on a stand (no rigid-body modes at zero frequency). */
40- onStand: boolean;
4137 }
4238
4339 /** Whether the profile vanishes, in which case only the transverse motion is assembled. */
@@ -51,7 +47,7 @@ export function isFlatPlate(pl: Plate): boolean {
5147 * system for a flat plate (the tests use it); by default a flat plate gets
5248 * the one-field system.
5349 */
54-export function assemble(pl: Plate, fields?: 1 | 3, sr: StandRates | null = standRates(pl)): LinearSystem {
50+export function assemble(pl: Plate, fields?: 1 | 3): LinearSystem {
5551 const g = pl.grid, m = pl.mat;
5652 const { nx, ny, hx, hy } = g;
5753 const nf: 1 | 3 = fields ?? (isFlatPlate(pl) ? 1 : 3);
@@ -181,16 +177,7 @@ export function assemble(pl: Plate, fields?: 1 | 3, sr: StandRates | null = stan
181177 }
182178 }
183179
184- // ---- Stand: the pads as springs to ground ----
185- if (sr) {
186- for (const k of sr.nodes) {
187- if (nodeDof[k] < 0) continue;
188- add(wdof(k), wdof(k), sr.kw[k]);
189- if (nf === 3) { add(nodeDof[k], nodeDof[k], sr.ku[k]); add(nodeDof[k] + 1, nodeDof[k] + 1, sr.ku[k]); }
190- }
191- }
192-
193- return { n, bw, band, mass, fields: nf, dofNode, nodeDof, onStand: !!sr && sr.nodes.length > 0 };
180+ return { n, bw, band, mass, fields: nf, dofNode, nodeDof };
194181 }
195182
196183 /** y = K x from the banded storage. */
src/model/reference.tsmodified+8−103View file
@@ -10,7 +10,6 @@ import {
1010 type Grid,
1111 type Material,
1212 type Plate,
13- type StandRates,
1413 ROOM_TEMPERATURE,
1514 bendingStiffness,
1615 cellThickness,
@@ -20,10 +19,8 @@ import {
2019 nodeMass,
2120 nodeWeight,
2221 shearStiffness,
23- standRates,
2422 } from "../params.ts";
2523 import { type SigmaFit } from "../loss.ts";
26-import { feltStepConstants } from "../felt.ts";
2724
2825 export interface Fields {
2926 u: Float64Array;
@@ -39,7 +36,7 @@ export function fields(n: number): Fields {
3936 export function uniformPlate(g: Grid, mat: Material, h: number, w0?: Float64Array): Plate {
4037 const hs = new Float64Array(g.n);
4138 for (let j = 0; j < g.ny; j++) for (let i = 0; i < g.nx; i++) if (insideNode(g, i, j)) hs[j * g.nx + i] = h;
42- return { name: "uniform", Lx: (g.nx - 1) * g.hx, Ly: (g.ny - 1) * g.hy, grid: g, mat, h: hs, w0: w0 ?? new Float64Array(g.n), stand: null, temperature: ROOM_TEMPERATURE };
39+ return { name: "uniform", Lx: (g.nx - 1) * g.hx, Ly: (g.ny - 1) * g.hy, grid: g, mat, h: hs, w0: w0 ?? new Float64Array(g.n), temperature: ROOM_TEMPERATURE };
4340 }
4441
4542 const G3 = 1 / Math.sqrt(3);
@@ -58,91 +55,6 @@ export function shapeDerivatives(gp: number, hx: number, hy: number, dNdx: Float
5855 dNdy[3] = ((1 - xi) / 4) * (2 / hy);
5956 }
6057
61-/**
62- * Relaxed (static) force of the pads at a node with transverse displacement
63- * w: the felt's relaxed stiffness kInf times the stiffness at the bounce
64- * frequency, the rest of the felt's force being carried by the arms of its
65- * chain (addPadMemory). Without lift-off the pads are a linear spring of
66- * rate kInf kw. With lift-off each pad pushes only while compressed: the
67- * bottom pad by (preloadBottom - w), the top pad by (preloadTop + w); the
68- * constant part of the two, which balances the weight of the plate, is
69- * subtracted so that w = 0 is the rest state.
70- */
71-export function padForce(sr: StandRates, k: number, w: number, liftoff: boolean): number {
72- const s = sr.chain.kInf;
73- if (!liftoff) return -s * sr.kw[k] * w;
74- const kb = s * sr.kwBottom[k];
75- const hasTop = sr.kw[k] > sr.kwBottom[k] + 1e-300;
76- const bottom = kb * Math.max(sr.preloadBottom - w, 0);
77- const top = hasTop ? kb * Math.max(sr.preloadTop + w, 0) : 0;
78- const rest = kb * sr.preloadBottom - (hasTop ? kb * sr.preloadTop : 0);
79- return bottom - top - rest;
80-}
81-
82-/** Number of pads in contact with the node at transverse displacement w (both pads without lift-off). */
83-export function padsInContact(sr: StandRates, k: number, w: number, liftoff: boolean): number {
84- const hasTop = sr.kw[k] > sr.kwBottom[k] + 1e-300;
85- if (!liftoff) return hasTop ? 2 : 1;
86- return (sr.preloadBottom - w > 0 ? 1 : 0) + (hasTop && sr.preloadTop + w > 0 ? 1 : 0);
87-}
88-
89-/**
90- * The extensions of the arms of the felt's chain at every node, for the
91- * three fields: pad[(field * nArms + arm) * n + k] with field 0, 1, 2 for
92- * u, v, w. Zero is the relaxed felt, which is the initial state; the layout
93- * is that of the GPU buffer.
94- */
95-export type PadState = Float64Array;
96-export function padState(pl: Plate, sr: StandRates | null): PadState {
97- return new Float64Array(3 * (sr ? sr.chain.tau.length : 0) * pl.grid.n);
98-}
99-
100-/**
101- * Adds to the accelerations `out` the force of the arms of the felt's
102- * chain, -k_node sum_i k_i p_i on each field, with the transverse force
103- * scaled by the number of pads in contact (the shear of the felt is applied
104- * whether or not a pad has lifted, as its static part is).
105- */
106-export function addPadMemory(pl: Plate, sr: StandRates | null, f: Fields, pad: PadState, liftoff: boolean, out: Fields): void {
107- if (!sr) return;
108- const g = pl.grid, n = g.n, ch = sr.chain, nArms = ch.tau.length;
109- for (const k of sr.nodes) {
110- const mass = nodeMass(pl, k % g.nx, Math.floor(k / g.nx));
111- if (mass === 0) continue;
112- let su = 0, sv = 0, sw = 0;
113- for (let a = 0; a < nArms; a++) {
114- su += ch.k[a] * pad[a * n + k];
115- sv += ch.k[a] * pad[(nArms + a) * n + k];
116- sw += ch.k[a] * pad[(2 * nArms + a) * n + k];
117- }
118- const nc = padsInContact(sr, k, f.w[k], liftoff);
119- out.u[k] -= (sr.ku[k] * su) / mass;
120- out.v[k] -= (sr.ku[k] * sv) / mass;
121- out.w[k] -= (nc * sr.kwBottom[k] * sw) / mass;
122- }
123-}
124-
125-/**
126- * Advances the arm extensions over one step of dt in which the fields moved
127- * by `dw` (the increments of u, v, w) from the displacements `f`: exact for
128- * a linear motion within the step. At a node that had left every pad at the
129- * start of the step the transverse arms relax at once.
130- */
131-export function stepPadMemory(pl: Plate, sr: StandRates | null, dt: number, dw: Fields, f: Fields, pad: PadState, liftoff: boolean): void {
132- if (!sr) return;
133- const n = pl.grid.n, ch = sr.chain, nArms = ch.tau.length;
134- const { decay, growth } = feltStepConstants(ch, dt);
135- for (const k of sr.nodes) {
136- const lifted = padsInContact(sr, k, f.w[k], liftoff) === 0;
137- for (let a = 0; a < nArms; a++) {
138- pad[a * n + k] = pad[a * n + k] * decay[a] + dw.u[k] * growth[a];
139- pad[(nArms + a) * n + k] = pad[(nArms + a) * n + k] * decay[a] + dw.v[k] * growth[a];
140- const iw = (2 * nArms + a) * n + k;
141- pad[iw] = lifted ? 0 : pad[iw] * decay[a] + dw.w[k] * growth[a];
142- }
143- }
144-}
145-
14658 /**
14759 * Undamped accelerations (out.u, out.v, out.w) of every node for the
14860 * displacement fields f, measured from the stress-free shape pl.w0. Bending
@@ -150,10 +62,9 @@ export function stepPadMemory(pl: Plate, sr: StandRates | null, dt: number, dw:
15062 * bending stiffness D(h) at the nodes and at the cells (natural free-edge
15163 * conditions); membrane forces come from bilinear cells with 2x2 Gauss
15264 * integration, the cell's mean thickness, and the Marguerre form of the von
153- * Karman strains; the pads of the stand add springs to ground over their
154- * annulus. Nodes outside the plate get zero.
65+ * Karman strains. Nodes outside the plate get zero.
15566 */
156-export function accelerations(pl: Plate, f: Fields, out: Fields, liftoff = false, sr: StandRates | null = standRates(pl)): void {
67+export function accelerations(pl: Plate, f: Fields, out: Fields): void {
15768 const g = pl.grid, m = pl.mat, w0 = pl.w0;
15869 const { nx, ny, hx, hy } = g;
15970 const { u, v, w } = f;
@@ -230,10 +141,6 @@ export function accelerations(pl: Plate, f: Fields, out: Fields, liftoff = false
230141 (Mb(i, j - 1) - 2 * Mb(i, j) + Mb(i, j + 1)) / hy2 +
231142 (2 * (1 - nu) * (Cc(i, j) - Cc(i - 1, j) - Cc(i, j - 1) + Cc(i - 1, j - 1))) / (hx * hy);
232143 F[0] = 0; F[1] = 0; F[2] = -hx * hy * bend;
233- if (sr) {
234- F[0] -= sr.chain.kInf * sr.ku[k] * u[k]; F[1] -= sr.chain.kInf * sr.ku[k] * v[k];
235- F[2] += padForce(sr, k, w[k], liftoff);
236- }
237144 cellForce(i, j, 0, F);
238145 cellForce(i - 1, j, 1, F);
239146 cellForce(i - 1, j - 1, 2, F);
@@ -245,7 +152,7 @@ export function accelerations(pl: Plate, f: Fields, out: Fields, liftoff = false
245152 }
246153 }
247154
248-/** The loss terms of the full simulation on the plate itself: the fitted sigma model on the transverse motion and an in-plane rate. The felt's loss is in its chain (addPadMemory). */
155+/** The loss terms of the full simulation: the fitted sigma model on the transverse motion and an in-plane rate. */
249156 export interface FullLoss {
250157 sigma0: number; // [1/s], applied as 2 sigma0 on all components
251158 sigma1: number; // [m^2/s], applied as 2 sigma1 times the Laplacian of the transverse velocity
@@ -285,8 +192,8 @@ export function addDamping(pl: Plate, loss: FullLoss, vel: Fields, out: Fields):
285192 }
286193 }
287194
288-/** Discrete potential energy (bending + membrane + the relaxed pads, linear) whose negative gradient is the force in `accelerations`. */
289-export function potentialEnergy(pl: Plate, f: Fields, sr: StandRates | null = standRates(pl)): number {
195+/** Discrete potential energy (bending + membrane) whose negative gradient is the force in `accelerations`. */
196+export function potentialEnergy(pl: Plate, f: Fields): number {
290197 const g = pl.grid, m = pl.mat, w0 = pl.w0;
291198 const { nx, ny, hx, hy } = g;
292199 const { u, v, w } = f;
@@ -337,9 +244,7 @@ export function potentialEnergy(pl: Plate, f: Fields, sr: StandRates | null = st
337244 }
338245 }
339246
340- let Us = 0;
341- if (sr) for (const k of sr.nodes) Us += 0.5 * sr.chain.kInf * (sr.ku[k] * (u[k] * u[k] + v[k] * v[k]) + sr.kw[k] * w[k] * w[k]);
342- return Ub + Um + Us;
247+ return Ub + Um;
343248 }
344249
345250 /**
@@ -372,7 +277,7 @@ export function maxAngularFrequency(g: Grid, m: Material, h: number, nSmall = 12
372277 for (let k = 0; k < n; k++) {
373278 x.u[k] *= eps / norm; x.v[k] *= eps / norm; x.w[k] *= eps / norm;
374279 }
375- accelerations(pl, x, y, false, null);
280+ accelerations(pl, x, y);
376281 let dot = 0;
377282 for (let k = 0; k < n; k++) dot += x.u[k] * y.u[k] + x.v[k] * y.v[k] + x.w[k] * y.w[k];
378283 lambda = -dot / (eps * eps);
src/params.tsmodified+2−139View file
@@ -1,10 +1,8 @@
11 // Simulation parameters (strike, listener, numerics), the grid with its mask,
22 // the plate as evaluated on the grid (material, thickness field, stress-free
3-// shape, stand, loss constants), and the derived quantities both solvers need. All
3+// shape, loss constants), and the derived quantities both solvers need. All
44 // values are SI. The design itself is produced by the script in design.ts.
55
6-import { type FeltChain, feltChain } from "./felt.ts";
7-
86 export interface SimParams {
97 // Strike
108 strikeFracX: number; // strike point as a fraction of the bounding-box width from the left edge
@@ -18,21 +16,9 @@ export interface SimParams {
1816 listenerDist: number; // [m]
1917 listenerAngle: number; // [deg] angle from the plate's normal (0 = straight above the plate)
2018 listenerAzimuth: number; // [deg] rotation about the normal (0 = toward +x)
21- // Stand (a parameter of the setup, not of the design script)
22- standOn: boolean; // false for a free plate
23- standX: number; // position of the rod relative to the center of the bounding box [m]
24- standY: number;
25- padInnerRadius: number; // the felt annulus [m]
26- padOuterRadius: number;
27- padThickness: number; // [m]
28- padE: number; // Young's modulus of the felt in compression [Pa]
29- padLossFactor: number;
30- topPad: boolean; // a second pad above the plate under the wing nut
31- preload: number; // compression of the top pad at rest [m]
3219 // Numerics and output
3320 spacing: number; // requested grid spacing [m]
3421 modeMaxFreq: number; // modes up to this frequency are kept in the linear model [Hz]
35- padLiftoff: boolean; // full simulation: the pads push only while compressed
3622 fs: number; // audio sample rate [Hz]
3723 fps: number; // movie frame rate, also the batch rate of the full simulation [Hz]
3824 duration: number; // recorded duration after the impact [s]
@@ -49,19 +35,8 @@ export const defaultParams: SimParams = {
4935 listenerDist: 0.5,
5036 listenerAngle: 30,
5137 listenerAzimuth: 0,
52- standOn: true,
53- standX: 0,
54- standY: 0,
55- padInnerRadius: 7e-3,
56- padOuterRadius: 20e-3,
57- padThickness: 12e-3,
58- padE: 0.5e6,
59- padLossFactor: 0.3,
60- topPad: true,
61- preload: 1e-3,
6238 spacing: 0.008,
6339 modeMaxFreq: 16000,
64- padLiftoff: true,
6540 fs: 48000,
6641 fps: 120,
6742 duration: 4,
@@ -69,11 +44,9 @@ export const defaultParams: SimParams = {
6944
7045 export const SPEED_OF_SOUND = 343; // [m/s]
7146 export const RHO_AIR = 1.2; // [kg/m^3]
72-export const GRAVITY = 9.81; // [m/s^2]
7347 export const CFL = 0.9; // fraction of the stability limit used for the time step of the full simulation
7448 export const HIST_MARGIN = 128; // ring buffer margin of the retarded-time listener sum [samples]
7549 export const STICK_RELEASE = 5e-3; // the stick is withdrawn once the contact has been open this long [s]
76-export const PAD_POISSON = 0.3; // Poisson ratio of the felt, which sets its shear stiffness from E
7750
7851 // ---------- Grid ----------
7952
@@ -192,7 +165,7 @@ export function countInside(g: Grid): number {
192165 return c;
193166 }
194167
195-// ---------- Material, stand, losses, plate ----------
168+// ---------- Material, losses, plate ----------
196169
197170 /** Thermal constants of the metal, from which the thermoelastic damping follows (loss.ts). */
198171 export interface Thermal {
@@ -214,40 +187,6 @@ export interface Material {
214187 /** Temperature of the plate [K], which enters the thermoelastic damping. */
215188 export const ROOM_TEMPERATURE = 293;
216189
217-/** A felt pad of the stand: an annulus of felt between the plate and the stand (or the wing nut). */
218-export interface Pad {
219- innerRadius: number; // [m]
220- outerRadius: number; // [m]
221- thickness: number; // [m]
222- E: number; // Young's modulus of the felt in compression [Pa]
223- lossFactor: number; // loss factor eta of the felt
224-}
225-
226-/**
227- * The stand: a rod through the center hole at (x, y) with a felt pad below
228- * the plate and, when topPad is set, the same pad above it under the wing
229- * nut, which compresses the top pad by `preload` at rest.
230- */
231-export interface Stand {
232- x: number;
233- y: number;
234- pad: Pad;
235- topPad: boolean;
236- preload: number; // compression of the top pad at rest [m]
237-}
238-
239-/** The stand described by the parameters, or null for a free plate. */
240-export function standFromParams(p: SimParams): Stand | null {
241- if (!p.standOn) return null;
242- return {
243- x: p.standX,
244- y: p.standY,
245- pad: { innerRadius: p.padInnerRadius, outerRadius: Math.max(p.padOuterRadius, p.padInnerRadius + 1e-4), thickness: p.padThickness, E: p.padE, lossFactor: p.padLossFactor },
246- topPad: p.topPad,
247- preload: p.preload,
248- };
249-}
250-
251190 /** The plate as evaluated on the grid. */
252191 export interface Plate {
253192 name: string;
@@ -259,7 +198,6 @@ export interface Plate {
259198 h: Float64Array;
260199 /** Stress-free height at the nodes [m], mean plane removed; zero outside the plate. */
261200 w0: Float64Array;
262- stand: Stand | null;
263201 /** Temperature [K], for the thermoelastic damping. */
264202 temperature: number;
265203 }
@@ -352,81 +290,6 @@ export function thicknessRange(pl: Plate): { min: number; max: number } {
352290 return { min: isFinite(lo) ? lo : 0, max: hi };
353291 }
354292
355-// ---------- Stand as seen by the solvers ----------
356-
357-/**
358- * The pads as per-node spring rates. `kw[k]` is the transverse stiffness of
359- * the pads at node k (both pads summed), `ku[k]` the in-plane (shear)
360- * stiffness, both in N/m; `kwBottom` is the bottom pad alone, which the full
361- * simulation needs for the lift-off rule. All three are the stiffness at the
362- * bounce frequency; `chain` is the viscoelastic model of the felt (felt.ts),
363- * normalized to that stiffness, whose `kInf` is the relaxed fraction that
364- * holds a static load. `eta` is the loss factor, and `preloadTop`,
365- * `preloadBottom` are the pad compressions at rest.
366- */
367-export interface StandRates {
368- kw: Float64Array;
369- ku: Float64Array;
370- kwBottom: Float64Array;
371- nodes: number[];
372- eta: number;
373- chain: FeltChain;
374- preloadTop: number;
375- preloadBottom: number;
376- /** Total transverse stiffness of the pads [N/m]. */
377- totalKw: number;
378- /** Bounce frequency of the plate on the pads [Hz]. */
379- bounceHz: number;
380-}
381-
382-/**
383- * Per-node pad stiffness of the stand, or null for a free plate. The pads
384- * act on the plate nodes within the annulus; because a coarse grid covers
385- * the annulus with a handful of nodes, the per-node rates are scaled so
386- * that their sum equals the stiffness of the exact annulus area.
387- */
388-export function standRates(pl: Plate): StandRates | null {
389- const st = pl.stand;
390- if (!st) return null;
391- const g = pl.grid;
392- const kw = new Float64Array(g.n), ku = new Float64Array(g.n), kwBottom = new Float64Array(g.n);
393- const nodes: number[] = [];
394- const pad = st.pad;
395- const kArea = pad.E / pad.thickness; // N/m^3
396- const gArea = pad.E / (2 * (1 + PAD_POISSON)) / pad.thickness;
397- const nPads = st.topPad ? 2 : 1;
398- const annulus = Math.PI * (pad.outerRadius ** 2 - pad.innerRadius ** 2);
399- let nodeArea = 0;
400- for (let j = 0; j < g.ny; j++) {
401- for (let i = 0; i < g.nx; i++) {
402- const wgt = nodeWeight(g, i, j);
403- if (wgt === 0) continue;
404- const r = Math.hypot(nodeX(g, i) - st.x, nodeY(g, j) - st.y);
405- if (r < pad.innerRadius || r > pad.outerRadius) continue;
406- nodes.push(j * g.nx + i);
407- nodeArea += g.hx * g.hy * wgt;
408- }
409- }
410- const empty = { kw, ku, kwBottom, nodes, eta: pad.lossFactor, chain: feltChain(0, 1), preloadTop: 0, preloadBottom: 0, totalKw: 0, bounceHz: 0 };
411- if (nodes.length === 0 || nodeArea === 0) return empty;
412- const scale = annulus / nodeArea;
413- let totalKw = 0;
414- for (const k of nodes) {
415- const dA = scale * g.hx * g.hy * nodeWeight(g, k % g.nx, Math.floor(k / g.nx));
416- kwBottom[k] = kArea * dA;
417- kw[k] = nPads * kArea * dA;
418- ku[k] = nPads * gArea * dA;
419- totalKw += kw[k];
420- }
421- const { mass } = massProperties(pl);
422- const bounceHz = Math.sqrt(totalKw / mass) / (2 * Math.PI);
423- const chain = feltChain(pad.lossFactor, bounceHz);
424- const preloadTop = st.topPad ? st.preload : 0;
425- // The bottom pad carries the top pad's preload plus the weight of the plate,
426- // held by the relaxed stiffness of the felt.
427- const preloadBottom = preloadTop + (mass * GRAVITY) / (kArea * annulus * chain.kInf);
428- return { kw, ku, kwBottom, nodes, eta: pad.lossFactor, chain, preloadTop, preloadBottom, totalKw, bounceHz };
429-}
430293
431294 // ---------- Patches ----------
432295
src/presets.tsmodified+11−18View file
@@ -9,8 +9,7 @@ export interface Preset {
99 source: string;
1010 }
1111
12-const flatRectangle = `% A flat rectangular sheet of stainless steel (turn the stand off in the
13-% Stand parameters to let it hang freely).
12+const flatRectangle = `% A flat rectangular sheet of stainless steel, hanging freely.
1413 % The script is an ordinary MATLAB program that fills in a struct named
1514 % design. Lengths are meters; the constants at the top are the only units.
1615 % A field can be a number or a function handle @(x, y, r, theta) of the
@@ -20,7 +19,7 @@ const flatRectangle = `% A flat rectangular sheet of stainless steel (turn the s
2019 % functions at the bottom; the other presets show how. The decay of the
2120 % sound is not prescribed: it follows from the loss factor of the metal, its
2221 % thermal constants (thermoelastic damping), the radiation of each mode into
23-% the air, and the felt of a stand.
22+% the air.
2423 inch = 0.0254; % m
2524 mm = 1e-3; % m
2625
@@ -38,8 +37,7 @@ design.thickness = 0.79*mm;
3837 `;
3938
4039 const aluminumSquare = `% A flat 12 inch square of 6061-T6 aluminum, 1.02 mm thick, with a center
41-% hole for the rod of a stand (the stand itself is set in the Stand
42-% parameters beside the script). Aluminum is a third the density of steel and a third as stiff, so
40+% hole. Aluminum is a third the density of steel and a third as stiff, so
4341 % the bending frequencies land near those of a steel plate of the same size
4442 % and thickness. Its high thermal expansion and conductivity give it the
4543 % strongest thermoelastic damping of the three metals here.
@@ -88,11 +86,7 @@ v = sag * (1 - t.^2) .* (1 + offset*t);
8886 end
8987 `;
9088
91-const flatDisc = `% A flat 16 inch bronze disc with a center hole for the rod of a stand. The
92-% stand itself (the felt pads, the preload, the position of the rod) is set
93-% in the Stand parameters beside the script, so that any design can be put
94-% on it or taken off it; the pads hold the disc level and let it bounce and
95-% rock.
89+const flatDisc = `% A flat 16 inch bronze disc with a center hole, hanging freely.
9690 inch = 0.0254; mm = 1e-3;
9791 D = 16*inch; hole = 13*mm;
9892
@@ -111,7 +105,7 @@ design.thickness = 1.0*mm;
111105 `;
112106
113107 const crash16 = `% A 16 inch flat crash: bronze, lathed to a taper from 1.2 mm at the
114-% center to 0.7 mm at the edge, with a center hole for the stand. Flat, so
108+% center to 0.7 mm at the edge, with a center hole. Flat, so
115109 % the linear modes are pure bending modes; the tapered profile is what
116110 % makes it a crash rather than a ride.
117111 inch = 0.0254; mm = 1e-3;
@@ -138,10 +132,9 @@ end
138132 `;
139133
140134 const ride20 = `% A 20 inch flat ride: thicker bronze, 2.2 mm at the center to 1.4 mm at
141-% the edge, with a center hole for the stand. The thickness raises every
142-% frequency in proportion, and the ride has a clearer ping and a longer
143-% decay than the crash; try a small preload in the Stand parameters for a
144-% loose wing nut.
135+% the edge, with a center hole. The thickness raises every frequency in
136+% proportion, and the ride has a clearer ping and a longer decay than the
137+% crash.
145138 inch = 0.0254; mm = 1e-3;
146139 D = 20*inch; R = D/2; hole = 13*mm;
147140
@@ -219,9 +212,9 @@ design.inside = @(x, y, r, theta) (x/(A/2)).^2 + (y/(B/2)).^2 <= 1;
219212 design.thickness = 1.0*mm;
220213 `;
221214
222-const thunderSheet = `% A thunder sheet: a large thin steel sheet, hung from a corner in reality
223-% (turn the stand off). Its modes
224-% are dense and low, and the full simulation turns a hard hit into a roar.
215+const thunderSheet = `% A thunder sheet: a large thin steel sheet, hung from a corner in reality.
216+% Its modes are dense and low, and the full simulation turns a hard hit into
217+% a roar.
225218 % Long at the default spacing; a coarser grid keeps it interactive.
226219 inch = 0.0254; mm = 1e-3;
227220
src/render/renderer.tsmodified+6−38View file
@@ -1,7 +1,7 @@
11 // Drawing of the plate: a 2D color map, or (WebGPU only) a lit 3D height
22 // field of the displayed field plus, optionally, the stress-free profile,
33 // colored by the height or by the thickness. Both views paint the strike
4-// point, the stand pad and the listener as markers. Without WebGPU the map is
4+// point and the listener as markers. Without WebGPU the map is
55 // drawn with the 2D canvas API and the 3D mode falls back to the map.
66 //
77 // On the GPU the plate is one indexed mesh: a static vertex buffer of node
@@ -32,8 +32,6 @@ export interface ViewOptions {
3232 export interface Markers {
3333 /** Plate coordinates [m]. */
3434 strike: { x: number; y: number; radius: number } | null;
35- /** Pad annulus of the stand [m]. */
36- stand: { x: number; y: number; innerRadius: number; outerRadius: number } | null;
3735 /** Relative to the plate center [m], z along the normal. */
3836 listener: { x: number; y: number; z: number } | null;
3937 }
@@ -42,8 +40,6 @@ type RGB = [number, number, number];
4240
4341 const BG: RGB = [0.11, 0.12, 0.14];
4442 const COLOR_STRIKE: RGB = [0.95, 0.78, 0.12];
45-const COLOR_STAND: RGB = [0.1, 0.1, 0.1];
46-const COLOR_POST: RGB = [0.55, 0.55, 0.58];
4743 const COLOR_LINE: RGB = [0.55, 0.55, 0.55];
4844 const COLOR_MIC_HEAD: RGB = [0.32, 0.32, 0.35];
4945 const COLOR_MIC_BODY: RGB = [0.72, 0.73, 0.76];
@@ -216,7 +212,7 @@ export class PlateRenderer {
216212 private plate: Plate | null = null;
217213 private field: Float32Array | Float64Array | null = null;
218214 private opts: ViewOptions = { mode: "2d", color: "height", showProfile: true, exaggeration: 1, range: 0 };
219- private markers: Markers = { strike: null, stand: null, listener: null };
215+ private markers: Markers = { strike: null, listener: null };
220216
221217 // Derived per plate
222218 private hRange = { min: 0, max: 1 };
@@ -537,19 +533,6 @@ export class PlateRenderer {
537533 return wsum > 0 ? s / wsum : fallback;
538534 }
539535
540- /** Mean displayed height of the plate nodes within `r` of (cx, cy), or 0 when there are none. */
541- private meanZWithin(cx: number, cy: number, r: number): number {
542- const g = this.plate!.grid;
543- let s = 0, c = 0;
544- for (let j = 0; j < g.ny; j++) {
545- for (let i = 0; i < g.nx; i++) {
546- if (!insideNode(g, i, j) || Math.hypot(nodeX(g, i) - cx, nodeY(g, j) - cy) > r) continue;
547- s += this.z[j * g.nx + i]; c++;
548- }
549- }
550- return c > 0 ? s / c : 0;
551- }
552-
553536 // ---------- Listener placement ----------
554537
555538 /** Where the listener is drawn in the 2D map: its in-plane direction, at a saturated distance from the center. */
@@ -580,19 +563,8 @@ export class PlateRenderer {
580563 const R = this.boxRadius;
581564 const lw = 0.012 * R; // width of the ring lines
582565 const lift = 0.004 * R; // markers float this far above the surface
583- const { strike, stand, listener } = this.markers;
584- const zStand = stand ? this.meanZWithin(stand.x, stand.y, stand.outerRadius) : 0;
585- const zAt = is3d ? (x: number, y: number) => this.surfaceZ(x, y, zStand) + lift : () => lift;
586- if (stand) {
587- const dark = pack(COLOR_STAND);
588- if (stand.innerRadius > 0) mesh.ring(stand.x, stand.y, stand.innerRadius, lw / 2, zAt, dark);
589- mesh.ring(stand.x, stand.y, stand.outerRadius, lw / 2, zAt, dark);
590- if (is3d) {
591- const ro = stand.outerRadius;
592- const rod = stand.innerRadius > 0 ? 0.6 * stand.innerRadius : 0.3 * ro;
593- mesh.cylinder([stand.x, stand.y, zStand - 2 * ro], [stand.x, stand.y, zStand + 0.5 * ro], rod, 20, pack(COLOR_POST));
594- }
595- }
566+ const { strike, listener } = this.markers;
567+ const zAt = is3d ? (x: number, y: number) => this.surfaceZ(x, y, 0) + lift : () => lift;
596568 if (strike) mesh.ring(strike.x, strike.y, Math.max(strike.radius, 2 * lw), lw / 2, zAt, pack(COLOR_STRIKE));
597569 if (listener) {
598570 if (is3d) {
@@ -600,7 +572,7 @@ export class PlateRenderer {
600572 if (l) {
601573 const [e1, e2] = frame(l.u);
602574 const head = 0.035 * R, body = 0.025 * R, len = 0.06 * R;
603- const zc = this.surfaceZ(0, 0, zStand);
575+ const zc = this.surfaceZ(0, 0, 0);
604576 const start: Vec3 = [0, 0, zc + lift];
605577 const end = add(l.p, scale(l.u, -head));
606578 if (dot(end, l.u) > dot(start, l.u)) mesh.cylinder(start, end, 0.003 * R, 6, pack(COLOR_LINE));
@@ -851,11 +823,7 @@ export class PlateRenderer {
851823 ctx.strokeStyle = css(color); ctx.lineWidth = lw;
852824 ctx.stroke();
853825 };
854- const { strike, stand } = this.markers;
855- if (stand) {
856- if (stand.innerRadius > 0) ring(stand.x, stand.y, stand.innerRadius, COLOR_STAND);
857- ring(stand.x, stand.y, stand.outerRadius, COLOR_STAND);
858- }
826+ const { strike } = this.markers;
859827 if (strike) ring(strike.x, strike.y, Math.max(strike.radius, 0.024 * R), COLOR_STRIKE);
860828 const l = this.listener2d();
861829 if (l) {
src/sim/selftest.tsmodified+6−25View file
@@ -1,10 +1,9 @@
11 // In-browser check of the WGSL step kernel against the float64 reference:
22 // one step from a smooth random state of a few thicknesses amplitude, with
3-// the losses and the stand (lift-off and the felt's chain included),
4-// compared as accelerations, and the update of the felt's state.
3+// the losses, compared as accelerations.
54
6-import { type Plate, type SimParams, meanThickness, nodeX, nodeY, insideNode, standRates } from "../params.ts";
7-import { accelerations, addDamping, addPadMemory, fields, fullLoss, padState, stepPadMemory } from "../model/reference.ts";
5+import { type Plate, type SimParams, meanThickness, nodeX, nodeY, insideNode } from "../params.ts";
6+import { accelerations, addDamping, fields, fullLoss } from "../model/reference.ts";
87 import { fitSigma, inPlaneRate } from "../loss.ts";
98 import { Simulation } from "./simulation.ts";
109
@@ -12,8 +11,6 @@ export interface SelfTestResult {
1211 n: number;
1312 /** Largest acceleration error over all nodes and components relative to the largest reference acceleration. */
1413 maxRel: number;
15- /** Largest error of the felt's arm extensions after the step relative to their largest value (0 on a free plate). */
16- padRel: number;
1714 dt: number;
1815 }
1916
@@ -40,7 +37,6 @@ export async function runSelfTest(device: GPUDevice, pl: Plate, params: SimParam
4037 const g = pl.grid, n = g.n;
4138 const hm = meanThickness(pl);
4239 const fit = fitSigma(pl);
43- const sr = standRates(pl);
4440 const loss = fullLoss(fit, inPlaneRate(pl));
4541 const sim = new Simulation(device, pl, params, loss);
4642 const amp = 2.5 * hm;
@@ -51,24 +47,17 @@ export async function runSelfTest(device: GPUDevice, pl: Plate, params: SimParam
5147 pos32[k] = f.u[k]; pos32[n + k] = f.v[k]; pos32[2 * n + k] = f.w[k];
5248 vel32[k] = vel.u[k]; vel32[n + k] = vel.v[k]; vel32[2 * n + k] = vel.w[k];
5349 }
54- // The felt's arms start with smooth random extensions of a fraction of the amplitude.
55- const pad = padState(pl, sr);
56- const nArms = sr ? sr.chain.tau.length : 0;
57- for (let a = 0; a < 3 * nArms; a++) pad.set(smoothField(pl, 0.3 * amp, 20 + a), a * n);
58- const pad32 = new Float32Array(pad);
5950 // Use the float32-rounded state as the reference input.
6051 for (let k = 0; k < n; k++) {
6152 f.u[k] = pos32[k]; f.v[k] = pos32[n + k]; f.w[k] = pos32[2 * n + k];
6253 vel.u[k] = vel32[k]; vel.v[k] = vel32[n + k]; vel.w[k] = vel32[2 * n + k];
6354 }
64- pad.set(pad32);
65- sim.writeState(pos32, vel32, pad32);
55+ sim.writeState(pos32, vel32);
6656 await sim.stepOnce();
6757 const st = await sim.readState();
6858 const ref = fields(n);
69- accelerations(pl, f, ref, params.padLiftoff, sr);
59+ accelerations(pl, f, ref);
7060 addDamping(pl, loss, vel, ref);
71- addPadMemory(pl, sr, f, pad, params.padLiftoff, ref);
7261 let maxRef = 0, maxErr = 0;
7362 const dt = sim.dt;
7463 for (let k = 0; k < n; k++) {
@@ -76,14 +65,6 @@ export async function runSelfTest(device: GPUDevice, pl: Plate, params: SimParam
7665 maxRef = Math.max(maxRef, Math.abs(ref.u[k]), Math.abs(ref.v[k]), Math.abs(ref.w[k]));
7766 maxErr = Math.max(maxErr, Math.abs(au - ref.u[k]), Math.abs(av - ref.v[k]), Math.abs(aw - ref.w[k]));
7867 }
79- // The felt's state after the step, from the GPU's own new velocities.
80- let padRef = 0, padErr = 0;
81- if (sr) {
82- const dw = fields(n);
83- for (let k = 0; k < n; k++) { dw.u[k] = dt * st.vel[k]; dw.v[k] = dt * st.vel[n + k]; dw.w[k] = dt * st.vel[2 * n + k]; }
84- stepPadMemory(pl, sr, dt, dw, f, pad, params.padLiftoff);
85- for (let i = 0; i < pad.length; i++) { padRef = Math.max(padRef, Math.abs(pad[i])); padErr = Math.max(padErr, Math.abs(pad[i] - st.pad[i])); }
86- }
8768 sim.destroy();
88- return { n, maxRel: maxErr / maxRef, padRel: padRef > 0 ? padErr / padRef : 0, dt };
69+ return { n, maxRel: maxErr / maxRef, dt };
8970 }
src/sim/shaders/common.wgslmodified+8−21View file
@@ -31,25 +31,17 @@ struct Params {
3131 patchJ1: u32,
3232 patchCX: f32, // mass-weighted center of the contact patch [m]
3333 patchCY: f32,
34- standOn: u32, // 1 when the plate sits on a stand, 0 when it is free
35- liftoff: u32, // 1 when the pads push only while compressed
36- preloadTop: f32, // compression of the top pad at rest [m] (0 without a top pad)
37- preloadBottom: f32, // compression of the bottom pad at rest [m]
38- kInf: f32, // relaxed stiffness of the felt, per unit of its stiffness at the bounce frequency
39- hasTop: u32, // 1 when there is a top pad
40- nArms: u32, // arms of the felt's chain (at most MAX_ARMS)
34+ pad0: u32,
4135 pad1: u32,
4236 pad2: u32,
4337 pad3: u32,
44- // The chain of the felt (felt.ts), four arms per vec4: decay factors
45- // exp(-dt / tau) in [0..2], growth factors (tau / dt)(1 - decay) in
46- // [3..5], and arm stiffnesses per unit of pad stiffness in [6..8].
47- chain: array<vec4f, 9>,
38+ pad4: u32,
39+ pad5: u32,
40+ pad6: u32,
41+ pad7: u32,
42+ pad8: u32,
43+ pad9: u32,
4844 };
49-const MAX_ARMS = 12u;
50-fn armDecay(a: u32) -> f32 { return P.chain[a >> 2u][a & 3u]; }
51-fn armGrowth(a: u32) -> f32 { return P.chain[3u + (a >> 2u)][a & 3u]; }
52-fn armK(a: u32) -> f32 { return P.chain[6u + (a >> 2u)][a & 3u]; }
5345
5446 // Per-step parameters, addressed with a dynamic uniform offset.
5547 struct StepParams {
@@ -72,13 +64,8 @@ fn nodeY(j: i32) -> f32 { return (f32(j) - 0.5 * f32(P.ny - 1u)) * P.hy; }
7264
7365 // Static per-node data, written once by the host: x = stress-free shape w0 [m],
7466 // y = quadrature and lumped-mass weight (the fraction of the node's four cells
75-// that belong to the plate; zero outside), z = thickness [m], w = bottom-pad
76-// stiffness at the node [N/m] (zero away from the pads; the top pad, when
77-// present, has the same rate).
67+// that belong to the plate; zero outside), z = thickness [m], w unused.
7868 @group(0) @binding(8) var<storage, read> nd: array<vec4f>;
79-// Extensions of the arms of the felt's chain: padQ[(field * nArms + arm) * n + k]
80-// with field 0, 1, 2 for u, v, w; zero is the relaxed felt.
81-@group(0) @binding(9) var<storage, read_write> padQ: array<f32>;
8269 fn nodeWeight(i: i32, j: i32) -> f32 {
8370 if (i < 0 || i >= i32(P.nx) || j < 0 || j >= i32(P.ny)) { return 0.0; }
8471 return nd[u32(j) * P.nx + u32(i)].y;
src/sim/shaders/step.wgslmodified+7−68View file
@@ -1,6 +1,5 @@
11 // One explicit time step (symplectic Euler) of the von Karman plate with a
2-// thickness field, a stress-free profile, the drumstick, and the pads of the
3-// stand.
2+// thickness field, a stress-free profile, and the drumstick.
43 //
54 // State layout (structure of arrays, n = nx * ny nodes):
65 // pos: [u | v | w | wLo] wLo carries the rounding error of the w update
@@ -168,29 +167,6 @@ fn cellForce(ci: i32, cj: i32, a: u32) -> vec3f {
168167 return vec3f(cellF[base], cellF[base + 1u], cellF[base + 2u]);
169168 }
170169
171-// Relaxed force of the pads on a node with bottom-pad rate kb (at the bounce
172-// frequency) and transverse displacement w (same rule as padForce in
173-// reference.ts); the rest of the felt's force is in the arms of its chain.
174-// The constant part that balances the weight of the plate is subtracted,
175-// so w = 0 is rest.
176-fn padForce(kb0: f32, w: f32) -> f32 {
177- let kb = P.kInf * kb0;
178- if (P.liftoff == 0u) {
179- if (P.hasTop != 0u) { return -2.0 * kb * w; }
180- return -kb * w;
181- }
182- var f = kb * max(P.preloadBottom - w, 0.0) - kb * P.preloadBottom;
183- if (P.hasTop != 0u) { f -= kb * max(P.preloadTop + w, 0.0) - kb * P.preloadTop; }
184- return f;
185-}
186-// Number of pads in contact with a node at transverse displacement w.
187-fn padsInContact(w: f32) -> f32 {
188- if (P.liftoff == 0u) { return select(1.0, 2.0, P.hasTop != 0u); }
189- var nc = select(0.0, 1.0, P.preloadBottom - w > 0.0);
190- if (P.hasTop != 0u && P.preloadTop + w > 0.0) { nc += 1.0; }
191- return nc;
192-}
193-
194170 @compute @workgroup_size(16, 16)
195171 fn main(@builtin(local_invocation_id) lid: vec3u, @builtin(workgroup_id) wid: vec3u) {
196172 let n = P.n;
@@ -259,27 +235,6 @@ fn main(@builtin(local_invocation_id) lid: vec3u, @builtin(workgroup_id) wid: ve
259235 F += cellForce(i, j - 1, 3u);
260236 let pOld = vec3f(posIn[k], posIn[n + k], posIn[2u * n + k]);
261237 let vel = vec3f(velIn[k], velIn[n + k], velIn[2u * n + k]);
262- // The pads: the relaxed spring to ground on the transverse motion
263- // (compression-only with lift-off) and the felt's shear on the in-plane
264- // motion, plus the arms of the felt's chain, whose extensions padQ carry
265- // the memory of the felt (felt.ts, addPadMemory in reference.ts).
266- let kb = nd[k].w;
267- let onPad = P.standOn != 0u && kb > 0.0;
268- var nc = 0.0;
269- if (onPad) {
270- var nPads = 1.0; if (P.hasTop != 0u) { nPads = 2.0; }
271- let ku = nPads * kb * 0.3846153846; // shear: G / E = 1 / (2 (1 + 0.3))
272- F.x -= P.kInf * ku * pOld.x;
273- F.y -= P.kInf * ku * pOld.y;
274- F.z += padForce(kb, pOld.z);
275- nc = padsInContact(pOld.z);
276- var arm = vec3f(0.0);
277- for (var a = 0u; a < P.nArms; a++) {
278- let base = a * n + k;
279- arm += armK(a) * vec3f(padQ[base], padQ[P.nArms * n + base], padQ[2u * P.nArms * n + base]);
280- }
281- F -= vec3f(ku, ku, nc * kb) * arm;
282- }
283238 let mass = P.rho * TH(i, j) * P.hx * P.hy * weight;
284239 var acc = F * (1.0 / mass);
285240
@@ -298,15 +253,12 @@ fn main(@builtin(local_invocation_id) lid: vec3u, @builtin(workgroup_id) wid: ve
298253 // contact has been released for STICK_RELEASE the stick is withdrawn.
299254 if ((SP.flags & 1u) != 0u) {
300255 let sIn = 8u * SP.parity;
301- if (P.standOn == 0u) {
302- // Free plate: inertia relief. The force of the previous step is applied
303- // to the whole plate as the rigid acceleration field that balances the
304- // patch load and its moments about the center, so the plate gains no net
305- // momentum (the thread takes it up in reality). A plate on a stand takes
306- // the reaction through the pads.
307- let Fprev = stick[sIn + 2u];
308- acc.z += Fprev * (1.0 / P.mass + nodeY(j) * P.patchCY / P.Ix + nodeX(i) * P.patchCX / P.Iy);
309- }
256+ // Inertia relief of the free plate. The force of the previous step is
257+ // applied to the whole plate as the rigid acceleration field that balances
258+ // the patch load and its moments about the center, so the plate gains no
259+ // net momentum (the thread takes it up in reality).
260+ let Fprev = stick[sIn + 2u];
261+ acc.z += Fprev * (1.0 / P.mass + nodeY(j) * P.patchCY / P.Ix + nodeX(i) * P.patchCX / P.Iy);
310262 if (inContactPatch(i, j)) {
311263 let sx = stick[sIn]; let sv = stick[sIn + 1u]; let sMax = stick[sIn + 3u];
312264 var released = stick[sIn + 4u];
@@ -344,19 +296,6 @@ fn main(@builtin(local_invocation_id) lid: vec3u, @builtin(workgroup_id) wid: ve
344296 let vNew = vel + P.dt * acc;
345297 if (SP.audioSlot >= 0) { hist[u32(SP.audioSlot) * n + k] = vNew.z; }
346298 var pNew = pOld + P.dt * vNew;
347- // The arms of the felt relax over the step while the node moves by dt vNew
348- // (exact for a linear motion within the step); a node that has left every
349- // pad relaxes its transverse arms at once.
350- if (onPad) {
351- let dw = P.dt * vNew;
352- for (var a = 0u; a < P.nArms; a++) {
353- let base = a * n + k;
354- let iu = base; let iv = P.nArms * n + base; let iw = 2u * P.nArms * n + base;
355- padQ[iu] = padQ[iu] * armDecay(a) + dw.x * armGrowth(a);
356- padQ[iv] = padQ[iv] * armDecay(a) + dw.y * armGrowth(a);
357- padQ[iw] = select(0.0, padQ[iw] * armDecay(a) + dw.z * armGrowth(a), nc > 0.0);
358- }
359- }
360299 // Compensated update of w: the increment dt * v can be far below the
361300 // float32 resolution of w, so the rounding error is carried separately.
362301 let y = P.dt * vNew.z + posIn[3u * n + k];
src/sim/simulation.tsmodified+11−41View file
@@ -11,11 +11,9 @@
1111 import {
1212 type Plate,
1313 type SimParams,
14- type StandRates,
1514 CFL,
1615 HIST_MARGIN,
1716 massProperties,
18- standRates,
1917 strikePatch,
2018 strikePoint,
2119 thicknessRange,
@@ -23,7 +21,6 @@ import {
2321 listenerTerms,
2422 } from "../params.ts";
2523 import { type FullLoss, maxAngularFrequency } from "../model/reference.ts";
26-import { FELT_ARMS, feltStepConstants } from "../felt.ts";
2724 import commonWgsl from "./shaders/common.wgsl?raw";
2825 import stepWgsl from "./shaders/step.wgsl?raw";
2926 import audioWgsl from "./shaders/audio.wgsl?raw";
@@ -53,7 +50,7 @@ export interface SimStatus {
5350 stick: StickState;
5451 }
5552
56-const PARAMS_LEN = 40 + 9 * 4; // Params: 40 scalars and the 9 vec4 of the felt's chain
53+const PARAMS_LEN = 40;
5754 const STEP_LEN = 12;
5855 const STICK_LEN = 16;
5956 /** Nodes per workgroup along each axis in step.wgsl. */
@@ -75,7 +72,6 @@ export class Simulation {
7572 readonly omegaMax: number;
7673 readonly histLen: number;
7774 readonly mass: number;
78- readonly stand: StandRates | null;
7975
8076 audio: Float32Array = new Float32Array(0);
8177 samples = 0;
@@ -101,7 +97,6 @@ export class Simulation {
10197 private hist: GPUBuffer;
10298 private stickBuf: GPUBuffer;
10399 private nd: GPUBuffer;
104- private padQ: GPUBuffer;
105100 private lw: GPUBuffer;
106101 private audioOut: GPUBuffer;
107102 private frameBuf: GPUBuffer;
@@ -138,18 +133,13 @@ export class Simulation {
138133 this.vel = [f32(3 * n), f32(3 * n)];
139134 this.hist = f32(this.histLen * n);
140135 this.stickBuf = f32(STICK_LEN);
141- this.stand = standRates(plate);
142136 const weights = nodeWeights(g);
143137 const packed = new Float32Array(4 * n);
144138 for (let k = 0; k < n; k++) {
145139 packed[4 * k] = plate.w0[k]; packed[4 * k + 1] = weights[k]; packed[4 * k + 2] = plate.h[k];
146- packed[4 * k + 3] = this.stand ? this.stand.kwBottom[k] : 0;
147140 }
148141 this.nd = f32(4 * n);
149142 device.queue.writeBuffer(this.nd, 0, packed);
150- const nArms = this.stand ? this.stand.chain.tau.length : 0;
151- if (nArms > FELT_ARMS) throw new Error("the felt's chain has more arms than the kernel holds");
152- this.padQ = f32(3 * nArms * n);
153143 this.lw = f32(2 * n);
154144 this.audioOut = f32(this.samplesPerBatch);
155145 this.frameBuf = f32(n);
@@ -180,17 +170,6 @@ export class Simulation {
180170 pu[20] = cp.nearest; pf[21] = params.contactStiffness; pf[22] = params.contactExponent; pf[23] = params.stickMass;
181171 pu[24] = i0; pu[25] = i1; pu[26] = j0; pu[27] = j1;
182172 pf[28] = cp.cx; pf[29] = cp.cy;
183- const st = this.stand;
184- pu[30] = st && st.nodes.length > 0 ? 1 : 0;
185- pu[31] = params.padLiftoff ? 1 : 0;
186- pf[32] = st ? st.preloadTop : 0; pf[33] = st ? st.preloadBottom : 0;
187- pf[34] = st ? st.chain.kInf : 1;
188- pu[35] = plate.stand?.topPad ? 1 : 0;
189- pu[36] = nArms;
190- if (st) {
191- const { decay, growth } = feltStepConstants(st.chain, this.dt);
192- for (let a = 0; a < nArms; a++) { pf[40 + a] = decay[a]; pf[52 + a] = growth[a]; pf[64 + a] = st.chain.k[a]; }
193- }
194173 this.paramsBuf = device.createBuffer({ size: pf.byteLength, usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST });
195174 device.queue.writeBuffer(this.paramsBuf, 0, pf);
196175
@@ -205,7 +184,7 @@ export class Simulation {
205184 const pipeline = (module: GPUShaderModule, layout: GPUBindGroupLayout) =>
206185 device.createComputePipeline({ layout: device.createPipelineLayout({ bindGroupLayouts: [layout] }), compute: { module, entryPoint: "main" } });
207186 const stepLayout = device.createBindGroupLayout({
208- entries: [entry(0, "uniform"), entry(1, "uniform", true), entry(2, "read-only-storage"), entry(3, "read-only-storage"), entry(4, "storage"), entry(5, "storage"), entry(6, "storage"), entry(7, "storage"), entry(8, "read-only-storage"), entry(9, "storage")],
187+ entries: [entry(0, "uniform"), entry(1, "uniform", true), entry(2, "read-only-storage"), entry(3, "read-only-storage"), entry(4, "storage"), entry(5, "storage"), entry(6, "storage"), entry(7, "storage"), entry(8, "read-only-storage")],
209188 });
210189 this.stepPipe = pipeline(shader(stepWgsl), stepLayout);
211190 const audioLayout = device.createBindGroupLayout({
@@ -229,7 +208,6 @@ export class Simulation {
229208 { binding: 6, resource: { buffer: this.hist } },
230209 { binding: 7, resource: { buffer: this.stickBuf } },
231210 { binding: 8, resource: { buffer: this.nd } },
232- { binding: 9, resource: { buffer: this.padQ } },
233211 ],
234212 }),
235213 );
@@ -406,7 +384,7 @@ export class Simulation {
406384
407385 // ---------- Testing support ----------
408386
409- /** Advances the state by exactly one free step (no drumstick, no audio) with the losses and the stand. */
387+ /** Advances the state by exactly one free step (no drumstick, no audio) with the losses. */
410388 async stepOnce(): Promise<void> {
411389 this.writeStepParams(0, { flags: 0, slot: -1, outIdx: -1 });
412390 this.device.queue.writeBuffer(this.stepBuf, 0, this.stepData, 0, this.stepStride);
@@ -421,39 +399,31 @@ export class Simulation {
421399 this.totalSteps++;
422400 await this.device.queue.onSubmittedWorkDone();
423401 }
424- /** Length of the felt's state: 3 fields times the arms of the chain times the nodes (0 on a free plate). */
425- get padStateLength(): number {
426- return 3 * (this.stand ? this.stand.chain.tau.length : 0) * this.plate.grid.n;
427- }
428- /** Overwrites the state: pos as [u | v | w | wLo], vel as [du | dv | dw], and the extensions of the felt's arms (relaxed when omitted). */
429- writeState(pos: Float32Array<ArrayBuffer>, vel: Float32Array<ArrayBuffer>, pad?: Float32Array<ArrayBuffer>): void {
402+ /** Overwrites the state: pos as [u | v | w | wLo], vel as [du | dv | dw]. */
403+ writeState(pos: Float32Array<ArrayBuffer>, vel: Float32Array<ArrayBuffer>): void {
430404 const n = this.plate.grid.n;
431405 if (pos.length !== 4 * n || vel.length !== 3 * n) throw new Error("state arrays have the wrong length");
432406 this.device.queue.writeBuffer(this.pos[this.cur], 0, pos);
433407 this.device.queue.writeBuffer(this.vel[this.cur], 0, vel);
434- const np = this.padStateLength;
435- if (pad !== undefined && pad.length !== np) throw new Error("the felt state has the wrong length");
436- if (np > 0) this.device.queue.writeBuffer(this.padQ, 0, pad ?? new Float32Array(np));
437408 }
438- /** Reads back the full state as [u | v | w | wLo], [du | dv | dw], and the felt's arm extensions. */
439- async readState(): Promise<{ pos: Float32Array; vel: Float32Array; pad: Float32Array }> {
440- const n = this.plate.grid.n, np = this.padStateLength;
441- const stg = this.device.createBuffer({ size: (7 * n + np) * 4, usage: GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST });
409+ /** Reads back the full state as [u | v | w | wLo] and [du | dv | dw]. */
410+ async readState(): Promise<{ pos: Float32Array; vel: Float32Array }> {
411+ const n = this.plate.grid.n;
412+ const stg = this.device.createBuffer({ size: 7 * n * 4, usage: GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST });
442413 const enc = this.device.createCommandEncoder();
443414 enc.copyBufferToBuffer(this.pos[this.cur], 0, stg, 0, 4 * n * 4);
444415 enc.copyBufferToBuffer(this.vel[this.cur], 0, stg, 4 * n * 4, 3 * n * 4);
445- if (np > 0) enc.copyBufferToBuffer(this.padQ, 0, stg, 7 * n * 4, np * 4);
446416 this.device.queue.submit([enc.finish()]);
447417 await stg.mapAsync(GPUMapMode.READ);
448418 const all = new Float32Array(stg.getMappedRange()).slice();
449419 stg.unmap();
450420 stg.destroy();
451- return { pos: all.subarray(0, 4 * n), vel: all.subarray(4 * n, 7 * n), pad: all.subarray(7 * n) };
421+ return { pos: all.subarray(0, 4 * n), vel: all.subarray(4 * n) };
452422 }
453423
454424 destroy(): void {
455425 this.inFlight = null;
456- for (const b of [...this.pos, ...this.vel, this.hist, this.stickBuf, this.nd, this.padQ, this.lw, this.audioOut, this.frameBuf, this.paramsBuf, this.stepBuf, ...this.staging]) b.destroy();
426+ for (const b of [...this.pos, ...this.vel, this.hist, this.stickBuf, this.nd, this.lw, this.audioOut, this.frameBuf, this.paramsBuf, this.stepBuf, ...this.staging]) b.destroy();
457427 }
458428 get fs(): number { return this.params.fs; }
459429 get fps(): number { return this.params.fps; }
src/tutorial/steps.tsmodified+26−112View file
@@ -10,8 +10,8 @@
1010 // and $$ delimiters and typeset by the panel.
1111
1212 import { type Design, profileStats } from "../design.ts";
13-import { type Plate, type SimParams, bendingConstant, bendingStiffness, inPlaneSpeed, massProperties, meanThickness, nodeX, nodeY, standRates, strikePoint, thicknessRange, dispersionFrequency, CFL, listenerPosition, GRAVITY } from "../params.ts";
14-import { LN1000, plateGeometry, estimatedRate, thermoelasticRelaxation, rateText } from "../loss.ts";
13+import { type Plate, type SimParams, bendingConstant, bendingStiffness, inPlaneSpeed, massProperties, meanThickness, strikePoint, thicknessRange, dispersionFrequency, CFL, listenerPosition } from "../params.ts";
14+import { plateGeometry, estimatedRate, thermoelasticRelaxation, rateText } from "../loss.ts";
1515 import { type LinearSystem } from "../model/assemble.ts";
1616 import { maxAngularFrequency } from "../model/reference.ts";
1717 import { type ModeSet, defaultSliceSize } from "../modal/modes.ts";
@@ -129,43 +129,6 @@ function firstFreqs(modes: ModeSet | null, n: number): number[] {
129129 return modes ? Array.from(modes.freqs.slice(0, n)) : [];
130130 }
131131
132-/** How much of each of the first `count` modes is a rigid translation or a rotation of the plate: sum_k m_k phi_m(k), sum_k m_k x_k phi_m(k), sum_k m_k y_k phi_m(k), each relative to the value of the exact rigid motion. */
133-function participation(pl: Plate, model: ModalModel, count: number): { trans: number[]; rotX: number[]; rotY: number[] } {
134- const g = pl.grid;
135- const { mass, cx, cy, Ix, Iy } = massProperties(pl);
136- const n = Math.min(count, model.count);
137- const trans = new Array<number>(n).fill(0), rotX = new Array<number>(n).fill(0), rotY = new Array<number>(n).fill(0);
138- for (let r = 0; r < model.nRows; r++) {
139- const k = model.rowNode[r];
140- const m = model.masses[k];
141- const x = nodeX(g, k % g.nx) - cx, y = nodeY(g, Math.floor(k / g.nx)) - cy;
142- for (let j = 0; j < n; j++) {
143- const phi = model.wShapes[r * model.count + j];
144- trans[j] += m * phi; rotX[j] += m * y * phi; rotY[j] += m * x * phi;
145- }
146- }
147- // A rigid translation, mass-normalized, has sum m phi = sqrt(M); a rotation about x has sum m y phi = sqrt(Ix).
148- return {
149- trans: trans.map((v) => Math.abs(v) / Math.sqrt(mass)),
150- rotX: rotX.map((v) => Math.abs(v) / Math.sqrt(Ix)),
151- rotY: rotY.map((v) => Math.abs(v) / Math.sqrt(Iy)),
152- };
153-}
154-
155-/** A label for each of the first modes of a plate on a stand (rock, bounce, or plate) with the share of the mode's kinetic energy in that rigid motion. */
156-function modeKinds(pl: Plate, model: ModalModel, count: number): { label: string; share: number }[] {
157- const p = participation(pl, model, count);
158- return p.trans.map((t, m) => {
159- const r = Math.max(p.rotX[m], p.rotY[m]);
160- if (t > 0.5) return { label: "bounce", share: t * t };
161- if (r > 0.5) return { label: "rock", share: r * r };
162- return { label: "plate", share: Math.max(t, r) ** 2 };
163- });
164-}
165-function modeLabels(pl: Plate, model: ModalModel, count: number): string[] {
166- return modeKinds(pl, model, count).map((k) => k.label);
167-}
168-
169132 /** Magnitude of the radiation gain of each of the first `count` modes for the current listener, times omega^2 (a pressure per unit modal amplitude). */
170133 function pressureGains(pl: Plate, model: ModalModel, p: SimParams, count: number): number[] {
171134 const sigma = modalRates(model);
@@ -222,7 +185,6 @@ const REFERENCE_NAME = "Flat steel rectangle";
222185 /** Loads a design, evaluates it at a spacing, computes its modes, and remembers the result under a key. */
223186 async function remember(host: TutorialHost, key: string, label: string, load: () => Promise<void>, spacing: number): Promise<void> {
224187 host.setParam("spacing", spacing);
225- host.setParam("standOn", false);
226188 await load();
227189 await host.computeModes();
228190 const c = host.ctx();
@@ -239,11 +201,11 @@ export const STEPS: TutorialStep[] = [
239201 title: "What the simulator computes",
240202 body: (c) => `
241203 <p>This tutorial goes through the physics and the mathematics behind cymbalsim one idea per page, using the application itself as the laboratory. Each page reads the design in the editor and shows its numbers, and the buttons under <i>Try it</i> load presets, change parameters, compute modes and strike the plate, so that every formula on a page can be checked against the model. The pages build on one another, so a first pass is best taken in order; afterwards the list at the top jumps to any page.</p>
242-<p>The object of study is a thin metal plate, described by the script in the sidebar: its material, its outline, its thickness, an optional curved rest shape, the stand it sits on, and how fast its vibrations die away. A drumstick hits it once. The sound is computed in two ways. The <i>linear model</i> assumes the motion is small, finds the natural modes of the plate, and writes the sound as a sum of decaying sine waves, one per mode; this takes a few seconds. The <i>full model</i> keeps the terms that matter when the motion is not small and steps the equations of motion in time on the graphics processor; this takes tens of seconds per second of sound. Both use the same grid, so for a light hit they agree, which the last chapter checks.</p>
243-<p>We start with the flat steel rectangle, the simplest design: a sheet hanging freely. The stand is a setting of its own, in the Stand parameters beside the script, and the button below turns it off.</p>
244-${c.plate ? live(`In the editor: <b>${c.plate.name}</b>, ${mm(c.plate.Lx, 0)} by ${mm(c.plate.Ly, 0)}, ${c.plate.stand ? "on a stand" : "free"}, on a grid of ${mm(c.spacing, 0)}.`) : need("No design has compiled yet.")}`,
204+<p>The object of study is a thin metal plate, described by the script in the sidebar: its material, its outline, its thickness, an optional curved rest shape, and how fast its vibrations die away. A drumstick hits it once. The sound is computed in two ways. The <i>linear model</i> assumes the motion is small, finds the natural modes of the plate, and writes the sound as a sum of decaying sine waves, one per mode; this takes a few seconds. The <i>full model</i> keeps the terms that matter when the motion is not small and steps the equations of motion in time on the graphics processor; this takes tens of seconds per second of sound. Both use the same grid, so for a light hit they agree, which the last chapter checks.</p>
205+<p>We start with the flat steel rectangle, the simplest design: a sheet hanging freely. Every plate here is free: nothing holds it but the thread it hangs from, which the model leaves out.</p>
206+${c.plate ? live(`In the editor: <b>${c.plate.name}</b>, ${mm(c.plate.Lx, 0)} by ${mm(c.plate.Ly, 0)}, on a grid of ${mm(c.spacing, 0)}.`) : need("No design has compiled yet.")}`,
245207 actions: [
246- { label: "Load the flat steel rectangle, free", run: async (h) => { h.setParam("spacing", 0.008); h.setParam("standOn", false); await h.loadPreset("flat-rectangle"); } },
208+ { label: "Load the flat steel rectangle", run: async (h) => { h.setParam("spacing", 0.008); await h.loadPreset("flat-rectangle"); } },
247209 ],
248210 },
249211
@@ -389,9 +351,9 @@ ${liveText ? live(liveText) : need("Load a design to see its grid.")}
389351 },
390352 actions: [
391353 { label: "Reference at 8 mm", run: (h) => remember(h, "ref", REFERENCE_NAME, () => h.loadPreset("flat-rectangle"), 0.008) },
392- { label: "Coarser: 12 mm", run: async (h) => { h.setParam("spacing", 0.012); h.setParam("standOn", false); await h.loadPreset("flat-rectangle"); await h.computeModes(); } },
393- { label: "Finer: 6 mm", run: async (h) => { h.setParam("spacing", 0.006); h.setParam("standOn", false); await h.loadPreset("flat-rectangle"); await h.computeModes(); } },
394- { label: "Back to 8 mm", run: async (h) => { h.setParam("spacing", 0.008); h.setParam("standOn", false); await h.loadPreset("flat-rectangle"); } },
354+ { label: "Coarser: 12 mm", run: async (h) => { h.setParam("spacing", 0.012); await h.loadPreset("flat-rectangle"); await h.computeModes(); } },
355+ { label: "Finer: 6 mm", run: async (h) => { h.setParam("spacing", 0.006); await h.loadPreset("flat-rectangle"); await h.computeModes(); } },
356+ { label: "Back to 8 mm", run: async (h) => { h.setParam("spacing", 0.008); await h.loadPreset("flat-rectangle"); } },
395357 ],
396358 },
397359 {
@@ -428,9 +390,9 @@ ${liveText ? live(liveText) : need("Compute the modes of a design to see the siz
428390 const t60 = rateText;
429391 const rows: string[][] = [];
430392 for (let m = 0; m < Math.min(8, c.model.count); m++) {
431- rows.push([String(m + 1), c.model.freqs[m].toFixed(1), t60(sigma[m]), t60(c.model.sigmaRad[m]), t60(c.model.sigmaMat[m]), t60(c.model.sigmaTe[m]), t60(c.model.padSigma[m])]);
393+ rows.push([String(m + 1), c.model.freqs[m].toFixed(1), t60(sigma[m]), t60(c.model.sigmaRad[m]), t60(c.model.sigmaMat[m]), t60(c.model.sigmaTe[m])]);
432394 }
433- tab = table(["Mode", "f (Hz)", "T60 (s)", "air alone", "metal alone", "thermal alone", "felt alone"], rows);
395+ tab = table(["Mode", "f (Hz)", "T60 (s)", "air alone", "metal alone", "thermal alone"], rows);
434396 }
435397 let est = "";
436398 if (c.plate) {
@@ -454,7 +416,6 @@ $$ \\ddot q + 2\\sigma\\dot q + \\omega^2 q = 0, \\qquad q(t) = e^{-\\sigma t}\\
454416 $$ \\eta_{te} = \\Delta\\,\\frac{\\omega\\tau}{1 + (\\omega\\tau)^2}, \\qquad \\Delta = \\frac{E\\alpha^2 T}{\\rho c_p}\\,\\frac{1+\\nu}{1-\\nu}, \\qquad \\tau = \\frac{h^2 \\rho c_p}{\\pi^2 k}, $$
455417 <p>from the thermal expansion $\\alpha$, the specific heat $c_p$, the conductivity $k$ and the temperature $T$; $\\tau$ is the time heat takes to cross the thickness, and the loss peaks at $\\omega\\tau = 1$. The model applies it to the bending energy of each mode with the local thickness. For a millimeter of metal the peak lies below 100 Hz and the effect is small; it matters for thin aluminum.</p>
456418 <p><i>Radiation.</i> The sound itself carries energy away. The power a mode radiates follows from its shape through the same Rayleigh integral that gives the listener's pressure, summed over all directions (the page <i>From the plate to the ear</i>), and the mode decays at $\\sigma = P/(2E)$. A plate with free edges radiates surprisingly little below its coincidence frequency, where the bending wavelength equals the wavelength of sound: its modes have no net volume velocity, since they are orthogonal to the rigid motions, and a free edge, where the moment and the shear vanish, lacks the uncancelled strip that makes a simply supported edge radiate. The efficiency rises steeply toward coincidence, so the radiation overtakes the metal only in the top octaves, and it is there that the high partials die first. The column <i>air alone</i> below shows how weak it is for the low modes.</p>
457-<p><i>The felt of a stand</i> adds its own loss to the modes that move the center (the page on the stand).</p>
458419 ${est ? live(est) : need("Load a design to see its losses.")}
459420 <p>Try it: raise the loss factor to $10^{-3}$, the value of a heavily worked alloy, and hear the ring shorten; remove it and hear the low modes ring on with only the air to stop them.</p>`;
460421 },
@@ -462,7 +423,7 @@ ${est ? live(est) : need("Load a design to see its losses.")}
462423 { label: "Compute the modes", run: (h) => h.computeModes() },
463424 { label: "Loss factor 1e-3", run: async (h) => { await h.loadScript(rectangleScript({ eta: 1e-3 })); await h.synthesize(); } },
464425 { label: "No internal friction", run: async (h) => { await h.loadScript(rectangleScript({ eta: 0 })); await h.synthesize(); } },
465- { label: "The preset", run: async (h) => { h.setParam("standOn", false); await h.loadPreset("flat-rectangle"); await h.synthesize(); } },
426+ { label: "The preset", run: async (h) => { await h.loadPreset("flat-rectangle"); await h.synthesize(); } },
466427 { label: "Play", run: (h) => h.play() },
467428 ],
468429 },
@@ -606,54 +567,10 @@ ${live(`Window of ${N} samples at ${fs / 1000} kHz: $\\Delta f$ = ${(fs / N).toF
606567 ],
607568 },
608569
609- // ----- Stands and cymbals -----
610- {
611- id: "stand",
612- chapter: "Stands and cymbals",
613- title: "The cymbal stand",
614- body: (c) => {
615- const pl = c.plate;
616- const sr = pl ? standRates(pl) : null;
617- let liveText = "";
618- if (pl && sr && pl.stand) {
619- const pad = pl.stand.pad;
620- const np = pl.stand.topPad ? 2 : 1;
621- const { mass, Ix, Iy } = massProperties(pl);
622- const Ktheta = (np * (pad.E / pad.thickness) * Math.PI * (pad.outerRadius ** 4 - pad.innerRadius ** 4)) / 4;
623- const I = Math.min(Ix, Iy);
624- const fRock = Math.sqrt(Ktheta / I) / (2 * Math.PI);
625- const R = Math.min(pl.Lx, pl.Ly) / 2;
626- const fUmbrella = (9.08 * bendingConstant(pl.mat, meanThickness(pl))) / (2 * Math.PI * R * R);
627- liveText = `${np} pad${np > 1 ? "s" : ""} of ${mm(pad.thickness, 0)} felt at $E_p$ = ${sig(pad.E / 1e6)} MPa between ${mm(pad.innerRadius, 0)} and ${mm(pad.outerRadius, 0)}; mass $M$ = ${sig(mass * 1e3)} g.<br>$K$ = ${sig(sr.totalKw / 1e3)} kN/m, so $f_b$ = ${sr.bounceHz.toFixed(1)} Hz for a rigid cymbal; $K_\\theta$ = ${sig(Ktheta)} N m, $I$ = ${sig(I)} kg m², so $f_r$ = ${fRock.toFixed(1)} Hz. Static compression of the bottom pad under the weight alone: ${mm(mass * GRAVITY / (sr.totalKw / np), 3)}.<br>First umbrella mode of a uniform free disc of this radius and mean thickness: ${fUmbrella.toFixed(1)} Hz.`;
628- if (c.model) {
629- const kinds = modeKinds(pl, c.model, 8);
630- const sigma = modalRates(c.model);
631- const rows = kinds.map((k, m) => [String(m + 1), c.model!.freqs[m].toFixed(1), k.label, k.share.toFixed(2), (LN1000 / sigma[m]).toFixed(2)]);
632- liveText += table(["Mode", "f (Hz)", "Kind", "Rigid share", "T60 (s)"], rows);
633- } else liveText += "<br>Compute the modes to find the rock and the bounce among them.";
634- }
635- return `
636-<p>A cymbal rests on a felt pad around the rod of its stand, usually with a second pad above it under a wing nut. Felt is a soft spring: a pad of thickness $t$ and Young's modulus $E_p$ has a stiffness per unit area of $E_p/t$ in compression, and an annulus of felt between the radii $r_1$ and $r_2$ resists a vertical motion of the cymbal with the stiffness</p>
637-$$ K = n_p\\,\\frac{E_p}{t}\\,\\pi\\left(r_2^2 - r_1^2\\right), $$
638-<p>with $n_p$ the number of pads. The cymbal as a whole then bounces on the felt at $f_b = \\frac{1}{2\\pi}\\sqrt{K/M}$, and it rocks, tilting about a diameter, against the torque of the annulus, whose stiffness is its second moment $K_\\theta = n_p (E_p/t)\\,\\pi(r_2^4 - r_1^4)/4$, at $f_r = \\frac{1}{2\\pi}\\sqrt{K_\\theta/I}$ with $I$ the moment of inertia of the cymbal about a diameter. These are modes that a free plate does not have, below every plate mode: the resonances of the stand, a few hertz for the rock and tens of hertz for the bounce. The felt is lossy, and its loss factor $\\eta$ adds to each mode a decay rate proportional to the fraction of the mode's energy stored in the felt, $\\sigma = \\eta\\,\\phi^T K_s\\phi/(2\\omega)$; the rock and the bounce die within a fraction of a second, while the plate modes, which hardly move the center, are barely affected. A loss factor that is the same at every frequency is not what a spring and a dashpot give (a dashpot's loss grows with frequency), so the felt is modeled as a viscoelastic solid, a chain of springs and dashpots with relaxation times spread over the band, which the full model steps in time and the linear model evaluates at each mode's frequency; its modulus is then the value at the bounce, and it is a little stiffer above and softer below.</p>
639-<p>In the linear model the pads are linear springs. The full model lets them do what real felt does, push only while compressed, so that a hard hit can lift the cymbal off the bottom pad. The rod and its sleeve are not modeled; the shear of the felt is what holds the cymbal sideways. The stand is a setting of the setup rather than of the design: the Stand parameters beside the script put any plate on it or take it off, and set the felt's radii, thickness, modulus and loss factor, the preload of the wing nut and the position of the rod.</p>
640-${liveText ? live(liveText) : need("Turn the stand on (the Stand parameters, or the button below) to see it.")}
641-<p>The rock and the bounce are found among the computed modes by their rigid share, the fraction of a mode's kinetic energy that is a rigid translation or a rigid tilt of the cymbal: $(\\sum_k m_k\\phi_k)^2/M$ for the translation, and the same with the tilt. The rock comes out near $f_r$: at a few hertz it lies far below the first plate mode, and the cymbal tilts nearly as a whole. The bounce comes out well below $f_b$, and the reason is instructive. The formula treats the cymbal as rigid, but a thin plate is not: the pads push on its center, and, loosely speaking, the compliance of the plate around the center adds in series with that of the felt, so the center yields and the rim follows late. The bounce mode is a rigid translation blended with the umbrella-shaped plate mode, one nodal circle and no nodal diameter, whose antinode is at the center where the pads act; and the higher plate modes that move the center take a share of the felt's damping, which shortens their decay. A thicker plate is stiffer at its center and its bounce comes closer to $f_b$: compare the 20 inch ride with the flat disc.</p>`;
642- },
643- actions: [
644- { label: "The flat bronze disc on the stand", run: async (h) => { h.setParam("spacing", 0.008); h.setParam("standOn", true); await h.loadPreset("flat-disc"); } },
645- { label: "The 20 inch ride on the stand", run: async (h) => { h.setParam("spacing", 0.008); h.setParam("standOn", true); await h.loadPreset("ride-20"); } },
646- { label: "Take it off the stand", run: (h) => h.setParam("standOn", false) },
647- { label: "Compute the modes", run: (h) => h.computeModes() },
648- { label: "Show the rock", run: (h) => { const c = h.ctx(); if (c.plate && c.model) { const l = modeLabels(c.plate, c.model, 8); const m = l.indexOf("rock"); if (m >= 0) h.showMode(m); } } },
649- { label: "Show the bounce", run: (h) => { const c = h.ctx(); if (c.plate && c.model) { const l = modeLabels(c.plate, c.model, 8); const m = l.indexOf("bounce"); if (m >= 0) h.showMode(m); } } },
650- { label: "Show the first plate mode", run: (h) => { const c = h.ctx(); if (c.plate && c.model) { const l = modeLabels(c.plate, c.model, 8); const m = l.indexOf("plate"); if (m >= 0) h.showMode(m); } } },
651- { label: "Synthesize and play", run: async (h) => { await h.synthesize(); h.play(); } },
652- ],
653- },
570+ // ----- Cymbals -----
654571 {
655572 id: "taper",
656- chapter: "Stands and cymbals",
573+ chapter: "Cymbals",
657574 title: "Thick at the bell, thin at the edge",
658575 body: (c) => {
659576 const pl = c.plate;
@@ -663,11 +580,8 @@ ${liveText ? live(liveText) : need("Turn the stand on (the Stand parameters, or
663580 const R = Math.min(pl.Lx, pl.Ly) / 2;
664581 const f = (h: number) => (5.25 * bendingConstant(pl.mat, h)) / (2 * Math.PI * R * R);
665582 liveText = `Thickness from ${mm(hr.min)} to ${mm(hr.max)}; radius ${mm(R, 0)}. First mode of a uniform free disc of this radius: ${f(hr.min).toFixed(1)} Hz at the thin value, ${f(hr.max).toFixed(1)} Hz at the thick one.`;
666- if (c.model) {
667- const labels = modeLabels(pl, c.model, 8);
668- const m = labels.indexOf("plate");
669- if (m >= 0) liveText += `<br>The first plate mode computed: mode ${m + 1} at ${c.model.freqs[m].toFixed(1)} Hz${pl.stand ? " (after the stand's own modes)" : ""}.`;
670- } else liveText += "<br>Compute the modes to compare.";
583+ if (c.model && c.model.count > 0) liveText += `<br>The first mode computed: ${c.model.freqs[0].toFixed(1)} Hz.`;
584+ else liveText += "<br>Compute the modes to compare.";
671585 }
672586 return `
673587 <p>A cymbal is lathed to a taper, and the model takes the thickness as a field $h(x,y)$: the mass per area $\\rho h$ and the bending stiffness $D \\propto h^3$ vary from node to node, and the discretized energy uses the local values. Where the plate is thin, bending waves are slow and short ($\\kappa \\propto h$), so a mode that reaches the edge has its wavelength compressed there and its frequency lowered; each frequency lies between those of uniform plates of the thinnest and the thickest value, weighted toward the thickness where the mode has most of its curvature. The taper also sorts the partials in space: the thin rim, with its short waves, carries the high ones, and the thick center the low ones.</p>
@@ -676,17 +590,17 @@ ${liveText ? live(liveText) : need("Load the 16 inch crash to see its taper.")}
676590 <p>Try it: color the view by thickness to see the taper, and step through the first plate modes to see how their curvature sits on the thin rim.</p>`;
677591 },
678592 actions: [
679- { label: "Load the 16 inch crash, on the stand", run: async (h) => { h.setParam("spacing", 0.008); h.setParam("standOn", true); await h.loadPreset("crash-16"); } },
593+ { label: "Load the 16 inch crash", run: async (h) => { h.setParam("spacing", 0.008); await h.loadPreset("crash-16"); } },
680594 { label: "Color by thickness", run: (h) => h.setViewOptions({ color: "thickness" }) },
681595 { label: "Color by height", run: (h) => h.setViewOptions({ color: "height" }) },
682596 { label: "Compute the modes", run: (h) => h.computeModes() },
683- { label: "Show the first plate mode", run: (h) => { const c = h.ctx(); if (c.plate && c.model) { const l = modeLabels(c.plate, c.model, 8); const m = l.indexOf("plate"); if (m >= 0) h.showMode(m); } } },
597+ { label: "Show the first mode", run: (h) => h.showMode(0) },
684598 { label: "Synthesize and play", run: async (h) => { await h.synthesize(); h.play(); } },
685599 ],
686600 },
687601 {
688602 id: "shell",
689- chapter: "Stands and cymbals",
603+ chapter: "Cymbals",
690604 title: "Curvature: a plate becomes a shell",
691605 body: (c) => {
692606 const flat = c.memo.flat12, bowed = c.memo.bowed12;
@@ -710,7 +624,7 @@ ${liveText ? live(liveText) : need("Load a design.")}`;
710624 { label: "Show mode 1", run: (h) => h.showMode(0) },
711625 { label: "Mode 2", run: (h) => h.showMode(1) },
712626 { label: "Mode 3", run: (h) => h.showMode(2) },
713- { label: "Back to 8 mm", run: async (h) => { h.setParam("spacing", 0.008); h.setParam("standOn", false); await h.loadPreset("flat-rectangle"); } },
627+ { label: "Back to 8 mm", run: async (h) => { h.setParam("spacing", 0.008); await h.loadPreset("flat-rectangle"); } },
714628 ],
715629 },
716630
@@ -741,7 +655,7 @@ ${liveText ? live(liveText) : need("Load a design.")}
741655 <p>Try it: a hard hit on the flat disc at 12 mm spacing for 1 s of sound, which takes some tens of seconds, then the same hit in the linear model. Compare the two spectrograms: the full model's lines bend downward as the sound decays and the space between them fills in; the linear model's stay put. Then a light hit, on which the two agree.</p>`;
742656 },
743657 actions: [
744- { label: "Set up: flat disc on the stand, 12 mm, 1 s, 10 m/s", run: async (h) => { h.setParam("spacing", 0.012); h.setParam("duration", 1); h.setParam("strikeVelocity", 10); h.setParam("standOn", true); await h.loadPreset("flat-disc"); h.setRunMode("full"); } },
658+ { label: "Set up: flat disc, 12 mm, 1 s, 10 m/s", run: async (h) => { h.setParam("spacing", 0.012); h.setParam("duration", 1); h.setParam("strikeVelocity", 10); await h.loadPreset("flat-disc"); h.setRunMode("full"); } },
745659 { label: "Run the full model", needsGpu: true, run: async (h) => { h.setRunMode("full"); await h.runFull(); } },
746660 { label: "Same hit, linear model", run: async (h) => { h.setRunMode("linear"); await h.synthesize(); } },
747661 { label: "Light hit: 0.2 m/s", run: (h) => h.setParam("strikeVelocity", 0.2) },
@@ -755,13 +669,13 @@ ${liveText ? live(liveText) : need("Load a design.")}
755669 chapter: "Closing",
756670 title: "What the model leaves out",
757671 body: () => `
758-<p>Every page above named an approximation, and it is worth collecting them. The plate is thin and its rotations moderate, which holds for a cymbal moved by millimeters. A curved rest shape is treated as a shallow shell, valid for a gentle bow but strained by the deep bell of a real crash; the interface reports the largest slope. The grid places the low modes within a percent or two and flattens the high partials by the amounts on <i>The grid</i>, and its staircase edge splits the mode pairs of a disc that the hammering of a real cymbal splits anyway. The decay times are guesses that sound plausible, not measurements, and the full model matches them only through a two-parameter fit. The stand is a linear felt with a loss factor, without the rod, the sleeve or friction, and gravity enters only as the static preload. The listener hears a plate in a wall, which overstates the bass of a plate in free air, and air loading of the plate itself is neglected.</p>
672+<p>Every page above named an approximation, and it is worth collecting them. The plate is thin and its rotations moderate, which holds for a cymbal moved by millimeters. A curved rest shape is treated as a shallow shell, valid for a gentle bow but strained by the deep bell of a real crash; the interface reports the largest slope. The grid places the low modes within a percent or two and flattens the high partials by the amounts on <i>The grid</i>, and its staircase edge splits the mode pairs of a disc that the hammering of a real cymbal splits anyway. The decay times are guesses that sound plausible, not measurements, and the full model matches them only through a two-parameter fit. Nothing supports the plate: neither the thread a sheet hangs from nor the stand a cymbal sits on is modeled. The listener hears a plate in a wall, which overstates the bass of a plate in free air, and air loading of the plate itself is neglected.</p>
759673 <p>What remains is a model in which every number can be traced to a formula on these pages, and in which the two ways of computing the sound, the modes of the linearized plate and the time stepping of the nonlinear one, share one discretization and agree where they should. The document THEORY.md in the repository gives the complete equations, the discretization, and the checks that were made, and README.md describes the script and the interface.</p>
760-<p>Some experiments to try from here: an ellipse, whose degenerate mode pairs split; a thunder sheet, whose modes are dense and low; a ride with a loose wing nut (a small preload in the Stand parameters) hit hard in the full model, which lifts off its pad; a design of your own, written in the script, with a thickness or a profile the presets do not have.</p>`,
674+<p>Some experiments to try from here: an ellipse, whose degenerate mode pairs split; a thunder sheet, whose modes are dense and low; the ride hit hard in the full model; a design of your own, written in the script, with a thickness or a profile the presets do not have.</p>`,
761675 actions: [
762- { label: "Steel ellipse, free", run: async (h) => { h.setParam("spacing", 0.008); h.setParam("standOn", false); await h.loadPreset("ellipse"); } },
763- { label: "Thunder sheet, free", run: async (h) => { h.setParam("spacing", 0.012); h.setParam("standOn", false); await h.loadPreset("thunder-sheet"); } },
764- { label: "20 inch ride on the stand", run: async (h) => { h.setParam("spacing", 0.008); h.setParam("standOn", true); await h.loadPreset("ride-20"); } },
676+ { label: "Steel ellipse", run: async (h) => { h.setParam("spacing", 0.008); await h.loadPreset("ellipse"); } },
677+ { label: "Thunder sheet", run: async (h) => { h.setParam("spacing", 0.012); await h.loadPreset("thunder-sheet"); } },
678+ { label: "20 inch ride", run: async (h) => { h.setParam("spacing", 0.008); await h.loadPreset("ride-20"); } },
765679 ],
766680 },
767681 ];
src/ui/modes.tsmodified+2−3View file
@@ -11,7 +11,6 @@ export interface ModeRow {
1111 radiation: number;
1212 material: number;
1313 thermoelastic: number;
14- felt: number;
1514 }
1615
1716 const fmtT60 = t60Text;
@@ -52,7 +51,7 @@ export class ModeBrowser {
5251 this.tableWrap.className = "modes-table";
5352 const table = document.createElement("table");
5453 const thead = document.createElement("thead");
55- thead.innerHTML = "<tr><th>#</th><th>f (Hz)</th><th title=\"Decay time of the mode\">T60 (s)</th><th title=\"Decay time from radiation into the air alone\">air</th><th title=\"Decay time from the metal's internal friction alone\">metal</th><th title=\"Decay time from thermoelastic damping alone\">thermal</th><th title=\"Decay time from the felt of the stand alone\">felt</th></tr>";
54+ thead.innerHTML = "<tr><th>#</th><th>f (Hz)</th><th title=\"Decay time of the mode\">T60 (s)</th><th title=\"Decay time from radiation into the air alone\">air</th><th title=\"Decay time from the metal's internal friction alone\">metal</th><th title=\"Decay time from thermoelastic damping alone\">thermal</th></tr>";
5655 this.tbody = document.createElement("tbody");
5756 table.append(thead, this.tbody);
5857 this.tableWrap.append(table);
@@ -107,7 +106,7 @@ export class ModeBrowser {
107106 tr.dataset.m = String(m);
108107 if (m === this.sel) tr.classList.add("selected");
109108 const f = r.freq < 1000 ? r.freq.toFixed(1) : r.freq.toFixed(0);
110- tr.innerHTML = `<td>${m + 1}</td><td>${f}</td><td>${fmtT60(r.t60)}</td><td>${fmtT60(r.radiation)}</td><td>${fmtT60(r.material)}</td><td>${fmtT60(r.thermoelastic)}</td><td>${fmtT60(r.felt)}</td>`;
109+ tr.innerHTML = `<td>${m + 1}</td><td>${f}</td><td>${fmtT60(r.t60)}</td><td>${fmtT60(r.radiation)}</td><td>${fmtT60(r.material)}</td><td>${fmtT60(r.thermoelastic)}</td>`;
111110 tr.addEventListener("click", () => this.select(this.sel === m ? null : m));
112111 frag.append(tr);
113112 });
tests/assemble.test.tsmodified+5−9View file
@@ -1,8 +1,8 @@
11 // The assembled linear stiffness against the reference forces, on flat and
2-// curved plates, free and on a stand.
2+// curved plates.
33 import { test } from "node:test";
44 import assert from "node:assert/strict";
5-import { makeGrid, withMask, type Plate, nodeMasses, nodeX, nodeY, standRates } from "../src/params.ts";
5+import { makeGrid, withMask, type Plate, nodeMasses, nodeX, nodeY } from "../src/params.ts";
66 import { uniformPlate, accelerations, fields } from "../src/model/reference.ts";
77 import { assemble, stiffnessMatvec, toGrid, fromGrid, isFlatPlate } from "../src/model/assemble.ts";
88
@@ -22,7 +22,7 @@ function checkAgainstReference(pl: Plate, fieldsN: 1 | 3, seed = 3): number {
2222 const eps = 1e-9;
2323 const fx = toGrid(sys, x.map((v) => v * eps), g.n);
2424 const out = fields(g.n);
25- accelerations(pl, fx, out, false, standRates(pl));
25+ accelerations(pl, fx, out);
2626 const fm = { u: out.u.map((a, k) => (-a * masses[k]) / eps), v: out.v.map((a, k) => (-a * masses[k]) / eps), w: out.w.map((a, k) => (-a * masses[k]) / eps) };
2727 const yref = fromGrid(sys, fm);
2828 let err = 0, nrm = 0;
@@ -37,19 +37,15 @@ test("K x equals the linearized reference force on a flat rectangle (one and thr
3737 assert.ok(checkAgainstReference(pl, 3) < 1e-6);
3838 });
3939
40-test("K x equals the linearized reference force on a bowed, tapered plate with a hole on a stand", () => {
40+test("K x equals the linearized reference force on a bowed, tapered plate with a hole", () => {
4141 const g = withMask(makeGrid(0.2, 0.2, 0.0125), (x, y) => x * x + y * y <= 0.1 ** 2 && x * x + y * y >= 0.01 ** 2);
4242 const w0 = new Float64Array(g.n);
4343 for (let j = 0; j < g.ny; j++) for (let i = 0; i < g.nx; i++) w0[j * g.nx + i] = 3e-3 * Math.sin(20 * nodeX(g, i)) * Math.cos(15 * nodeY(g, j));
44- // A lossless felt: with a loss factor the reference's static pad force is the relaxed
45- // stiffness, and the rest of the linear model's stiffness at the bounce sits in the
46- // arms of the felt's chain (reference.test.ts checks those).
47- const pl: Plate = { ...uniformPlate(g, bronze, 1e-3, w0), stand: { x: 0, y: 0, pad: { innerRadius: 0.012, outerRadius: 0.03, thickness: 0.01, E: 1e6, lossFactor: 0 }, topPad: true, preload: 1e-3 } };
44+ const pl: Plate = uniformPlate(g, bronze, 1e-3, w0);
4845 for (let k = 0; k < g.n; k++) if (pl.h[k] > 0) pl.h[k] = 1e-3 * (1 + 0.4 * Math.abs(nodeX(g, k % g.nx)) / 0.1);
4946 assert.ok(!isFlatPlate(pl));
5047 const sys = assemble(pl);
5148 assert.equal(sys.fields, 3);
52- assert.ok(sys.onStand);
5349 const rel = checkAgainstReference(pl, 3);
5450 assert.ok(rel < 1e-6, `relative error ${rel}`);
5551 });
tests/design.test.tsmodified+1−30View file
@@ -2,7 +2,7 @@ import { test } from "node:test";
22 import assert from "node:assert/strict";
33 import { compileDesign, evaluatePlate, DesignError, removePlane, isFlat } from "../src/design.ts";
44 import { PRESETS } from "../src/presets.ts";
5-import { countInside, insideNode, nodeX, nodeY, nodeWeight, standRates, massProperties, standFromParams, defaultParams, type Stand } from "../src/params.ts";
5+import { countInside, insideNode, nodeX, nodeY, nodeWeight, massProperties } from "../src/params.ts";
66 import { ROOM_TEMPERATURE } from "../src/params.ts";
77
88 const base = `
@@ -99,35 +99,6 @@ test("the profile's mean plane is removed", () => {
9999 assert.ok(f.every((v) => Math.abs(v) < 1e-12));
100100 });
101101
102-test("the stand comes from the parameters, not the script, and its pads sum to the annulus stiffness", () => {
103- assert.throws(() => compileDesign(base + "design.stand.preload = 1e-3;"), /design.stand is no longer part of the script/);
104- const d = compileDesign(base + `
105-design.bounds.width = 0.3; design.bounds.height = 0.3;
106-design.inside = @(x, y, r, theta) (r <= 0.15) & (r >= 0.0065);
107-`);
108- const stand: Stand = { x: 0, y: 0, pad: { innerRadius: 7e-3, outerRadius: 20e-3, thickness: 10e-3, E: 1e6, lossFactor: 0.2 }, topPad: true, preload: 0.5e-3 };
109- const pl = evaluatePlate(d, 0.005, stand);
110- const sr = standRates(pl)!;
111- const annulus = Math.PI * (0.02 ** 2 - 0.007 ** 2);
112- const kPad = (1e6 / 0.01) * annulus;
113- assert.ok(Math.abs(sr.totalKw / (2 * kPad) - 1) < 1e-9, `total ${sr.totalKw} vs ${2 * kPad}`);
114- const { mass } = massProperties(pl);
115- // The weight rests on the relaxed felt, which is softer than the felt at the bounce.
116- assert.ok(sr.chain.kInf < 1 && sr.chain.kInf > 0);
117- assert.ok(Math.abs(sr.preloadBottom - (0.5e-3 + (mass * 9.81) / (kPad * sr.chain.kInf))) < 1e-12);
118- assert.equal(sr.preloadTop, 0.5e-3);
119- assert.ok(Math.abs(sr.bounceHz - Math.sqrt(sr.totalKw / mass) / (2 * Math.PI)) < 1e-9);
120- const pl2 = evaluatePlate(d, 0.01, { ...stand, topPad: false });
121- const sr2 = standRates(pl2)!;
122- assert.equal(sr2.preloadTop, 0);
123- assert.ok(sr2.totalKw > 0);
124- // The parameters describe the same stand; off gives a free plate.
125- const fromParams = standFromParams({ ...defaultParams, padThickness: 10e-3, padE: 1e6, padLossFactor: 0.2, preload: 0.5e-3 });
126- assert.deepEqual(fromParams, stand);
127- assert.equal(standFromParams({ ...defaultParams, standOn: false }), null);
128- assert.equal(evaluatePlate(d, 0.01).stand, null);
129-});
130-
131102 test("the loss constants of the material are read and checked", () => {
132103 const d = compileDesign(base + "design.material.lossFactor = 2e-4; design.material.alpha = 12e-6; design.material.cp = 460; design.material.k = 50; design.temperature = 300;");
133104 assert.equal(d.material.eta, 2e-4);
tests/felt.test.tsdeleted+0−65View file
@@ -1,65 +0,0 @@
1-// The viscoelastic felt: its chain of Maxwell arms has the loss factor it
2-// was built for across the audio band, the storage stiffness it was
3-// normalized to, and an exact time step.
4-import { test } from "node:test";
5-import assert from "node:assert/strict";
6-import { feltChain, feltModulus, feltStepConstants, FELT_ARMS } from "../src/felt.ts";
7-
8-test("the loss factor of the chain is flat across the band and the storage is one at the reference", () => {
9- for (const eta of [0.05, 0.3, 0.6]) {
10- for (const fRef of [15, 96, 400]) {
11- const ch = feltChain(eta, fRef);
12- assert.equal(ch.tau.length, FELT_ARMS);
13- const ref = feltModulus(ch, 2 * Math.PI * fRef);
14- assert.ok(Math.abs(ref.storage - 1) < 1e-12);
15- assert.ok(Math.abs(ref.loss / ref.storage / eta - 1) < 0.02, `eta ${eta} at the reference ${fRef} Hz: ${ref.loss / ref.storage}`);
16- let worst = 0, prev = 0;
17- for (let f = 5; f <= 10000; f *= 1.1) {
18- const m = feltModulus(ch, 2 * Math.PI * f);
19- worst = Math.max(worst, Math.abs(m.loss / m.storage / eta - 1));
20- assert.ok(m.storage > prev, "the storage stiffness rises with frequency");
21- prev = m.storage;
22- }
23- const tol = eta > 0.5 ? 0.1 : 0.06;
24- assert.ok(worst < tol, `eta ${eta}, reference ${fRef} Hz: loss factor off by up to ${(100 * worst).toFixed(1)}% between 5 Hz and 10 kHz`);
25- // Relaxed and glassy stiffness bracket the reference.
26- const glassy = ch.kInf + ch.k.reduce((a, b) => a + b, 0);
27- assert.ok(ch.kInf < 1 && glassy > 1);
28- }
29- }
30- const none = feltChain(0, 100);
31- assert.equal(none.tau.length, 0);
32- assert.equal(none.kInf, 1);
33- assert.deepEqual(feltModulus(none, 1000), { storage: 1, loss: 0 });
34-});
35-
36-test("the exact step of an arm reproduces its response to a ramp and to a sinusoid", () => {
37- const ch = feltChain(0.3, 50);
38- const a = 6; // an arm in the middle of the band
39- const tau = ch.tau[a];
40- // A ramp w = s t from a relaxed arm: p(t) = s tau (1 - exp(-t / tau)), whatever the step.
41- for (const dt of [tau / 50, tau / 3, 2 * tau]) {
42- const { decay, growth } = feltStepConstants(ch, dt);
43- const s = 1e-3;
44- let p = 0;
45- const steps = 20;
46- for (let i = 0; i < steps; i++) p = p * decay[a] + s * dt * growth[a];
47- const exact = s * tau * (1 - Math.exp(-(steps * dt) / tau));
48- assert.ok(Math.abs(p / exact - 1) < 1e-12, `dt = ${dt / tau} tau: ${p} vs ${exact}`);
49- }
50- // A sinusoid w = A sin(omega t), stepped finely: the work per cycle against
51- // the arm's force -k p is pi k A^2 omega tau / (1 + omega^2 tau^2).
52- const omega = 1 / tau, A = 1e-4, steps = 2000;
53- const dt = (2 * Math.PI) / omega / steps;
54- const { decay, growth } = feltStepConstants(ch, dt);
55- let p = 0, work = 0;
56- for (let i = 0; i < 6 * steps; i++) {
57- const w0 = A * Math.sin(omega * i * dt), w1 = A * Math.sin(omega * (i + 1) * dt);
58- const f0 = ch.k[a] * p;
59- p = p * decay[a] + (w1 - w0) * growth[a];
60- const f1 = ch.k[a] * p;
61- if (i >= 5 * steps) work += 0.5 * (f0 + f1) * (w1 - w0);
62- }
63- const expected = (Math.PI * ch.k[a] * A * A * (omega * tau)) / (1 + omega * omega * tau * tau);
64- assert.ok(Math.abs(work / expected - 1) < 1e-3, `work per cycle ${work} vs ${expected}`);
65-});
tests/gpu.test.tsmodified+3−47View file
@@ -57,10 +57,9 @@ test("the WGSL step agrees with the float64 reference on several presets", async
5757 const text = await page.evaluate(() => (window as unknown as { cymbalsim: { selftest: () => Promise<string> } }).cymbalsim.selftest());
5858 console.log(text);
5959 for (const line of text.split("\n")) {
60- const m = /max relative error ([0-9.e+-]+), felt state error ([0-9.e+-]+)/.exec(line);
60+ const m = /max relative error ([0-9.e+-]+)/.exec(line);
6161 assert.ok(m, line);
6262 assert.ok(Number(m![1]) < 2e-4, line);
63- assert.ok(Number(m![2]) < 1e-5, line);
6463 }
6564 assert.deepEqual(errors, []);
6665 });
@@ -75,7 +74,7 @@ test("the linear model computes the modes and synthesizes a strike", async () =>
7574 st = await status(page);
7675 assert.ok(st.modes && st.modes.count > 100, `modes ${st.modes?.count}`);
7776 assert.equal(st.modes!.missing, 0);
78- assert.ok(st.modes!.first[0] > 1 && st.modes!.first[0] < 200, `first mode ${st.modes!.first[0]} Hz (the rock or the bounce on the stand)`);
77+ assert.ok(st.modes!.first[0] > 1 && st.modes!.first[0] < 200, `first mode ${st.modes!.first[0]} Hz`);
7978 // A second synthesis reuses the cached modes.
8079 const t0 = Date.now();
8180 await page.evaluate(() => (window as unknown as { cymbalsim: { runLinear: () => Promise<void> } }).cymbalsim.runLinear());
@@ -230,7 +229,7 @@ test("a script with an error is reported", async () => {
230229 });
231230 });
232231
233-test("the full simulation runs a short strike on a free plate and on a stand", async () => {
232+test("the full simulation runs a short strike on a flat and on a tapered plate", async () => {
234233 await withPage("preset=bowed-rectangle&mode=full&duration=0.05&spacing=0.012", async (page, errors) => {
235234 if (!(await hasGpu(page))) { console.log("no WebGPU in this Chromium; skipping"); return; }
236235 await page.evaluate(() => (window as unknown as { cymbalsim: { runFull: () => Promise<void> } }).cymbalsim.runFull());
@@ -324,49 +323,6 @@ test("a light hit sounds the same in the linear and the full model", async () =>
324323 });
325324 });
326325
327-test("the felt damps a light hit on the aluminum square alike in the linear and the full model", async () => {
328- // The square on the default stand puts a share of many modes' energy into the felt.
329- // What the stand adds to the decay of each band, measured as the band's drop between
330- // 0.25 and 0.55 s with the stand on minus the same drop with the stand off, must agree
331- // between the models: the difference cancels the fitted plate loss of the full model
332- // against the per-mode loss of the linear one, and leaves the felt. (The first 0.25 s
333- // are left out: the modes that hold a large share of their stiffness in the felt die
334- // within it, and the linear model's first-order rates over-estimate their damping.)
335- const fs = 48000, N = 8192;
336- const bands: [number, number][] = [[50, 200], [200, 600], [600, 2000], [2000, 6000]];
337- const bandEnergy = (x: number[], t0: number, lo: number, hi: number) => {
338- const n0 = Math.round(t0 * fs);
339- const re = new Float64Array(N), im = new Float64Array(N);
340- for (let i = 0; i < N; i++) re[i] = x[n0 + i] * (0.5 - 0.5 * Math.cos((2 * Math.PI * i) / N));
341- fftInPlace(re, im);
342- let e = 0;
343- for (let k = Math.round((lo / fs) * N); k < Math.round((hi / fs) * N); k++) e += re[k] * re[k] + im[k] * im[k];
344- return e;
345- };
346- const drops = (x: number[]) => bands.map(([lo, hi]) => 10 * Math.log10(bandEnergy(x, 0.25, lo, hi) / bandEnergy(x, 0.55, lo, hi)));
347- const result: Record<string, { linear: number[]; full: number[] }> = {};
348- let skipped = false;
349- for (const stand of ["1", "0"]) {
350- await withPage(`preset=aluminum-square&duration=0.8&strikeVelocity=0.05&spacing=0.012&fps=60&standOn=${stand}`, async (page, errors) => {
351- if (!(await hasGpu(page))) { console.log("no WebGPU in this Chromium; skipping"); skipped = true; return; }
352- const grab = () => page.evaluate(() => { const a = (window as unknown as { cymbalsim: { app: { rec: { audio: Float32Array; samples: number } } } }).cymbalsim.app; return Array.from(a.rec.audio.subarray(0, a.rec.samples)); });
353- await page.evaluate(() => (window as unknown as { cymbalsim: { runLinear: () => Promise<void> } }).cymbalsim.runLinear());
354- const linear = await grab();
355- await page.evaluate(() => (window as unknown as { cymbalsim: { runFull: () => Promise<void> } }).cymbalsim.runFull());
356- const full = await grab();
357- assert.equal(linear.length, full.length);
358- result[stand] = { linear: drops(linear), full: drops(full) };
359- assert.deepEqual(errors, []);
360- });
361- if (skipped) return;
362- }
363- bands.forEach(([lo, hi], b) => {
364- const dl = result["1"].linear[b] - result["0"].linear[b], df = result["1"].full[b] - result["0"].full[b];
365- console.log(`${lo} to ${hi} Hz: decay over 0.3 s added by the stand, linear ${dl.toFixed(1)} dB (${result["0"].linear[b].toFixed(1)} free), full ${df.toFixed(1)} dB (${result["0"].full[b].toFixed(1)} free)`);
366- assert.ok(Math.abs(dl - df) < 4, `${lo} to ${hi} Hz: the stand adds ${dl} dB in the linear model and ${df} dB in the full model`);
367- });
368-});
369-
370326 function fftInPlace(re: Float64Array, im: Float64Array): void {
371327 const n = re.length;
372328 for (let i = 1, j = 0; i < n; i++) {
tests/loss.test.tsmodified+1−1View file
@@ -135,7 +135,7 @@ test("simply supported shapes reproduce Wallace's low-frequency radiation effici
135135 const v = Math.sin((m * Math.PI * i * g.hx) / a) * Math.sin((n * Math.PI * j * g.hy) / b);
136136 rows.push(kk); vals.push(v); norm += rhoH * g.hx * g.hy * w * v * v;
137137 }
138- const one = { nGrid: g.n, nRows: rows.length, count: 1, freqs: new Float64Array([f]), wShapes: Float32Array.from(vals.map((v) => v / Math.sqrt(norm))), wRow: new Int32Array(0), rowNode: Int32Array.from(rows), sigmaMat: new Float64Array(1), sigmaTe: new Float64Array(1), sigmaRad: new Float64Array(1), padSigma: new Float64Array(1), masses: new Float64Array(0) };
138+ const one = { nGrid: g.n, nRows: rows.length, count: 1, freqs: new Float64Array([f]), wShapes: Float32Array.from(vals.map((v) => v / Math.sqrt(norm))), wRow: new Int32Array(0), rowNode: Int32Array.from(rows), sigmaMat: new Float64Array(1), sigmaTe: new Float64Array(1), sigmaRad: new Float64Array(1), masses: new Float64Array(0) };
139139 const efficiency = (radiationRateDirect(one, dp, 0) * rhoH) / (RHO_AIR * c); // sigma = rho_air c sigma_bar / (rho h)
140140 const odd = (q: number) => q % 2 === 1;
141141 const wallace = odd(m) && odd(n) ? (32 * k * k * a * b) / (Math.PI ** 5 * m * m * n * n)
tests/modal.test.tsmodified+0−14View file
@@ -164,17 +164,3 @@ test("free circular plate frequencies are close to Leissa's values", async () =>
164164 }
165165 });
166166
167-test("a plate on a stand has no rigid modes and bounces at the pad frequency", async () => {
168- const g = withMask(makeGrid(0.2, 0.2, 0.01), (x, y) => x * x + y * y <= 0.1 ** 2 && x * x + y * y >= 0.01 ** 2);
169- const pl: Plate = { ...uniformPlate(g, steel, 2e-3), stand: { x: 0, y: 0, pad: { innerRadius: 0.012, outerRadius: 0.03, thickness: 0.01, E: 0.3e6, lossFactor: 0.3 }, topPad: true, preload: 1e-3 } };
170- const sys = assemble(pl);
171- const ms = await computeModes(sys, 400, localRunner(sys), { sliceSize: 20 });
172- const { standRates, massProperties } = await import("../src/params.ts");
173- const sr = standRates(pl)!;
174- const { mass } = massProperties(pl);
175- const bounce = Math.sqrt(sr.totalKw / mass) / (2 * Math.PI);
176- // The lowest modes are the two rocking modes and the bounce; the bounce is close to sqrt(K/M).
177- assert.ok(ms.freqs[0] > 1, "no zero-frequency mode");
178- const near = Array.from(ms.freqs).some((f) => Math.abs(f / bounce - 1) < 0.15);
179- assert.ok(near, `bounce ${bounce} Hz among ${Array.from(ms.freqs.slice(0, 5)).map((f) => f.toFixed(1)).join(", ")}`);
180-});
tests/reference.test.tsmodified+6−92View file
@@ -1,13 +1,12 @@
11 // The float64 reference discretization: symmetry and rigid modes of the
22 // linear stiffness, forces as the gradient of the energy, dissipativity,
3-// thickness scaling, the stand, and the free-plate frequencies.
3+// thickness scaling, and the free-plate frequencies.
44 import { test } from "node:test";
55 import assert from "node:assert/strict";
6-import { makeGrid, withMask, type Grid, type Plate, nodeMasses, standRates, insideNode, nodeX, nodeY, bendingStiffness } from "../src/params.ts";
7-import { uniformPlate, accelerations, addDamping, addPadMemory, padState, stepPadMemory, potentialEnergy, fields, maxAngularFrequency, type Fields } from "../src/model/reference.ts";
6+import { makeGrid, withMask, type Grid, type Plate, nodeMasses, insideNode, nodeX, nodeY, bendingStiffness } from "../src/params.ts";
7+import { uniformPlate, accelerations, addDamping, potentialEnergy, fields, maxAngularFrequency, type Fields } from "../src/model/reference.ts";
88 import { assemble, stiffnessEntry } from "../src/model/assemble.ts";
99 import { fullLoss } from "../src/model/reference.ts";
10-import { feltModulus } from "../src/felt.ts";
1110
1211 const steel = { E: 200e9, nu: 0.3, rho: 7850, eta: 0, thermal: null };
1312
@@ -30,10 +29,6 @@ function bowed(g: Grid, sag: number): Float64Array {
3029 for (let j = 0; j < g.ny; j++) for (let i = 0; i < g.nx; i++) w0[j * g.nx + i] = sag * (1 - (2 * nodeX(g, i) / Lx) ** 2) + 0.3 * sag * Math.sin(3 * nodeY(g, j) / Lx);
3130 return w0;
3231 }
33-function withStand(pl: Plate): Plate {
34- return { ...pl, stand: { x: 0, y: 0, pad: { innerRadius: 0.005, outerRadius: 0.03, thickness: 0.01, E: 1e6, lossFactor: 0.3 }, topPad: true, preload: 1e-3 } };
35-}
36-
3732 test("the linear stiffness is symmetric with exactly three rigid modes per block (uniform, tapered, disc with hole)", () => {
3833 const cases: { name: string; pl: Plate; fields: 1 | 3 }[] = [];
3934 const g1 = makeGrid(0.12, 0.09, 0.015);
@@ -85,7 +80,7 @@ test("the linear stiffness is symmetric with exactly three rigid modes per block
8580 import * as assembleModule from "../src/model/assemble.ts";
8681 function require_assemble(): typeof assembleModule { return assembleModule; }
8782
88-test("forces are the exact gradient of the discrete energy (flat and bowed, tapered, free and on a stand)", () => {
83+test("forces are the exact gradient of the discrete energy (flat and bowed, tapered)", () => {
8984 const g = makeGrid(0.1, 0.08, 0.01);
9085 const variants: { name: string; pl: Plate }[] = [];
9186 variants.push({ name: "flat free", pl: uniformPlate(g, steel, 1e-3) });
@@ -93,19 +88,17 @@ test("forces are the exact gradient of the discrete energy (flat and bowed, tape
9388 const tapered = uniformPlate(g, steel, 1e-3, bowed(g, 2e-3));
9489 for (let k = 0; k < g.n; k++) if (tapered.h[k] > 0) tapered.h[k] = 1e-3 * (1 + 0.4 * (k % g.nx) / g.nx);
9590 variants.push({ name: "bowed tapered", pl: tapered });
96- variants.push({ name: "bowed on a stand", pl: withStand(uniformPlate(g, steel, 1e-3, bowed(g, 3e-3))) });
9791 for (const v of variants) {
9892 const pl = v.pl;
99- const sr = standRates(pl);
10093 const f = randomFields(pl, 2.5e-3, 7);
10194 const masses = nodeMasses(pl);
10295 const acc = fields(g.n);
103- accelerations(pl, f, acc, false, sr);
96+ accelerations(pl, f, acc);
10497 const dir = randomFields(pl, 1, 11);
10598 // Directional derivative of U along dir versus -F . dir.
10699 const eps = 1e-7;
107100 const shift = (s: number): Fields => ({ u: f.u.map((x, k) => x + s * dir.u[k]), v: f.v.map((x, k) => x + s * dir.v[k]), w: f.w.map((x, k) => x + s * dir.w[k]) });
108- const dU = (potentialEnergy(pl, shift(eps), sr) - potentialEnergy(pl, shift(-eps), sr)) / (2 * eps);
101+ const dU = (potentialEnergy(pl, shift(eps)) - potentialEnergy(pl, shift(-eps))) / (2 * eps);
109102 let fd = 0;
110103 for (let k = 0; k < g.n; k++) fd -= masses[k] * (acc.u[k] * dir.u[k] + acc.v[k] * dir.v[k] + acc.w[k] * dir.w[k]);
111104 assert.ok(Math.abs(dU - fd) < 1e-6 * Math.abs(dU), `${v.name}: dU ${dU} vs -F.dir ${fd}`);
@@ -127,47 +120,6 @@ test("the loss terms only remove energy", () => {
127120 }
128121 });
129122
130-test("the felt's chain, stepped in time, dissipates the work of its loss stiffness per cycle", () => {
131- const g = withMask(makeGrid(0.2, 0.2, 0.01), (x, y) => x * x + y * y <= 0.1 ** 2 && x * x + y * y >= 0.01 ** 2);
132- const pl = withStand(uniformPlate(g, steel, 1e-3));
133- const sr = standRates(pl)!;
134- const masses = nodeMasses(pl);
135- // A rigid translation w = A sin(omega t) of the plate, stepped finely; the work done
136- // against the arms over a cycle, once the transient of the arms has passed, is
137- // pi L(omega) K A^2 with L the loss stiffness per unit of K (eta at the bounce).
138- for (const f of [sr.bounceHz, 2000]) {
139- const omega = 2 * Math.PI * f, A = 1e-4;
140- const steps = 400, dt = 1 / (f * steps);
141- const pad = padState(pl, sr);
142- const x = fields(g.n), dw = fields(g.n), acc = fields(g.n);
143- const setW = (t: number, out: Fields) => { for (let k = 0; k < g.n; k++) out.w[k] = insideNode(g, k % g.nx, Math.floor(k / g.nx)) ? A * Math.sin(omega * t) : 0; };
144- let work = 0;
145- const cycles = 4;
146- for (let s = 0; s < cycles * steps; s++) {
147- const t0 = s * dt, t1 = (s + 1) * dt;
148- setW(t0, x); setW(t1, dw);
149- for (let k = 0; k < g.n; k++) dw.w[k] -= x.w[k];
150- acc.u.fill(0); acc.v.fill(0); acc.w.fill(0);
151- addPadMemory(pl, sr, x, pad, false, acc);
152- // Work of the plate against the arms over the step, with the force at mid-step (trapezoid).
153- let F0 = 0;
154- for (let k = 0; k < g.n; k++) F0 += masses[k] * acc.w[k];
155- stepPadMemory(pl, sr, dt, dw, x, pad, false);
156- setW(t1, x);
157- acc.u.fill(0); acc.v.fill(0); acc.w.fill(0);
158- addPadMemory(pl, sr, x, pad, false, acc);
159- let F1 = 0;
160- for (let k = 0; k < g.n; k++) F1 += masses[k] * acc.w[k];
161- if (s >= (cycles - 1) * steps) work -= 0.5 * (F0 + F1) * dw.w[sr.nodes[0]];
162- }
163- const L = feltModulus(sr.chain, omega).loss;
164- const expected = Math.PI * L * sr.totalKw * A * A;
165- assert.ok(work > 0);
166- assert.ok(Math.abs(work / expected - 1) < 0.02, `at ${f.toFixed(0)} Hz: work per cycle ${work} vs ${expected}`);
167- if (f === sr.bounceHz) assert.ok(Math.abs(L / sr.eta - 1) < 0.02, `loss stiffness at the bounce ${L} vs eta ${sr.eta}`);
168- }
169-});
170-
171123 test("doubling the thickness doubles the bending frequencies and leaves the membrane frequencies alone", () => {
172124 const g = makeGrid(0.1, 0.1, 0.01);
173125 const w1 = maxAngularFrequency(g, steel, 1e-3), w2 = maxAngularFrequency(g, steel, 2e-3);
@@ -182,41 +134,3 @@ test("doubling the thickness doubles the bending frequencies and leaves the memb
182134 assert.ok(Math.abs(bendingStiffness(steel, 2e-3) / bendingStiffness(steel, 1e-3) - 8) < 1e-12);
183135 });
184136
185-test("the stand resists a rigid translation with its total stiffness and a tilt with the pads' second moment", () => {
186- const g = withMask(makeGrid(0.2, 0.2, 0.01), (x, y) => x * x + y * y <= 0.1 ** 2 && x * x + y * y >= 0.01 ** 2);
187- const pl = withStand(uniformPlate(g, steel, 1e-3));
188- const sr = standRates(pl)!;
189- const masses = nodeMasses(pl);
190- const f = fields(g.n), acc = fields(g.n);
191- const dz = 1e-4;
192- for (let k = 0; k < g.n; k++) f.w[k] = insideNode(g, k % g.nx, Math.floor(k / g.nx)) ? dz : 0;
193- accelerations(pl, f, acc, false, sr);
194- let F = 0;
195- for (let k = 0; k < g.n; k++) F += masses[k] * acc.w[k];
196- // The relaxed felt holds kInf of the stiffness at the bounce; the arms of the chain hold the rest.
197- const kInf = sr.chain.kInf;
198- assert.ok(kInf > 0.2 && kInf < 0.8, `relaxed fraction ${kInf}`);
199- assert.ok(Math.abs(-F / dz / (kInf * sr.totalKw) - 1) < 1e-9, `translation stiffness ${-F / dz} vs ${kInf * sr.totalKw}`);
200- // Tilt about y: w = theta x; restoring moment = theta * sum_k kw_k x_k^2 (bending of a plane is
201- // zero; the von Karman membrane strain of a rotation is theta^2 / 2, so a tiny tilt keeps it negligible).
202- const th = 1e-7;
203- for (let k = 0; k < g.n; k++) f.w[k] = insideNode(g, k % g.nx, Math.floor(k / g.nx)) ? th * nodeX(g, k % g.nx) : 0;
204- accelerations(pl, f, acc, false, sr);
205- let M = 0, I2 = 0;
206- for (let k = 0; k < g.n; k++) { M += masses[k] * acc.w[k] * nodeX(g, k % g.nx); I2 += kInf * sr.kw[k] * nodeX(g, k % g.nx) ** 2; }
207- assert.ok(Math.abs(-M / th / I2 - 1) < 1e-6, `tilt stiffness ${-M / th} vs ${I2}`);
208- // With lift-off, a small motion within the preload is the same linear spring; a large upward motion unloads the bottom pad.
209- for (let k = 0; k < g.n; k++) f.w[k] = insideNode(g, k % g.nx, Math.floor(k / g.nx)) ? 1e-5 : 0;
210- const lin = fields(g.n), lo = fields(g.n);
211- accelerations(pl, f, lin, false, sr);
212- accelerations(pl, f, lo, true, sr);
213- for (let k = 0; k < g.n; k++) assert.ok(Math.abs(lin.w[k] - lo.w[k]) < 1e-9 * Math.max(1, Math.abs(lin.w[k])));
214- for (let k = 0; k < g.n; k++) f.w[k] = insideNode(g, k % g.nx, Math.floor(k / g.nx)) ? 0.01 : 0;
215- accelerations(pl, f, lo, true, sr);
216- let Flo = 0;
217- for (let k = 0; k < g.n; k++) Flo += masses[k] * lo.w[k];
218- // Bottom pad unloaded: only the top pad (compressed by preload + 10 mm) pushes down, minus the rest balance.
219- const kb = (kInf * sr.totalKw) / 2;
220- const expected = -(kb * (sr.preloadTop + 0.01) - kb * sr.preloadTop) - kb * sr.preloadBottom;
221- assert.ok(Math.abs(Flo / expected - 1) < 1e-6, `lift-off force ${Flo} vs ${expected}`);
222-});
tests/synth.test.tsmodified+2−2View file
@@ -29,7 +29,7 @@ test("a single mode under an impulse follows the closed-form damped sinusoid", a
2929 // Keep only mode 0 by slicing the model.
3030 const one = sliceModel(model, 0, 1);
3131 // Give the mode a decay time of exactly 1 s through the material term alone.
32- one.sigmaMat[0] = LN1000; one.sigmaTe[0] = 0; one.sigmaRad[0] = 0; one.padSigma[0] = 0;
32+ one.sigmaMat[0] = LN1000; one.sigmaTe[0] = 0; one.sigmaRad[0] = 0;
3333 const sigma = modalRates(one);
3434 const g = strikeGains(one, pl, p);
3535 const c = simulateContact(one, sigma, g, p, p.fps);
@@ -54,7 +54,7 @@ test("the radiation gain of a uniform piston equals the analytic Rayleigh sum",
5454 const { model, pl } = await smallModel();
5555 const p: SimParams = { ...defaultParams, listenerDist: 1.0, listenerAngle: 0 };
5656 // A fake model with one mode of uniform shape (mass-normalized to 1/sqrt(M)).
57- const uniform: ModalModel = { ...sliceModel(model, 0, 1), freqs: new Float64Array([1000]), padSigma: new Float64Array([0]) };
57+ const uniform: ModalModel = { ...sliceModel(model, 0, 1), freqs: new Float64Array([1000]) };
5858 let M = 0;
5959 for (let k = 0; k < uniform.nGrid; k++) M += uniform.masses[k];
6060 const phi = 1 / Math.sqrt(M);
tests/tutorial.test.tsmodified+3−13View file
@@ -8,7 +8,7 @@ import { test } from "node:test";
88 import assert from "node:assert/strict";
99 import { compileDesign, evaluatePlate } from "../src/design.ts";
1010 import { presetById } from "../src/presets.ts";
11-import { defaultParams, standFromParams, type SimParams } from "../src/params.ts";
11+import { defaultParams, type SimParams } from "../src/params.ts";
1212 import { assemble } from "../src/model/assemble.ts";
1313 import { computeModes, localRunner } from "../src/modal/modes.ts";
1414 import { modalModel, synthesize } from "../src/modal/synth.ts";
@@ -31,7 +31,7 @@ class TestHost implements TutorialHost {
3131 played = 0;
3232 private async build(source: string): Promise<void> {
3333 const d = compileDesign(source);
34- const pl = evaluatePlate(d, this.params.spacing, standFromParams(this.params));
34+ const pl = evaluatePlate(d, this.params.spacing);
3535 this.c = { ...this.c, design: d, plate: pl, params: { ...this.params }, spacing: this.params.spacing, sys: null, modes: null, model: null, lastContact: null };
3636 }
3737 ctx(): TutorialContext { return { ...this.c, params: { ...this.params }, runMode: this.runMode }; }
@@ -39,7 +39,7 @@ class TestHost implements TutorialHost {
3939 async loadScript(source: string): Promise<void> { await this.build(source); }
4040 setParam<K extends keyof SimParams>(key: K, value: SimParams[K]): void {
4141 this.params[key] = value;
42- if ((key === "spacing" || String(key).startsWith("stand") || String(key).startsWith("pad") || key === "topPad" || key === "preload") && this.c.design) void this.build(this.c.design.source);
42+ if (key === "spacing" && this.c.design) void this.build(this.c.design.source);
4343 }
4444 async computeModes(): Promise<void> {
4545 const pl = this.c.plate;
@@ -85,16 +85,6 @@ test("every step renders with an empty context, with a plate, and with its modes
8585 const c = host.ctx();
8686 assert.ok(c.modes && c.modes.count > 10);
8787 for (const s of STEPS) check(s.body(c), `${s.id} (modes)`);
88-
89- // The stand pages with a plate on a stand.
90- host.setParam("standOn", true);
91- await host.loadPreset("flat-disc");
92- await host.computeModes();
93- const cs = host.ctx();
94- const stand = STEPS[stepIndex("stand")].body(cs);
95- check(stand, "stand (disc)");
96- assert.match(stand, /bounce/);
97- assert.match(stand, /rock/);
9888 });
9989
10090 test("the actions run against an in-process host", async () => {