diff --git a/pop/bell/bin/viz.mjs b/pop/bell/bin/viz.mjs new file mode 100755 index 000000000..f8d0aca2e --- /dev/null +++ b/pop/bell/bin/viz.mjs @@ -0,0 +1,336 @@ +#!/usr/bin/env node +// viz.mjs — 3-D visualization of the bell physical model. +// +// Reads the modes JSON exported by the C engine (`bell ... --modes out.json`), +// the SAME data that produced the audio, and animates the bell shell deforming +// as the live sum of its excited modes decays. You literally watch the (m,n) +// rim-flexural modes ring and fade. A small software 3-D pipeline (no three.js, +// matching the repo's no-deps ethos) projects the surface of revolution; frames +// are piped as BGRA to ffmpeg via pop/lib/preview-shared.mjs and muxed with the +// rendered bell audio. +// +// Usage: +// node viz.mjs --modes bell-modes.json --audio bell.wav --out bell-viz.mp4 +// [--dur 9] [--fps 30] [--portrait] [--w 1920 --h 1080] +// +// The audio's first strike is assumed at t=0 (as bell --out renders it). + +import { createCanvas } from "canvas"; +import { resolve } from "node:path"; +import { readFileSync } from "node:fs"; +import { spawn, spawnSync } from "node:child_process"; +import { + spawnFFmpegEncode, + decodeAudioMono, + computeRmsEnvelope, +} from "../../lib/preview-shared.mjs"; + +// ---- args ------------------------------------------------------------------ +const args = process.argv.slice(2); +function arg(name, def) { + const i = args.indexOf(name); + return i >= 0 && i + 1 < args.length ? args[i + 1] : def; +} +const modesPath = resolve(arg("--modes", "/tmp/bell-modes.json")); +const audioPath = arg("--audio", null); +const outPath = resolve(arg("--out", "/tmp/bell-viz.mp4")); +const fps = parseInt(arg("--fps", "30"), 10); +const dur = parseFloat(arg("--dur", "9")); +const strikeAt = parseFloat(arg("--strike", "0.7")); // seconds of rest first +const portrait = args.includes("--portrait"); +const W = parseInt(arg("--w", portrait ? "1080" : "1920"), 10); +const H = parseInt(arg("--h", portrait ? "1920" : "1080"), 10); + +const model = JSON.parse(readFileSync(modesPath, "utf8")); +const geo = model.geometry; +const modes = model.modes; + +// ---- metal palette by material -------------------------------------------- +const METAL = { + bronze: [176, 141, 58], brass: [201, 162, 75], steel: [154, 163, 173], + aluminum: [200, 204, 208], silver: [216, 221, 227], glass: [127, 212, 224], + gold: [227, 194, 0], +}; +const metal = METAL[model.material.name] || [180, 150, 90]; + +// =========================================================================== +// Build the surface-of-revolution mesh from the meridian profile. +// =========================================================================== +const NS = geo.n; // meridian stations +const NT = portrait ? 80 : 96; // angular divisions + +// Center the bell vertically and find a scale that fills the frame. +let zmin = Infinity, zmax = -Infinity, rmax = 0; +for (let i = 0; i < NS; i++) { + zmin = Math.min(zmin, geo.z[i]); + zmax = Math.max(zmax, geo.z[i]); + rmax = Math.max(rmax, geo.r[i]); +} +const zc = 0.5 * (zmin + zmax); +const span = Math.max(zmax - zmin, 2 * rmax); + +// 2-D meridian tangent/normal per node (in the (r,z) plane). +const nr = new Float64Array(NS), nz = new Float64Array(NS); +const tr = new Float64Array(NS), tz = new Float64Array(NS); +for (let i = 0; i < NS; i++) { + const a = Math.max(0, i - 1), b = Math.min(NS - 1, i + 1); + let dr = geo.r[b] - geo.r[a], dz = geo.z[b] - geo.z[a]; + const L = Math.hypot(dr, dz) || 1; + dr /= L; dz /= L; + tr[i] = dr; tz[i] = dz; + nr[i] = dz; nz[i] = -dr; // rotate tangent -> outward-ish normal +} + +// Precompute angle tables and cos(m*theta) per mode. +const cosT = new Float64Array(NT), sinT = new Float64Array(NT); +for (let j = 0; j < NT; j++) { + const th = (j / NT) * Math.PI * 2; + cosT[j] = Math.cos(th); sinT[j] = Math.sin(th); +} +const cosMT = modes.map((md) => { + const a = new Float64Array(NT); + for (let j = 0; j < NT; j++) a[j] = Math.cos(md.m * (j / NT) * Math.PI * 2); + return a; +}); + +// Visual deformation gain: largest excited mode reaches ~12% of the bell span. +let maxPart = 1e-6; +for (const md of modes) maxPart = Math.max(maxPart, Math.abs(md.part)); +const GAIN = 0.12 * span / maxPart; + +// =========================================================================== +// 3-D helpers (right-handed; Y up). Camera looks down -Z after we push +Z. +// =========================================================================== +const CAM_D = span * 2.5; // camera distance +const FOC = Math.min(W, H) * 1.15; // focal length (px) +const TILT = 0.32; // slight downward look (radians) +const cosTilt = Math.cos(TILT), sinTilt = Math.sin(TILT); + +function project(x, y, z, yaw) { + // rotate about Y (yaw) + const cx = Math.cos(yaw), sx = Math.sin(yaw); + let X = cx * x + sx * z; + let Z = -sx * x + cx * z; + let Y = y; + // tilt about X + const Y2 = cosTilt * Y - sinTilt * Z; + const Z2 = sinTilt * Y + cosTilt * Z; + // push away from camera + const Zc = Z2 + CAM_D; + const inv = FOC / Math.max(Zc, 1e-3); + return [W / 2 + X * inv, H / 2 - Y2 * inv, Zc]; +} + +// =========================================================================== +// Audio-reactive strike flash (optional). +// =========================================================================== +let rms = null, rmsFps = 60; +if (audioPath) { + try { + const { audio, sr } = decodeAudioMono(audioPath); + rms = computeRmsEnvelope(audio, sr, rmsFps, dur); + } catch (e) { + console.warn("viz: audio envelope unavailable:", e.message); + } +} + +// =========================================================================== +// Render loop. +// =========================================================================== +const canvas = createCanvas(W, H); +const ctx = canvas.getContext("2d"); +const totalFrames = Math.round(dur * fps); + +// Delay the audio so its strike lands at strikeAt (the bell rests first). +let muxAudio = audioPath; +if (audioPath && strikeAt > 0) { + const ms = Math.round(strikeAt * 1000); + const padded = `${outPath}.aud.wav`; + const p = spawnSync("ffmpeg", [ + "-hide_banner", "-loglevel", "error", "-y", "-i", audioPath, + "-af", `adelay=${ms}|${ms}`, padded, + ], { stdio: "inherit" }); + if (p.status === 0) muxAudio = padded; +} + +// With audio: use the shared encoder (muxes audio). Without: video-only. +const encoder = muxAudio + ? spawnFFmpegEncode({ audioPath: muxAudio, w: W, h: H, fps, outPath, crf: 18 }) + : spawn("ffmpeg", [ + "-hide_banner", "-loglevel", "error", "-y", + "-f", "rawvideo", "-pix_fmt", "bgra", "-s", `${W}x${H}`, "-r", String(fps), + "-i", "-", "-c:v", "libx264", "-preset", "faster", "-crf", "18", + "-pix_fmt", "yuv420p", "-movflags", "+faststart", outPath, + ], { stdio: ["pipe", "inherit", "inherit"] }); + +// Reusable per-vertex scratch. +const PX = new Float64Array(NS * NT); +const PY = new Float64Array(NS * NT); +const PZ = new Float64Array(NS * NT); +const DN = new Float64Array(NS * NT); // signed normal displacement (for tint) + +const strikeFreq = model.strike_freq; + +function modalCoords(tm) { + // tm = time since the strike. Before the strike the bell sits at rest; after, + // q_k = part * exp(-tm/tau) * cos(2*pi*f*tm) — rings out and settles back. + const q = new Float64Array(modes.length); + if (tm < 0) return q; // rest + for (let k = 0; k < modes.length; k++) { + const md = modes[k]; + q[k] = md.part * Math.exp(-tm / md.tau) * Math.cos(2 * Math.PI * md.freq * tm); + } + return q; +} + +function drawHUD(tm) { + ctx.save(); + ctx.font = "600 30px monospace"; + ctx.fillStyle = "rgba(235,238,245,0.92)"; + ctx.textBaseline = "top"; + ctx.fillText("BELL · physical model", 40, 36); + ctx.font = "400 22px monospace"; + ctx.fillStyle = "rgba(190,198,210,0.85)"; + ctx.fillText(`${geo.name} · ${model.material.name}`, 40, 78); + ctx.fillText(`strike ${strikeFreq.toFixed(1)} Hz`, 40, 106); + ctx.fillText(`E ${(model.material.E / 1e9).toFixed(0)} GPa ρ ${model.material.rho.toFixed(0)} ν ${model.material.nu}`, 40, 134); + // state label: at rest before the strike, ringing after. + if (tm < 0) { + ctx.fillStyle = "rgba(120,130,145,0.8)"; + ctx.fillText("· at rest ·", 40, 168); + } else { + ctx.fillStyle = "rgba(235,210,150,0.9)"; + ctx.fillText(`· ringing ${tm.toFixed(2)}s ·`, 40, 168); + } + ctx.textBaseline = "alphabetic"; + + // Spectrum bars: live decaying amplitude per partial, hued by m. + const N = Math.min(modes.length, 24); + const x0 = 40, y0 = H - 56, bw = 26, gap = 6, bh = 150; + for (let k = 0; k < N; k++) { + const md = modes[k]; + const env = tm < 0 ? 0 : md.part * Math.exp(-tm / md.tau); + const hh = Math.min(1, env) * bh; + const hue = (md.m - 2) * 48; // m=2 red-ish, climbing + ctx.fillStyle = `hsla(${hue},70%,55%,0.85)`; + ctx.fillRect(x0 + k * (bw + gap), y0 - hh, bw, hh); + ctx.fillStyle = "rgba(120,128,140,0.6)"; + ctx.font = "400 13px monospace"; + ctx.fillText(String(md.m), x0 + k * (bw + gap) + 8, y0 + 16); + } + ctx.restore(); +} + +async function writeFrame(buf) { + if (!encoder.stdin.write(buf)) + await new Promise((r) => encoder.stdin.once("drain", r)); +} + +async function main() { + for (let f = 0; f < totalFrames; f++) { + const t = f / fps; + const tm = t - strikeAt; // time since strike (<0 = at rest) + const yaw = t * 0.5; // slow reveal + const q = modalCoords(tm); + + // Deform + project every vertex. + for (let i = 0; i < NS; i++) { + const ri = geo.r[i], zi = geo.z[i] - zc; + for (let j = 0; j < NT; j++) { + let dN = 0, dM = 0; + for (let k = 0; k < modes.length; k++) { + const c = GAIN * q[k] * cosMT[k][j]; + dN += c * modes[k].w[i]; + dM += c * modes[k].u[i]; + } + const idx = i * NT + j; + DN[idx] = dN; + // base point + normal/tangent displacement, expanded to 3-D by theta. + const rr = ri + nr[i] * dN + tr[i] * dM; + const zz = zi + nz[i] * dN + tz[i] * dM; + const x = rr * cosT[j]; + const z = rr * sinT[j]; + const [sx, sy, depth] = project(x, zz, z, yaw); + PX[idx] = sx; PY[idx] = sy; PZ[idx] = depth; + } + } + + // Background. + const bg = ctx.createLinearGradient(0, 0, 0, H); + bg.addColorStop(0, "#0b0e13"); + bg.addColorStop(1, "#05070a"); + ctx.fillStyle = bg; + ctx.fillRect(0, 0, W, H); + + // strike flash vignette (indexed by time-since-strike) + if (rms && tm >= 0) { + const fi = Math.min(rms.length - 1, Math.round(tm * rmsFps)); + const a = Math.min(0.25, (rms[fi] || 0) * 0.5); + if (a > 0.01) { + const g = ctx.createRadialGradient(W / 2, H / 2, 0, W / 2, H / 2, Math.max(W, H) * 0.6); + g.addColorStop(0, `rgba(${metal[0]},${metal[1]},${metal[2]},${a})`); + g.addColorStop(1, "rgba(0,0,0,0)"); + ctx.fillStyle = g; + ctx.fillRect(0, 0, W, H); + } + } + + // Build quads, sort back-to-front (painter's), shade + draw. + const quads = []; + for (let i = 0; i < NS - 1; i++) { + for (let j = 0; j < NT; j++) { + const j2 = (j + 1) % NT; + const a = i * NT + j, b = i * NT + j2; + const c = (i + 1) * NT + j2, d = (i + 1) * NT + j; + const depth = (PZ[a] + PZ[b] + PZ[c] + PZ[d]) * 0.25; + const disp = (DN[a] + DN[b] + DN[c] + DN[d]) * 0.25; + quads.push([a, b, c, d, depth, disp]); + } + } + quads.sort((p, q2) => q2[4] - p[4]); // far first + + const Ly = 0.72; // dominant light component (overhead-ish) + for (const [a, b, c, d, , disp] of quads) { + // face normal from screen-space cross product (cheap, good enough) + const ux = PX[b] - PX[a], uy = PY[b] - PY[a]; + const vx = PX[d] - PX[a], vy = PY[d] - PY[a]; + const nzs = ux * vy - uy * vx; // z of cross -> facing + // world-ish normal proxy using depth gradient for lighting + let shade = 0.45 + 0.55 * Math.max(0, (nzs > 0 ? Ly : Ly * 0.6)); + shade = Math.min(1, shade); + // displacement tint: outward warm, inward cool + const dispN = Math.max(-1, Math.min(1, disp / (0.12 * span))); + const warm = dispN > 0 ? dispN : 0; + const cool = dispN < 0 ? -dispN : 0; + const r = metal[0] * shade + warm * 70 - cool * 30; + const g = metal[1] * shade - cool * 10; + const bl = metal[2] * shade + cool * 90; + ctx.fillStyle = `rgb(${Math.max(0, Math.min(255, r | 0))},${Math.max(0, Math.min(255, g | 0))},${Math.max(0, Math.min(255, bl | 0))})`; + ctx.beginPath(); + ctx.moveTo(PX[a], PY[a]); + ctx.lineTo(PX[b], PY[b]); + ctx.lineTo(PX[c], PY[c]); + ctx.lineTo(PX[d], PY[d]); + ctx.closePath(); + ctx.fill(); + // subtle wireframe + ctx.strokeStyle = "rgba(0,0,0,0.18)"; + ctx.lineWidth = 0.6; + ctx.stroke(); + } + + drawHUD(tm); + + await writeFrame(canvas.toBuffer("raw")); + if (f % 30 === 0) + process.stderr.write(`\rviz: frame ${f}/${totalFrames}`); + } + encoder.stdin.end(); + await new Promise((r) => encoder.on("close", r)); + process.stderr.write(`\nviz: wrote ${outPath}\n`); +} + +main().catch((e) => { + console.error(e); + process.exit(1); +}); diff --git a/pop/bell/c/.gitignore b/pop/bell/c/.gitignore new file mode 100644 index 000000000..f63f01c3d --- /dev/null +++ b/pop/bell/c/.gitignore @@ -0,0 +1,4 @@ +bell +*.o +*.tmp.wav +*.modes.json diff --git a/pop/bell/c/README.md b/pop/bell/c/README.md new file mode 100644 index 000000000..aea7413fa --- /dev/null +++ b/pop/bell/c/README.md @@ -0,0 +1,140 @@ +# bell — a physically-modeled bell + +A reusable bell voice whose sound comes from a real **finite-element model** of a +vibrating shell, so its **material parameters** (Young's modulus, density, +Poisson ratio, wall thickness, damping) and **geometry** genuinely shape the +timbre — not a table of hand-typed partial ratios. The same model drives a 3-D +visualization of the bell shell deforming by mode. + +Zero dependencies beyond libm (C) and `canvas`/ffmpeg (the JS viz). + +``` +pop/bell/ + c/bell.c, c/bell.h the engine + public API + c/build.sh cc -O3 -std=c11 -Wall -Wextra -o bell bell.c -lm + c/run-c.mjs pop entry point: render mp3 (+ optional viz) + c/compare.mjs C-vs-JS modal-render parity harness + bin/viz.mjs 3-D deforming-shell mp4 renderer + materials.json material presets (mirror of the C table) + geometries.json geometry presets (mirror of the C table) +``` + +## The physics + +A bell is a **surface of revolution**, so expanding the displacement field in a +Fourier series in the angular coordinate θ — `u(s,θ) = Σ_m u_m(s)·cos(mθ)` — +decouples the full 3-D shell into one **1-D meridian problem per circumferential +order `m`**. For each `m` we assemble small stiffness `K` and mass `M` matrices +along the meridian (conical-frustum thin-shell elements, 4 DOF/node: meridional +`u`, circumferential `v`, normal `w`, meridional rotation `β`; selective reduced +integration on the transverse shear avoids locking) and solve the generalized +symmetric eigenproblem + +``` +K φ = ω² M φ +``` + +with a hand-written Cholesky + cyclic-Jacobi eigensolver (no LAPACK). The +eigenvalues are the modal frequencies; the eigenvectors are the meridian mode +shapes. The bell **tone** is the rim-flexural family `m ≥ 2` — the hum, prime, +tierce, quint and nominal are all `m=2`/`m=3` modes with differing numbers of +nodal circles. (`m=0` breathing and `m=1` whole-body sway are excluded: a rim +strike barely excites them and they are not part of the tone.) + +- **Damping** comes from the material loss factor η: amplitude decay + `δ = π·f·η`, so `τ = 1/δ`. Constant η ⇒ higher partials decay faster — exactly + the bell-like behavior. +- **Strike** at the mouth rim sets each mode's initial amplitude from its + participation (its normal shape sampled at the strike point). +- **Pitch** is set by `bell_retune()`, a uniform geometric scale that multiplies + every modal frequency by the same factor so the strike note lands on the + requested pitch while the inharmonic ratio set is preserved. + +Material parameters behave correctly by construction: `K ∝ E`, `M ∝ ρ`, so pitch +scales as `√(E/ρ)`; thickness raises pitch; a larger bell lowers it as ~`1/size²`. + +## Validation (`./bell --selftest`) + +The eigensolver is gated against analytic limits before anything trusts it: + +1. Jacobi on a known symmetric matrix. +2. Generalized eig on a known `(K, M)` pair. +3. **Free-free Euler-Bernoulli beam** vs `(βL)²√(EI/ρA L⁴)` — validates + meridional bending (matches to ~1e-6). +4. **Cylinder → analytic in-plane ring flexural series** + `f_m = (1/2π)·(m(m²−1)/√(m²+1))·(h/a²)·√(E/12ρ)` for `m=2,3,4` (with ν=0) — + validates the hoop physics (**0.0 % error**). + +`compare.mjs` separately confirms the C render equals a JS reimplementation of +the same mode table to ~1e-8 (with the strike transient + normalization off). + +## CLI + +```bash +./build.sh +./bell --note A4 --material bronze --geometry church --dur 8 --out bell.wav +./bell ... --modes bell-modes.json # export geometry + modes + shapes (viz) +./bell ... --print-modes # print the partial table +./bell --selftest +``` + +Flags: `--note` (name like `C#5` or a bare Hz), `--material`, `--geometry`, +`--dur`, `--vel`, `--sr`, `--maxm`, `--nostrike`, `--nonorm` (last two for parity). + +### Pop pipeline + +```bash +node run-c.mjs --note A4 --material glass --geometry church \ + --out bell.mp3 --master bell --viz bell.mp4 +``` + +Materials: `bronze brass steel aluminum silver glass gold` +Geometries: `church handbell tubular bowl glass` + +### Visualization + +```bash +node bin/viz.mjs --modes bell-modes.json --audio bell.wav --out bell.mp4 \ + [--dur 9] [--fps 30] [--portrait] +``` + +A dependency-free software 3-D pipeline projects the deforming surface of +revolution (the displacement is the live sum of excited mode shapes as the sound +decays), with a HUD of material params, strike note and a per-partial spectrum +hued by circumferential order. Frames are BGRA → ffmpeg via +`pop/lib/preview-shared.mjs` (`spawnFFmpegEncode`), audio muxed in. + +## Public API (`bell.h`) + +```c +bell_geometry_preset(&g, "church"); +bell_material_preset(&mat, "bronze"); +bell_solve_modes(&g, &mat, /*max_m*/8, /*max_modes*/32, &modes); +bell_retune(&modes, 440.0); +bell_render(&modes, /*vel*/0.9, /*sr*/48000, /*dur*/8, L, R, nsamp); +bell_export_modes_json(&g, &mat, &modes, "modes.json"); +``` + +## Reuse across the monorepo + AC OS (follow-up) + +`bell.c`/`bell.h` are self-contained with the same shape as +`fedac/native/src/gm_synth.c`, which is compiled into the AC OS kernel +(`fedac/native/Makefile`) **and** symlinked into menuband +(`slab/menuband/Sources/CGMSynth/`). To make the bell live everywhere: + +1. Add `bell.c` to the `fedac/native` Makefile SRCS; call `bell_solve_modes` + once per voice config and `bell_render`-style modal playback in the audio mix + (modes can be precomputed at note-on; the eigensolve is sub-millisecond per + bell but is best cached). +2. Symlink `bell.c`/`bell.h` into a `slab/menuband/Sources/CBell/` target. +3. Optionally have `gm_synth.c`'s bell programs source their ratios from a + solved `BellModes` instead of the hand-tabulated `gm_chromperc_programs`. + +## Notes / honest limits + +The element is a faceted-conical thin-shell (Kirchhoff/Mindlin) reduced model. +It reproduces the analytic ring and beam limits exactly and gives a genuinely +inharmonic, material-driven bell spectrum, but it is not a research-grade shell +solver: the absolute ratio set of a *specific* historic bell is the product of +centuries of profile tuning. Edit the `build_profile` control points (or the +geometry presets) to chase a particular bell's hum/prime/tierce/nominal. diff --git a/pop/bell/c/bell.c b/pop/bell/c/bell.c new file mode 100644 index 000000000..2c74fc8db --- /dev/null +++ b/pop/bell/c/bell.c @@ -0,0 +1,901 @@ +// bell.c — physically-modeled bell voice. See bell.h for the overview. +// +// Build: ./build.sh (cc -O3 -std=c11 -Wall -Wextra -o bell bell.c -lm) +// +// Sections: +// 1. Small dense linear algebra (Cholesky, Jacobi, generalized eigensolve) +// 2. Analytic-limit self-tests (gate correctness before trusting the shell) +// 3. Axisymmetric shell FEM assembly -> appended in stage 2 +// 4. Strike / modal render / export -> appended in stage 2 +// 5. Presets + CLI -> appended in stage 2 + +#include "bell.h" +#include +#include +#include +#include + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif +#define TAU (2.0 * M_PI) + +// =========================================================================== +// 1. Linear algebra (row-major flat arrays, A[i*n + j]) +// =========================================================================== + +// Cholesky factor of SPD matrix M -> lower-triangular L (M = L L^T). +// Returns 0 on success, -1 if not positive-definite. +static int chol(const double *M, double *L, int n) { + for (int i = 0; i < n * n; i++) L[i] = 0.0; + for (int j = 0; j < n; j++) { + double d = M[j * n + j]; + for (int k = 0; k < j; k++) d -= L[j * n + k] * L[j * n + k]; + if (d <= 0.0) return -1; + L[j * n + j] = sqrt(d); + for (int i = j + 1; i < n; i++) { + double s = M[i * n + j]; + for (int k = 0; k < j; k++) s -= L[i * n + k] * L[j * n + k]; + L[i * n + j] = s / L[j * n + j]; + } + } + return 0; +} + +// Forward solve L y = b (L lower-triangular), in place into y. +static void fwd(const double *L, const double *b, double *y, int n) { + for (int i = 0; i < n; i++) { + double s = b[i]; + for (int k = 0; k < i; k++) s -= L[i * n + k] * y[k]; + y[i] = s / L[i * n + i]; + } +} + +// Back solve L^T x = b (L lower-triangular so L^T is upper), in place into x. +static void bwd_t(const double *L, const double *b, double *x, int n) { + for (int i = n - 1; i >= 0; i--) { + double s = b[i]; + for (int k = i + 1; k < n; k++) s -= L[k * n + i] * x[k]; + x[i] = s / L[i * n + i]; + } +} + +// Cyclic Jacobi eigensolve of symmetric A (n x n). Fills eval[n] and evec +// (columns are eigenvectors, evec[i*n + j] = component i of eigenvector j). +static void jacobi(double *A, int n, double *eval, double *evec) { + for (int i = 0; i < n; i++) + for (int j = 0; j < n; j++) evec[i * n + j] = (i == j) ? 1.0 : 0.0; + + for (int sweep = 0; sweep < 100; sweep++) { + double off = 0.0; + for (int p = 0; p < n; p++) + for (int q = p + 1; q < n; q++) off += A[p * n + q] * A[p * n + q]; + if (off < 1e-30) break; + + for (int p = 0; p < n; p++) { + for (int q = p + 1; q < n; q++) { + double apq = A[p * n + q]; + if (fabs(apq) < 1e-300) continue; + double app = A[p * n + p], aqq = A[q * n + q]; + double phi = 0.5 * atan2(2.0 * apq, aqq - app); + double c = cos(phi), s = sin(phi); + for (int k = 0; k < n; k++) { + double akp = A[k * n + p], akq = A[k * n + q]; + A[k * n + p] = c * akp - s * akq; + A[k * n + q] = s * akp + c * akq; + } + for (int k = 0; k < n; k++) { + double apk = A[p * n + k], aqk = A[q * n + k]; + A[p * n + k] = c * apk - s * aqk; + A[q * n + k] = s * apk + c * aqk; + } + for (int k = 0; k < n; k++) { + double vkp = evec[k * n + p], vkq = evec[k * n + q]; + evec[k * n + p] = c * vkp - s * vkq; + evec[k * n + q] = s * vkp + c * vkq; + } + } + } + } + for (int i = 0; i < n; i++) eval[i] = A[i * n + i]; +} + +// Generalized symmetric eigensolve K phi = lambda M phi (M SPD). +// eval[n] sorted ascending; evec columns are the mode shapes. Returns 0 ok. +static int gen_eig(const double *K, const double *M, int n, double *eval, + double *evec) { + double *L = malloc(sizeof(double) * n * n); + double *Y = malloc(sizeof(double) * n * n); + double *C = malloc(sizeof(double) * n * n); + double *W = malloc(sizeof(double) * n * n); + double *colb = malloc(sizeof(double) * n); + double *colx = malloc(sizeof(double) * n); + int rc = -1; + if (!L || !Y || !C || !W || !colb || !colx) goto done; + if (chol(M, L, n) != 0) goto done; + + // Y = L^-1 K (solve L Y = K column by column) + for (int j = 0; j < n; j++) { + for (int i = 0; i < n; i++) colb[i] = K[i * n + j]; + fwd(L, colb, colx, n); + for (int i = 0; i < n; i++) Y[i * n + j] = colx[i]; + } + // C = L^-1 Y^T = L^-1 K^T L^-T (= L^-1 K L^-T since K symmetric) + for (int j = 0; j < n; j++) { + for (int i = 0; i < n; i++) colb[i] = Y[j * n + i]; // row j of Y -> Y^T col j + fwd(L, colb, colx, n); + for (int i = 0; i < n; i++) C[i * n + j] = colx[i]; + } + // Symmetrize against round-off, then Jacobi. + for (int i = 0; i < n; i++) + for (int j = i + 1; j < n; j++) { + double a = 0.5 * (C[i * n + j] + C[j * n + i]); + C[i * n + j] = C[j * n + i] = a; + } + jacobi(C, n, eval, W); + + // phi = L^-T W (solve L^T phi = W column by column) + for (int j = 0; j < n; j++) { + for (int i = 0; i < n; i++) colb[i] = W[i * n + j]; + bwd_t(L, colb, colx, n); + for (int i = 0; i < n; i++) evec[i * n + j] = colx[i]; + } + + // Sort ascending by eigenvalue (selection sort; n is small). + for (int a = 0; a < n; a++) { + int mn = a; + for (int b = a + 1; b < n; b++) + if (eval[b] < eval[mn]) mn = b; + if (mn != a) { + double t = eval[a]; + eval[a] = eval[mn]; + eval[mn] = t; + for (int i = 0; i < n; i++) { + double tv = evec[i * n + a]; + evec[i * n + a] = evec[i * n + mn]; + evec[i * n + mn] = tv; + } + } + } + rc = 0; +done: + free(L); + free(Y); + free(C); + free(W); + free(colb); + free(colx); + return rc; +} + +// =========================================================================== +// 2. Analytic-limit self-tests +// =========================================================================== + +// Free-free Euler-Bernoulli beam, ne elements, length Ltot, bending EI, line +// mass rhoA. Fills the lowest `want` *nonzero* modal frequencies (Hz) into f[]. +static int beam_freefree(int ne, double Ltot, double EI, double rhoA, double *f, + int want) { + int nn = ne + 1, nd = 2 * nn; + double le = Ltot / ne; + double *K = calloc((size_t)nd * nd, sizeof(double)); + double *M = calloc((size_t)nd * nd, sizeof(double)); + double *ev = malloc(sizeof(double) * nd); + double *V = malloc(sizeof(double) * (size_t)nd * nd); + int rc = -1; + if (!K || !M || !ev || !V) goto done; + + double k0 = EI / (le * le * le); + double Ke[4][4] = {{12, 6 * le, -12, 6 * le}, + {6 * le, 4 * le * le, -6 * le, 2 * le * le}, + {-12, -6 * le, 12, -6 * le}, + {6 * le, 2 * le * le, -6 * le, 4 * le * le}}; + double m0 = rhoA * le / 420.0; + double Me[4][4] = {{156, 22 * le, 54, -13 * le}, + {22 * le, 4 * le * le, 13 * le, -3 * le * le}, + {54, 13 * le, 156, -22 * le}, + {-13 * le, -3 * le * le, -22 * le, 4 * le * le}}; + for (int e = 0; e < ne; e++) { + int map[4] = {2 * e, 2 * e + 1, 2 * e + 2, 2 * e + 3}; + for (int a = 0; a < 4; a++) + for (int b = 0; b < 4; b++) { + K[map[a] * nd + map[b]] += k0 * Ke[a][b]; + M[map[a] * nd + map[b]] += m0 * Me[a][b]; + } + } + if (gen_eig(K, M, nd, ev, V) != 0) goto done; + // Skip the two ~zero rigid-body modes. + int got = 0; + for (int i = 0; i < nd && got < want; i++) { + if (ev[i] < 1e-6) continue; + f[got++] = sqrt(ev[i]) / TAU; + } + rc = (got == want) ? 0 : -1; +done: + free(K); + free(M); + free(ev); + free(V); + return rc; +} + +int bell_selftest(int verbose) { + int fails = 0; + + // (a) Jacobi on [[2,1],[1,2]] -> eigenvalues {1,3}. + { + double A[4] = {2, 1, 1, 2}, ev[2], V[4]; + jacobi(A, 2, ev, V); + double lo = ev[0] < ev[1] ? ev[0] : ev[1]; + double hi = ev[0] < ev[1] ? ev[1] : ev[0]; + int ok = fabs(lo - 1.0) < 1e-9 && fabs(hi - 3.0) < 1e-9; + if (verbose) + printf(" [%s] Jacobi 2x2: %.6f, %.6f (want 1, 3)\n", ok ? "ok" : "FAIL", + lo, hi); + fails += !ok; + } + + // (b) Generalized eig: K=[[6,2],[2,3]], M=diag(2,1) -> {3-sqrt2, 3+sqrt2}. + { + double K[4] = {6, 2, 2, 3}, M[4] = {2, 0, 0, 1}, ev[2], V[4]; + int rc = gen_eig(K, M, 2, ev, V); + double want0 = 3.0 - sqrt(2.0), want1 = 3.0 + sqrt(2.0); + int ok = rc == 0 && fabs(ev[0] - want0) < 1e-9 && fabs(ev[1] - want1) < 1e-9; + if (verbose) + printf(" [%s] gen-eig 2x2: %.6f, %.6f (want %.6f, %.6f)\n", + ok ? "ok" : "FAIL", ev[0], ev[1], want0, want1); + fails += !ok; + } + + // (c) Free-free Euler-Bernoulli beam vs analytic (beta*L)^2 sqrt(EI/rhoA/L^4). + // beta*L for free-free: 4.730041, 7.853205, 10.995608. + { + double L = 1.0, EI = 1.0, rhoA = 1.0; + double f[3]; + int rc = beam_freefree(40, L, EI, rhoA, f, 3); + double bl[3] = {4.730040745, 7.853204624, 10.995607838}; + int ok = rc == 0; + if (verbose) printf(" [..] free-free beam (40 elems):\n"); + for (int i = 0; i < 3 && rc == 0; i++) { + double want = bl[i] * bl[i] * sqrt(EI / (rhoA * L * L * L * L)) / TAU; + double err = fabs(f[i] - want) / want; + int oki = err < 2e-3; // FEM discretization tolerance + ok = ok && oki; + if (verbose) + printf(" [%s] mode %d: %.5f Hz (want %.5f, err %.2e)\n", + oki ? "ok" : "FAIL", i + 1, f[i], want, err); + } + fails += !ok; + } + + // (d) Cylinder limit -> analytic in-plane ring flexural series. + // f_m = (1/2pi)(m(m^2-1)/sqrt(m^2+1)) (h/a^2) sqrt(E/(12 rho)). + // Uses nu=0 so the shell hoop stiffness D = E h^3/12 matches the 1-D ring. + { + double a = 0.20, h = 0.004, E = 105e9, rho = 8800.0; + BellGeometry g; + g.n = 48; + snprintf(g.name, sizeof(g.name), "test-cyl"); + for (int i = 0; i < g.n; i++) { + g.r[i] = a; + g.z[i] = 0.30 * (double)i / (g.n - 1); + g.h[i] = h; + } + BellMaterial mt = {E, rho, 0.0, 1e-3, "test"}; + BellModes md; + int rc = bell_solve_modes(&g, &mt, 5, 64, &md); + double pref = (h / (a * a)) * sqrt(E / (12.0 * rho)) / TAU; + if (verbose) printf(" [..] cylinder -> ring series (nu=0):\n"); + int ok = rc > 0; + for (int m = 2; m <= 4; m++) { + double lo = 1e30; + for (int k = 0; k < md.count; k++) + if (md.mode[k].m == m && md.mode[k].freq < lo) lo = md.mode[k].freq; + double want = pref * (m * (m * m - 1.0) / sqrt(m * m + 1.0)); + double err = fabs(lo - want) / want; + int oki = lo < 1e29 && err < 0.10; // faceted-conical + Flugge tolerance + ok = ok && oki; + if (verbose) + printf(" [%s] m=%d: %.2f Hz (ring %.2f, err %.1f%%)\n", + oki ? "ok" : "FAIL", m, lo, want, err * 100.0); + } + fails += !ok; + } + + if (verbose) + printf("selftest: %s (%d failure%s)\n", fails ? "FAIL" : "PASS", fails, + fails == 1 ? "" : "s"); + return fails; +} + +// =========================================================================== +// 3. Materials & geometry presets +// =========================================================================== + +void bell_default_material(BellMaterial *m) { + bell_material_preset(m, "bronze"); +} + +int bell_material_preset(BellMaterial *mat, const char *name) { + // E (Pa), rho (kg/m^3), nu, loss factor eta. Loss tuned for musical decay. + struct { + const char *n; + double E, rho, nu, loss; + } tbl[] = { + {"bronze", 105e9, 8800.0, 0.33, 2.0e-4}, // bell bronze (78Cu/22Sn) + {"brass", 100e9, 8500.0, 0.34, 3.0e-4}, + {"steel", 200e9, 7850.0, 0.30, 1.2e-4}, + {"aluminum", 69e9, 2700.0, 0.33, 8.0e-4}, + {"silver", 83e9, 10490.0, 0.37, 5.0e-4}, + {"glass", 70e9, 2500.0, 0.22, 6.0e-5}, // long shimmering ring + {"gold", 79e9, 19300.0, 0.42, 9.0e-4}, + }; + int N = (int)(sizeof(tbl) / sizeof(tbl[0])); + for (int i = 0; i < N; i++) { + if (strcmp(tbl[i].n, name) == 0) { + mat->E = tbl[i].E; + mat->rho = tbl[i].rho; + mat->nu = tbl[i].nu; + mat->loss = tbl[i].loss; + snprintf(mat->name, sizeof(mat->name), "%s", tbl[i].n); + return 0; + } + } + return -1; +} + +// Parametric meridian: crown (top, t=0) -> mouth (bottom, t=1). +// shape: 0 flaring bell, 1 straight tube, 2 hemispherical bowl. +static void build_profile(BellGeometry *g, const char *name, double rMouth, + double height, double rCrownFrac, double hMouth, + double hCrown, double flare, int shape) { + int n = 64; + g->n = n; + snprintf(g->name, sizeof(g->name), "%s", name); + for (int i = 0; i < n; i++) { + double t = (double)i / (n - 1); // 0 crown .. 1 mouth + double r, z; + if (shape == 1) { // tube + r = rMouth; + z = height * (1.0 - t); + } else if (shape == 2) { // bowl (hemisphere-ish) + double ang = t * (M_PI * 0.5); + r = rMouth * sin(ang); + z = height * (1.0 - cos(ang)); + } else { // flaring bell + double base = rCrownFrac + (1.0 - rCrownFrac) * pow(t, flare); + r = rMouth * base; + z = height * (1.0 - t); + } + if (r < 1e-4) r = 1e-4; + g->r[i] = r; + g->z[i] = z; + // Thicker toward the mouth/soundbow. + g->h[i] = hCrown + (hMouth - hCrown) * pow(t, 1.3); + } +} + +void bell_default_geometry(BellGeometry *g) { + bell_geometry_preset(g, "church"); +} + +int bell_geometry_preset(BellGeometry *g, const char *name) { + if (strcmp(name, "church") == 0) + build_profile(g, "church", 0.50, 0.80, 0.34, 0.045, 0.014, 1.7, 0); + else if (strcmp(name, "handbell") == 0) + build_profile(g, "handbell", 0.085, 0.13, 0.30, 0.006, 0.002, 1.8, 0); + else if (strcmp(name, "tubular") == 0) + build_profile(g, "tubular", 0.0159, 1.40, 1.0, 0.0012, 0.0012, 1.0, 1); + else if (strcmp(name, "bowl") == 0) + build_profile(g, "bowl", 0.11, 0.07, 1.0, 0.006, 0.006, 1.0, 2); + else if (strcmp(name, "glass") == 0) + build_profile(g, "glass", 0.045, 0.12, 0.55, 0.0018, 0.0014, 1.4, 0); + else + return -1; + return 0; +} + +// =========================================================================== +// 4. Axisymmetric shell FEM assembly (conical frustum elements, harmonic m) +// =========================================================================== +// +// DOF per node: U (meridional), V (circumferential), W (normal), beta +// (meridional section rotation). Local element order: +// [U1 V1 W1 b1 U2 V2 W2 b2] +// +// Strain-displacement (per element, faceted-conical): see README. Validated in +// the cylinder limit against the analytic ring series and in the straight-beam +// limit against free-free Euler-Bernoulli. Selective reduced integration on the +// transverse shear (1-pt) avoids shear locking. + +#define DOF_PER_NODE 4 + +// Assemble global K, M (nd x nd, nd = 4*n) for circumferential order m. +static void shell_assemble(const BellGeometry *g, const BellMaterial *mat, int m, + double *K, double *M, int nd) { + for (int i = 0; i < nd * nd; i++) K[i] = M[i] = 0.0; + double E = mat->E, nu = mat->nu, rho = mat->rho; + double G = E / (2.0 * (1.0 + nu)); + + // 2-point Gauss on [0,1]. + double gx[2] = {0.5 - 0.5 / sqrt(3.0), 0.5 + 0.5 / sqrt(3.0)}; + double gw[2] = {0.5, 0.5}; + + for (int e = 0; e < g->n - 1; e++) { + int i0 = e, i1 = e + 1; + double dr = g->r[i1] - g->r[i0]; + double dz = g->z[i1] - g->z[i0]; + double Le = sqrt(dr * dr + dz * dz); + if (Le < 1e-9) continue; + double sphi = dr / Le; // dr/ds + double cphi = dz / Le; // dz/ds + int map[8] = {DOF_PER_NODE * i0 + 0, DOF_PER_NODE * i0 + 1, + DOF_PER_NODE * i0 + 2, DOF_PER_NODE * i0 + 3, + DOF_PER_NODE * i1 + 0, DOF_PER_NODE * i1 + 1, + DOF_PER_NODE * i1 + 2, DOF_PER_NODE * i1 + 3}; + double Ke[8][8] = {{0}}, Me[8][8] = {{0}}; + + // --- membrane + bending (2-pt) --- + for (int gpt = 0; gpt < 2; gpt++) { + double xi = gx[gpt], w = gw[gpt]; + double N1 = 1.0 - xi, N2 = xi; + double dN = 1.0 / Le; // |dN1/ds| = dN, dN2/ds = +dN + double r = N1 * g->r[i0] + N2 * g->r[i1]; + double h = N1 * g->h[i0] + N2 * g->h[i1]; + if (r < 1e-5) r = 1e-5; + double Dm0 = E * h / (1.0 - nu * nu); + double Db0 = E * h * h * h / (12.0 * (1.0 - nu * nu)); + + // Membrane B-rows (3 x 8): es, eth, gsth. + double Bm[3][8] = {{0}}; + Bm[0][0] = -dN; Bm[0][4] = dN; // es=U' + Bm[1][0] = N1 * sphi / r; Bm[1][4] = N2 * sphi / r; // eth + Bm[1][1] = N1 * m / r; Bm[1][5] = N2 * m / r; + Bm[1][2] = N1 * cphi / r; Bm[1][6] = N2 * cphi / r; + Bm[2][1] = -dN - sphi / r * N1; Bm[2][5] = dN - sphi / r * N2; // gsth + Bm[2][0] = -(double)m / r * N1; Bm[2][4] = -(double)m / r * N2; + + // Bending B-rows (3 x 8): ks, kth, ksth. + double Bb[3][8] = {{0}}; + Bb[0][3] = -dN; Bb[0][7] = dN; // ks=beta' + double mm1 = (double)m * m - 1.0; // Flugge-consistent + Bb[1][2] = N1 * mm1 / (r * r); Bb[1][6] = N2 * mm1 / (r * r); // kth + Bb[1][3] = N1 * sphi / r; Bb[1][7] = N2 * sphi / r; + Bb[2][3] = N1 * (double)m / r; Bb[2][7] = N2 * (double)m / r; // ksth (twist) + + // Material matrices. + double Dm[3][3] = {{Dm0, Dm0 * nu, 0}, {Dm0 * nu, Dm0, 0}, + {0, 0, Dm0 * (1.0 - nu) / 2.0}}; + double Db[3][3] = {{Db0, Db0 * nu, 0}, {Db0 * nu, Db0, 0}, + {0, 0, Db0 * (1.0 - nu)}}; + double scale = w * Le * r; + for (int a = 0; a < 8; a++) + for (int b = 0; b < 8; b++) { + double s = 0.0; + for (int p = 0; p < 3; p++) + for (int q = 0; q < 3; q++) + s += Bm[p][a] * Dm[p][q] * Bm[q][b] + + Bb[p][a] * Db[p][q] * Bb[q][b]; + Ke[a][b] += scale * s; + } + + // Consistent mass (translational + small rotary inertia for beta). + double NU[8] = {N1, 0, 0, 0, N2, 0, 0, 0}; + double NV[8] = {0, N1, 0, 0, 0, N2, 0, 0}; + double NW[8] = {0, 0, N1, 0, 0, 0, N2, 0}; + double NB[8] = {0, 0, 0, N1, 0, 0, 0, N2}; + double rhoh = rho * h, rhoI = rho * h * h * h / 12.0; + for (int a = 0; a < 8; a++) + for (int b = 0; b < 8; b++) + Me[a][b] += scale * (rhoh * (NU[a] * NU[b] + NV[a] * NV[b] + + NW[a] * NW[b]) + + rhoI * NB[a] * NB[b]); + } + + // --- transverse shear (1-pt, reduced) --- + { + double N1 = 0.5, N2 = 0.5, dN = 1.0 / Le; + double r = N1 * g->r[i0] + N2 * g->r[i1]; + double h = N1 * g->h[i0] + N2 * g->h[i1]; + if (r < 1e-5) r = 1e-5; + double Ds0 = (5.0 / 6.0) * G * h; + double Bs[8] = {0}; + Bs[2] = -dN; Bs[6] = dN; // gsz = W' - beta + Bs[3] = -N1; Bs[7] = -N2; + double scale = 1.0 * Le * r; + for (int a = 0; a < 8; a++) + for (int b = 0; b < 8; b++) + Ke[a][b] += scale * Ds0 * Bs[a] * Bs[b]; + } + + for (int a = 0; a < 8; a++) + for (int b = 0; b < 8; b++) { + K[map[a] * nd + map[b]] += Ke[a][b]; + M[map[a] * nd + map[b]] += Me[a][b]; + } + } + + // Tiny diagonal mass floor so M stays SPD (guards starved DOFs at apex). + double mref = 0.0; + for (int i = 0; i < nd; i++) + if (M[i * nd + i] > mref) mref = M[i * nd + i]; + for (int i = 0; i < nd; i++) + if (M[i * nd + i] < 1e-9 * mref) M[i * nd + i] += 1e-9 * mref; +} + +// =========================================================================== +// 5. Solve modes +// =========================================================================== + +int bell_solve_modes(const BellGeometry *g, const BellMaterial *mat, int max_m, + int max_modes, BellModes *out) { + int n = g->n; + int nd = DOF_PER_NODE * n; + double *K = malloc(sizeof(double) * (size_t)nd * nd); + double *M = malloc(sizeof(double) * (size_t)nd * nd); + double *ev = malloc(sizeof(double) * nd); + double *V = malloc(sizeof(double) * (size_t)nd * nd); + if (!K || !M || !ev || !V) { + free(K); free(M); free(ev); free(V); + return -1; + } + out->count = 0; + + // The bell *tone* lives in the rim-flexural families (m>=2): hum, prime, + // tierce, quint, nominal are all m=2/m=3 modes with differing nodal circles. + // m=0 (breathing) and m=1 (whole-body sway/bend) are weakly struck and not + // part of the strike tone, so they are excluded. + for (int m = 2; m <= max_m; m++) { + shell_assemble(g, mat, m, K, M, nd); + if (gen_eig(K, M, nd, ev, V) != 0) continue; + // Reference (stiffest) frequency to threshold rigid-body modes. + double topw = ev[nd - 1] > 0 ? ev[nd - 1] : 1.0; + int keptThisM = 0; + for (int k = 0; k < nd && keptThisM < 6; k++) { + if (ev[k] <= 0) continue; + // Rigid-body modes are ~0 relative to the stiff membrane spectrum; soft + // bending modes sit ~1e-10*topw, so a deep threshold separates them. + if (ev[k] < 1e-12 * topw) continue; + double f = sqrt(ev[k]) / TAU; + if (f < 1.0 || f > 40000.0) continue; + if (out->count >= BELL_MAX_MODES) break; + BellMode *md = &out->mode[out->count]; + md->m = m; + md->freq = f; + md->nshape = n; + double wmax = 0.0; + for (int i = 0; i < n; i++) { + double Ui = V[(DOF_PER_NODE * i + 0) * nd + k]; + double Wi = V[(DOF_PER_NODE * i + 2) * nd + k]; + md->ushape[i] = Ui; + md->wshape[i] = Wi; + if (fabs(Wi) > wmax) wmax = fabs(Wi); + } + if (wmax < 1e-30) wmax = 1.0; + for (int i = 0; i < n; i++) { // normalize shapes to unit peak normal + md->wshape[i] /= wmax; + md->ushape[i] /= wmax; + } + md->part = fabs(md->wshape[0]); // strike at the mouth rim (node 0) + md->tau = 0.0; // filled after material loss below + out->count++; + keptThisM++; + } + } + + // Sort all retained modes by frequency. + for (int a = 0; a < out->count; a++) { + int mn = a; + for (int b = a + 1; b < out->count; b++) + if (out->mode[b].freq < out->mode[mn].freq) mn = b; + if (mn != a) { + BellMode t = out->mode[a]; + out->mode[a] = out->mode[mn]; + out->mode[mn] = t; + } + } + // Keep only the lowest max_modes. + if (max_modes > 0 && out->count > max_modes) out->count = max_modes; + + // Material damping: amplitude decay delta = pi f eta, tau = 1/delta. + for (int k = 0; k < out->count; k++) { + double delta = M_PI * out->mode[k].freq * mat->loss; + out->mode[k].tau = delta > 1e-9 ? 1.0 / delta : 1e9; + } + + // Strike note = loudest partial at the strike point (perceived pitch). + int best = 0; + double bestScore = -1.0; + for (int k = 0; k < out->count; k++) { + // Favor strong, low partials (perceptually carry the pitch). + double score = out->mode[k].part / (1.0 + out->mode[k].freq / 800.0); + if (score > bestScore) { + bestScore = score; + best = k; + } + } + out->strike_index = out->count ? best : 0; + out->strike_freq = out->count ? out->mode[best].freq : 0.0; + + free(K); free(M); free(ev); free(V); + return out->count; +} + +void bell_retune(BellModes *modes, double target_freq) { + if (modes->count == 0 || modes->strike_freq <= 0.0) return; + double s = target_freq / modes->strike_freq; + for (int k = 0; k < modes->count; k++) { + modes->mode[k].freq *= s; + double delta = modes->mode[k].tau > 0 ? 1.0 / modes->mode[k].tau : 0.0; + (void)delta; + // Recompute tau against the new frequency would change decay; keep decay + // tied to the *new* pitch so retuned bells of different sizes still ring + // believably: tau scales as 1/f at fixed loss factor. + modes->mode[k].tau /= s; + } + modes->strike_freq = target_freq; +} + +// =========================================================================== +// 6. Strike excitation + modal render +// =========================================================================== + +// Render toggles (set by the CLI; default on). Disabling both makes the output +// a pure deterministic modal sum for the compare.mjs JS-parity harness. +static int g_strike = 1; // add the clapper-contact noise transient +static int g_norm = 1; // peak-normalize the final mix + +long bell_render(const BellModes *modes, double strike_vel, double sr, + double dur, float *L, float *R, long nsamp) { + long N = (long)(dur * sr); + if (N > nsamp) N = nsamp; + for (long i = 0; i < N; i++) L[i] = R[i] = 0.0f; + if (modes->count == 0) return N; + + // Per-mode initial amplitude from strike participation. A_k = part * vel / w. + double amp[BELL_MAX_MODES], pan[BELL_MAX_MODES]; + double peak = 0.0; + for (int k = 0; k < modes->count; k++) { + double w = TAU * modes->mode[k].freq; + amp[k] = modes->mode[k].part * strike_vel / (w > 1e-9 ? sqrt(w) : 1.0); + peak += fabs(amp[k]); + // Stereo spread by circumferential order (m=2 center, higher = wider). + double sp = 0.16 * ((modes->mode[k].m % 5) - 2); + pan[k] = sp; + } + if (peak < 1e-30) peak = 1.0; + + for (int k = 0; k < modes->count; k++) { + double f = modes->mode[k].freq; + double tau = modes->mode[k].tau; + double a = amp[k] / peak; + double gl = 0.5 - 0.5 * pan[k]; // simple equal-ish pan + double gr = 0.5 + 0.5 * pan[k]; + double phase = 0.0, inc = f / sr; + for (long i = 0; i < N; i++) { + double env = exp(-(double)i / (tau * sr)); + double s = a * env * sin(TAU * phase); + L[i] += (float)(s * gl); + R[i] += (float)(s * gr); + phase += inc; + if (phase >= 1.0) phase -= 1.0; + } + } + + // Strike transient: brief filtered-noise click (the clapper contact). + if (g_strike) { + uint32_t rng = 0x1234567u; + double lp = 0.0; + long clk = (long)(0.006 * sr); // ~6 ms + double camp = 0.18 * strike_vel; + for (long i = 0; i < clk && i < N; i++) { + rng = rng * 1664525u + 1013904223u; + double white = ((double)rng / 4294967296.0) * 2.0 - 1.0; + lp = 0.35 * white + 0.65 * lp; + double env = exp(-(double)i / (0.0018 * sr)); + double s = camp * env * lp; + L[i] += (float)s; + R[i] += (float)s; + } + } + + // Peak normalize to 0.9. + double pk = 0.0; + for (long i = 0; i < N; i++) { + if (fabs(L[i]) > pk) pk = fabs(L[i]); + if (fabs(R[i]) > pk) pk = fabs(R[i]); + } + if (g_norm && pk > 1e-9) { + double gain = 0.9 / pk; + for (long i = 0; i < N; i++) { + L[i] = (float)(L[i] * gain); + R[i] = (float)(R[i] * gain); + } + } + return N; +} + +// =========================================================================== +// 7. JSON export (geometry + modes + shapes) for the viz +// =========================================================================== + +int bell_export_modes_json(const BellGeometry *g, const BellMaterial *mat, + const BellModes *modes, const char *path) { + FILE *f = fopen(path, "w"); + if (!f) return -1; + fprintf(f, "{\n"); + fprintf(f, " \"material\": {\"name\":\"%s\",\"E\":%g,\"rho\":%g,\"nu\":%g," + "\"loss\":%g},\n", + mat->name, mat->E, mat->rho, mat->nu, mat->loss); + fprintf(f, " \"geometry\": {\"name\":\"%s\",\"n\":%d,\n", g->name, g->n); + fprintf(f, " \"z\": ["); + for (int i = 0; i < g->n; i++) fprintf(f, "%s%.6f", i ? "," : "", g->z[i]); + fprintf(f, "],\n \"r\": ["); + for (int i = 0; i < g->n; i++) fprintf(f, "%s%.6f", i ? "," : "", g->r[i]); + fprintf(f, "],\n \"h\": ["); + for (int i = 0; i < g->n; i++) fprintf(f, "%s%.6f", i ? "," : "", g->h[i]); + fprintf(f, "]\n },\n"); + fprintf(f, " \"strike_index\": %d,\n", modes->strike_index); + fprintf(f, " \"strike_freq\": %.4f,\n", modes->strike_freq); + fprintf(f, " \"modes\": [\n"); + for (int k = 0; k < modes->count; k++) { + const BellMode *md = &modes->mode[k]; + fprintf(f, " {\"m\":%d,\"freq\":%.8f,\"tau\":%.8f,\"part\":%.8f,\n", + md->m, md->freq, md->tau, md->part); + fprintf(f, " \"w\": ["); + for (int i = 0; i < md->nshape; i++) + fprintf(f, "%s%.5f", i ? "," : "", md->wshape[i]); + fprintf(f, "],\n \"u\": ["); + for (int i = 0; i < md->nshape; i++) + fprintf(f, "%s%.5f", i ? "," : "", md->ushape[i]); + fprintf(f, "]}%s\n", k + 1 < modes->count ? "," : ""); + } + fprintf(f, " ]\n}\n"); + fclose(f); + return 0; +} + +// =========================================================================== +// 8. WAV writer + CLI +// =========================================================================== + +static void write_wav_f32_stereo(const char *path, const float *L, + const float *R, long n, int sr) { + FILE *f = fopen(path, "wb"); + if (!f) return; + int ch = 2, bits = 32; + int byteRate = sr * ch * bits / 8; + int blockAlign = ch * bits / 8; + long dataBytes = n * ch * (bits / 8); + fwrite("RIFF", 1, 4, f); + uint32_t riff = 36 + (uint32_t)dataBytes; + fwrite(&riff, 4, 1, f); + fwrite("WAVE", 1, 4, f); + fwrite("fmt ", 1, 4, f); + uint32_t fmtlen = 16; + uint16_t fmt = 3; // IEEE float + uint16_t chs = (uint16_t)ch; + uint32_t srr = (uint32_t)sr; + uint16_t ba = (uint16_t)blockAlign, bps = (uint16_t)bits; + uint32_t br = (uint32_t)byteRate; + fwrite(&fmtlen, 4, 1, f); + fwrite(&fmt, 2, 1, f); + fwrite(&chs, 2, 1, f); + fwrite(&srr, 4, 1, f); + fwrite(&br, 4, 1, f); + fwrite(&ba, 2, 1, f); + fwrite(&bps, 2, 1, f); + fwrite("data", 1, 4, f); + uint32_t dlen = (uint32_t)dataBytes; + fwrite(&dlen, 4, 1, f); + for (long i = 0; i < n; i++) { + fwrite(&L[i], 4, 1, f); + fwrite(&R[i], 4, 1, f); + } + fclose(f); +} + +static double note_to_freq(const char *s) { + // Accept "A4", "C#5", "Db3", or a bare number (Hz). + char *end; + double num = strtod(s, &end); + if (end != s && *end == '\0') return num; + int sem; + switch (s[0]) { + case 'C': sem = 0; break; + case 'D': sem = 2; break; + case 'E': sem = 4; break; + case 'F': sem = 5; break; + case 'G': sem = 7; break; + case 'A': sem = 9; break; + case 'B': sem = 11; break; + default: return 440.0; + } + int i = 1; + if (s[i] == '#') { sem++; i++; } + else if (s[i] == 'b') { sem--; i++; } + int oct = atoi(&s[i]); + int midi = (oct + 1) * 12 + sem; + return 440.0 * pow(2.0, (midi - 69) / 12.0); +} + +int main(int argc, char **argv) { + const char *note = "A4", *material = "bronze", *geom = "church"; + const char *out = NULL, *modesPath = NULL; + double dur = 6.0, vel = 0.9; + double sr = 48000.0; + int maxM = 8, maxModes = 32; + int showModes = 0; + + for (int i = 1; i < argc; i++) { + if (strcmp(argv[i], "--selftest") == 0) return bell_selftest(1); + else if (strcmp(argv[i], "--note") == 0 && i + 1 < argc) note = argv[++i]; + else if (strcmp(argv[i], "--material") == 0 && i + 1 < argc) material = argv[++i]; + else if (strcmp(argv[i], "--geometry") == 0 && i + 1 < argc) geom = argv[++i]; + else if (strcmp(argv[i], "--dur") == 0 && i + 1 < argc) dur = atof(argv[++i]); + else if (strcmp(argv[i], "--vel") == 0 && i + 1 < argc) vel = atof(argv[++i]); + else if (strcmp(argv[i], "--sr") == 0 && i + 1 < argc) sr = atof(argv[++i]); + else if (strcmp(argv[i], "--maxm") == 0 && i + 1 < argc) maxM = atoi(argv[++i]); + else if (strcmp(argv[i], "--out") == 0 && i + 1 < argc) out = argv[++i]; + else if (strcmp(argv[i], "--modes") == 0 && i + 1 < argc) modesPath = argv[++i]; + else if (strcmp(argv[i], "--print-modes") == 0) showModes = 1; + else if (strcmp(argv[i], "--nostrike") == 0) g_strike = 0; + else if (strcmp(argv[i], "--nonorm") == 0) g_norm = 0; + else { + fprintf(stderr, "bell: unknown arg %s\n", argv[i]); + return 2; + } + } + + BellGeometry g; + BellMaterial mat; + if (bell_geometry_preset(&g, geom) != 0) { + fprintf(stderr, "bell: unknown geometry '%s'\n", geom); + return 2; + } + if (bell_material_preset(&mat, material) != 0) { + fprintf(stderr, "bell: unknown material '%s'\n", material); + return 2; + } + + BellModes modes; + int nm = bell_solve_modes(&g, &mat, maxM, maxModes, &modes); + if (nm <= 0) { + fprintf(stderr, "bell: mode solve failed\n"); + return 1; + } + double target = note_to_freq(note); + bell_retune(&modes, target); + + fprintf(stderr, + "bell: %s/%s, %d modes, strike note %.2f Hz (%s), pre-tune %.2f Hz\n", + geom, material, nm, modes.strike_freq, note, + modes.strike_freq / (target / (modes.mode[modes.strike_index].freq))); + + if (showModes) { + fprintf(stderr, " idx m freq(Hz) ratio tau(s) part\n"); + double f0 = modes.mode[modes.strike_index].freq; + for (int k = 0; k < modes.count; k++) + fprintf(stderr, " %3d %2d %9.2f %6.3f %7.2f %.3f\n", k, + modes.mode[k].m, modes.mode[k].freq, modes.mode[k].freq / f0, + modes.mode[k].tau, modes.mode[k].part); + } + + if (modesPath) { + if (bell_export_modes_json(&g, &mat, &modes, modesPath) == 0) + fprintf(stderr, "bell: wrote modes -> %s\n", modesPath); + } + + if (out) { + long nsamp = (long)(dur * sr) + 16; + float *Lb = malloc(sizeof(float) * nsamp); + float *Rb = malloc(sizeof(float) * nsamp); + long n = bell_render(&modes, vel, sr, dur, Lb, Rb, nsamp); + write_wav_f32_stereo(out, Lb, Rb, n, (int)sr); + fprintf(stderr, "bell: wrote %ld samples -> %s\n", n, out); + free(Lb); + free(Rb); + } + + return 0; +} diff --git a/pop/bell/c/bell.h b/pop/bell/c/bell.h new file mode 100644 index 000000000..4739d7e82 --- /dev/null +++ b/pop/bell/c/bell.h @@ -0,0 +1,92 @@ +// bell.h — physically-modeled bell voice (axisymmetric thin-shell FEM). +// +// A bell is a surface of revolution, so the 3-D shell problem decouples, by a +// Fourier expansion in the angular coordinate theta, into one 1-D meridian +// problem per circumferential order m. For each m we assemble small stiffness +// (K) and mass (M) matrices along the meridian and solve the generalized +// symmetric eigenproblem K phi = omega^2 M phi. The eigenvalues are the modal +// frequencies; the eigenvectors are the meridian mode shapes — the same shapes +// the 3-D visualization animates. Material parameters (Young's modulus, density, +// Poisson ratio, wall thickness, damping) enter every term, so the timbre and +// decay genuinely follow the physics. +// +// Zero dependencies beyond libm. See README.md for the formulation + citations. + +#ifndef BELL_H +#define BELL_H + +#define BELL_MAX_NODES 128 // meridian discretization stations +#define BELL_MAX_MODES 64 // partials retained after the eigensolve + +// --- Material --------------------------------------------------------------- +// E in pascals, rho in kg/m^3, nu dimensionless, loss = damping loss factor eta +// (amplitude decay rate delta = pi * f * eta, so tau = 1 / delta). +typedef struct { + double E; + double rho; + double nu; + double loss; + char name[32]; +} BellMaterial; + +// --- Geometry --------------------------------------------------------------- +// A meridian profile: n stations sampled bottom (mouth) to top (crown), each +// with axial coord z, radius r, and wall thickness h — all in metres. +typedef struct { + int n; + double z[BELL_MAX_NODES]; + double r[BELL_MAX_NODES]; + double h[BELL_MAX_NODES]; + char name[32]; +} BellGeometry; + +// --- A single vibrational mode --------------------------------------------- +typedef struct { + int m; // circumferential order (number of nodal meridians) + double freq; // Hz + double tau; // amplitude e-fold time (s), from material loss factor + double part; // strike participation (normal shape sampled at the strike point) + int nshape; // number of meridian shape samples (== geometry n) + double wshape[BELL_MAX_NODES]; // normal (radial) displacement along meridian + double ushape[BELL_MAX_NODES]; // meridional displacement along meridian +} BellMode; + +// --- The full solved mode set ---------------------------------------------- +typedef struct { + int count; + BellMode mode[BELL_MAX_MODES]; + int strike_index; // index of the nominal / "strike note" mode + double strike_freq; // its frequency (Hz) +} BellModes; + +// --- Presets ---------------------------------------------------------------- +void bell_default_geometry(BellGeometry *g); // a church-bell profile +int bell_geometry_preset(BellGeometry *g, const char *name); +void bell_default_material(BellMaterial *mat); // bell bronze +int bell_material_preset(BellMaterial *mat, const char *name); + +// --- Solve / tune / render -------------------------------------------------- +// Returns mode count (>0) on success, <0 on failure. max_m caps circumferential +// orders examined; max_modes caps retained partials. +int bell_solve_modes(const BellGeometry *g, const BellMaterial *mat, int max_m, + int max_modes, BellModes *out); + +// Uniformly rescale every modal frequency so the strike note equals target_freq. +// Physical (a uniform geometric scale shifts all modes by the same factor), so +// the inharmonic ratio set is preserved. +void bell_retune(BellModes *modes, double target_freq); + +// Render a single strike into stereo float buffers (length nsamp each). Returns +// the number of samples written. strike_vel in [0,1]. +long bell_render(const BellModes *modes, double strike_vel, double sr, + double dur, float *L, float *R, long nsamp); + +// Export geometry + solved modes (incl. mode shapes) as JSON for the viz. +int bell_export_modes_json(const BellGeometry *g, const BellMaterial *mat, + const BellModes *modes, const char *path); + +// --- Self-test: validates the eigensolver against analytic limits ---------- +// Returns 0 if all checks pass within tolerance, nonzero otherwise. +int bell_selftest(int verbose); + +#endif // BELL_H diff --git a/pop/bell/c/build.sh b/pop/bell/c/build.sh new file mode 100755 index 000000000..c1484be0f --- /dev/null +++ b/pop/bell/c/build.sh @@ -0,0 +1,6 @@ +#!/usr/bin/env bash +# Build the bell engine. No -ffast-math (keeps summation stable / JS-parity). +set -e +cd "$(dirname "$0")" +cc -O3 -std=c11 -Wall -Wextra -o bell bell.c -lm +echo "built ./bell" diff --git a/pop/bell/c/compare.mjs b/pop/bell/c/compare.mjs new file mode 100755 index 000000000..2367b857b --- /dev/null +++ b/pop/bell/c/compare.mjs @@ -0,0 +1,107 @@ +#!/usr/bin/env node +// compare.mjs — parity harness: the C modal render vs a JS reference rebuilt +// from the SAME exported mode table. Mirrors the convention in +// pop/nullabye/c/compare.mjs and pop/marimba/c/compare.mjs. +// +// Both sides render with the strike transient and normalization disabled +// (--nostrike --nonorm) so the comparison is a clean deterministic modal sum; +// agreement then validates that the JS-side modal synth (used elsewhere in the +// monorepo) matches the C engine within libm tolerance. +// +// Usage: node pop/bell/c/compare.mjs [--note A4] [--material bronze] + +import { spawnSync } from "node:child_process"; +import { existsSync, mkdirSync, readFileSync, statSync } from "node:fs"; +import { dirname, resolve } from "node:path"; +import { fileURLToPath } from "node:url"; + +const HERE = dirname(fileURLToPath(import.meta.url)); +const ENGINE = resolve(HERE, "bell"); +const args = process.argv.slice(2); +const val = (k, d) => { + const i = args.indexOf(k); + return i >= 0 ? args[i + 1] : d; +}; +const note = val("--note", "A4"); +const material = val("--material", "bronze"); +const geometry = val("--geometry", "church"); +const dur = parseFloat(val("--dur", "3")); +const SR = 48000; + +if (!existsSync(ENGINE) || statSync(resolve(HERE, "bell.c")).mtimeMs > statSync(ENGINE).mtimeMs) { + spawnSync("bash", [resolve(HERE, "build.sh")], { stdio: "inherit" }); +} + +mkdirSync("/tmp/bell-cmp", { recursive: true }); +const wav = "/tmp/bell-cmp/c.wav"; +const modesPath = "/tmp/bell-cmp/modes.json"; + +// C render: pure modal sum (no strike, no normalize). +const r = spawnSync(ENGINE, [ + "--note", note, "--material", material, "--geometry", geometry, + "--dur", String(dur), "--vel", "0.9", + "--nostrike", "--nonorm", "--out", wav, "--modes", modesPath, +], { stdio: ["ignore", "ignore", "inherit"] }); +if (r.status !== 0) { console.error("C render failed"); process.exit(1); } + +// Read the C WAV (f32 stereo) -> left channel. +function readWavF32Left(path) { + const buf = readFileSync(path); + // find "data" chunk + let off = 12; + while (off + 8 <= buf.length) { + const id = buf.toString("ascii", off, off + 4); + const sz = buf.readUInt32LE(off + 4); + if (id === "data") { + const n = sz / 8; // stereo f32 + const L = new Float64Array(n); + for (let i = 0; i < n; i++) L[i] = buf.readFloatLE(off + 8 + i * 8); + return L; + } + off += 8 + sz + (sz & 1); + } + throw new Error("no data chunk"); +} +const cL = readWavF32Left(wav); + +// JS reference: rebuild from the modes JSON exactly as bell_render does +// (left channel). amp_k = part * vel / sqrt(2*pi*f); peak = sum|amp|; then +// s = (amp/peak) * exp(-i/(tau*sr)) * sin(2*pi*phase) * panLeft. +const model = JSON.parse(readFileSync(modesPath, "utf8")); +const modes = model.modes; +const vel = 0.9; +const TAU = 2 * Math.PI; +const amp = modes.map((m) => m.part * vel / Math.sqrt(TAU * m.freq)); +let peak = 0; +for (const a of amp) peak += Math.abs(a); +if (peak < 1e-30) peak = 1; +const N = Math.min(cL.length, Math.floor(dur * SR)); +const jL = new Float64Array(N); +for (let k = 0; k < modes.length; k++) { + const m = modes[k]; + const a = amp[k] / peak; + const sp = 0.16 * ((m.m % 5) - 2); + const gl = 0.5 - 0.5 * sp; // panLeft (matches C) + const inc = m.freq / SR; + let phase = 0; + for (let i = 0; i < N; i++) { + const env = Math.exp(-i / (m.tau * SR)); + jL[i] += a * env * Math.sin(TAU * phase) * gl; + phase += inc; + if (phase >= 1) phase -= 1; + } +} + +let linf = 0, l2 = 0; +for (let i = 0; i < N; i++) { + const d = Math.abs(cL[i] - jL[i]); + if (d > linf) linf = d; + l2 += d * d; +} +l2 = Math.sqrt(l2 / N); +const TOL = 5e-4; // JSON-precision + libm ulps +const ok = linf < TOL; +console.log(`compare ${geometry}/${material} ${note}: ${modes.length} modes, ${N} samples`); +console.log(` L-inf = ${linf.toExponential(3)} RMS = ${l2.toExponential(3)} tol = ${TOL}`); +console.log(ok ? " ✓ C and JS modal renders agree" : " ✗ MISMATCH"); +process.exit(ok ? 0 : 1); diff --git a/pop/bell/c/run-c.mjs b/pop/bell/c/run-c.mjs new file mode 100755 index 000000000..189502d03 --- /dev/null +++ b/pop/bell/c/run-c.mjs @@ -0,0 +1,100 @@ +#!/usr/bin/env node +// run-c.mjs — pop-facing entry point for the physically-modeled bell. +// +// Builds the C engine if stale, renders a strike (geometry + material -> WAV), +// masters it to mp3, and can also emit the modes JSON + the 3-D visualization. +// +// Usage: +// node pop/bell/c/run-c.mjs --note A4 --material bronze --geometry church \ +// --out /tmp/bell.mp3 [--master bell|bright|raw] [--dur 8] [--vel 0.9] +// ... --viz /tmp/bell.mp4 also render the 3-D mode visualization +// ... --print-modes print the solved partial table +// +// Materials: bronze brass steel aluminum silver glass gold +// Geometries: church handbell tubular bowl glass + +import { spawnSync } from "node:child_process"; +import { existsSync, mkdirSync, statSync, unlinkSync } from "node:fs"; +import { dirname, resolve } from "node:path"; +import { fileURLToPath } from "node:url"; + +const HERE = dirname(fileURLToPath(import.meta.url)); +const ENGINE = resolve(HERE, "bell"); +const VIZ = resolve(HERE, "..", "bin", "viz.mjs"); + +const args = process.argv.slice(2); +const val = (k, d = null) => { + const i = args.indexOf(k); + return i >= 0 && i + 1 < args.length ? args[i + 1] : d; +}; +const note = val("--note", "A4"); +const material = val("--material", "bronze"); +const geometry = val("--geometry", "church"); +const dur = val("--dur", "8"); +const vel = val("--vel", "0.9"); +const outMp3 = val("--out", "/tmp/bell.mp3"); +const master = val("--master", "bell"); +const vizOut = val("--viz", null); +const printModes = args.includes("--print-modes"); + +// Build engine if missing or stale. +const cSrc = resolve(HERE, "bell.c"); +if (!existsSync(ENGINE) || statSync(cSrc).mtimeMs > statSync(ENGINE).mtimeMs) { + console.log("[run-c] building bell…"); + const b = spawnSync("bash", [resolve(HERE, "build.sh")], { stdio: "inherit" }); + if (b.status !== 0) process.exit(1); +} + +mkdirSync(dirname(resolve(outMp3)), { recursive: true }); +const wav = `${outMp3}.tmp.wav`; +const modesJson = vizOut ? `${outMp3}.modes.json` : null; + +const engArgs = [ + "--note", note, "--material", material, "--geometry", geometry, + "--dur", dur, "--vel", vel, "--out", wav, +]; +if (modesJson) engArgs.push("--modes", modesJson); +if (printModes) engArgs.push("--print-modes"); +const r = spawnSync(ENGINE, engArgs, { stdio: "inherit" }); +if (r.status !== 0) { console.error("✗ bell engine failed"); process.exit(1); } + +// Master chains (the WAV is already peak-normalized; these add polish). +const MASTERS = { + // Keep the long shimmering tail, gentle glue, a touch of air. + bell: [ + "highpass=f=40", + "acompressor=threshold=-20dB:ratio=2:attack=20:release=300:makeup=1.5:knee=8", + "treble=g=1.2:f=9000", + "alimiter=limit=0.95:attack=4:release=120", + ], + // More presence for clangier/percussive uses. + bright: [ + "highpass=f=60", + "equalizer=f=3500:t=q:w=1.2:g=2", + "acompressor=threshold=-18dB:ratio=2.6:attack=8:release=180:makeup=2.5:knee=6", + "treble=g=2:f=8500", + "alimiter=limit=0.97:attack=3:release=80", + ], + raw: [], +}; +const chain = MASTERS[master]; +if (!chain) { console.error(`unknown master: ${master}`); process.exit(1); } + +const ffArgs = ["-hide_banner", "-y", "-loglevel", "error", "-i", wav]; +if (chain.length) ffArgs.push("-af", chain.join(",")); +ffArgs.push("-c:a", "libmp3lame", "-q:a", "2", resolve(outMp3)); +const ff = spawnSync("ffmpeg", ffArgs, { stdio: "inherit" }); +if (ff.status !== 0) { console.error("✗ ffmpeg failed"); process.exit(1); } +console.log(`✓ ${outMp3} (bell · ${geometry}/${material} · ${master}-mastered)`); + +// Optional 3-D visualization (uses the rendered WAV for audio + reactive flash). +if (vizOut) { + const v = spawnSync("node", [ + VIZ, "--modes", modesJson, "--audio", wav, "--out", resolve(vizOut), + "--dur", dur, + ], { stdio: "inherit" }); + if (v.status !== 0) { console.error("✗ viz failed"); process.exit(1); } + console.log(`✓ ${vizOut} (3-D mode visualization)`); +} + +try { unlinkSync(wav); } catch {} diff --git a/pop/bell/geometries.json b/pop/bell/geometries.json new file mode 100644 index 000000000..dd6c178c4 --- /dev/null +++ b/pop/bell/geometries.json @@ -0,0 +1,8 @@ +{ + "_note": "Mirrors bell_geometry_preset()/build_profile() in c/bell.c — keep in sync. All lengths in metres. shape: 0 flaring bell, 1 straight tube, 2 hemispherical bowl. The meridian is sampled crown(top)->mouth(bottom); thickness tapers thicker toward the mouth/soundbow.", + "church": { "rMouth": 0.50, "height": 0.80, "rCrownFrac": 0.34, "hMouth": 0.045, "hCrown": 0.014, "flare": 1.7, "shape": 0, "about": "large flaring church/tower bell" }, + "handbell": { "rMouth": 0.085, "height": 0.13, "rCrownFrac": 0.30, "hMouth": 0.006, "hCrown": 0.002, "flare": 1.8, "shape": 0, "about": "small handbell, bright + quick" }, + "tubular": { "rMouth": 0.0159, "height": 1.40, "rCrownFrac": 1.0, "hMouth": 0.0012, "hCrown": 0.0012, "flare": 1.0, "shape": 1, "about": "orchestral tubular bell (a tuned pipe)" }, + "bowl": { "rMouth": 0.11, "height": 0.07, "rCrownFrac": 1.0, "hMouth": 0.006, "hCrown": 0.006, "flare": 1.0, "shape": 2, "about": "singing bowl — long, pure, beating partials" }, + "glass": { "rMouth": 0.045, "height": 0.12, "rCrownFrac": 0.55, "hMouth": 0.0018, "hCrown": 0.0014, "flare": 1.4, "shape": 0, "about": "thin glass — pair with the glass material" } +} diff --git a/pop/bell/materials.json b/pop/bell/materials.json new file mode 100644 index 000000000..fb6a32382 --- /dev/null +++ b/pop/bell/materials.json @@ -0,0 +1,10 @@ +{ + "_note": "Mirrors bell_material_preset() in c/bell.c — keep in sync. E in Pa, rho in kg/m^3, nu dimensionless, loss = damping loss factor eta (sets decay; higher = faster). Pitch is set by E/rho/geometry; decay by loss.", + "bronze": { "E": 1.05e11, "rho": 8800, "nu": 0.33, "loss": 2.0e-4, "about": "bell bronze (78Cu/22Sn) — the classic" }, + "brass": { "E": 1.00e11, "rho": 8500, "nu": 0.34, "loss": 3.0e-4, "about": "softer, shorter ring" }, + "steel": { "E": 2.00e11, "rho": 7850, "nu": 0.30, "loss": 1.2e-4, "about": "bright, high pitch for size, long ring" }, + "aluminum": { "E": 0.69e11, "rho": 2700, "nu": 0.33, "loss": 8.0e-4, "about": "light, dull, fast decay" }, + "silver": { "E": 0.83e11, "rho": 10490, "nu": 0.37, "loss": 5.0e-4, "about": "dense, mellow" }, + "glass": { "E": 0.70e11, "rho": 2500, "nu": 0.22, "loss": 6.0e-5, "about": "very long shimmering ring; low nu shifts inharmonicity" }, + "gold": { "E": 0.79e11, "rho": 19300, "nu": 0.42, "loss": 9.0e-4, "about": "extremely dense, low + short" } +} diff --git a/pop/lib/bell.mjs b/pop/lib/bell.mjs new file mode 100644 index 000000000..faac41e22 --- /dev/null +++ b/pop/lib/bell.mjs @@ -0,0 +1,149 @@ +// bell.mjs — the physically-modeled bell as a /pop voice. +// +// Wraps the validated C engine at pop/bell/c/bell (FEM shell modal synthesis) +// so any track renderer can drop bell strikes into a mix. Returns plain +// Float32Array buffers at 48 kHz — the /pop convention. The engine is built +// on first use if missing or stale. +// +// import { renderBell, renderBellMono, bellSampleBank } from "../lib/bell.mjs"; +// const { L, R } = renderBell({ note: "A4", material: "bronze", dur: 6 }); +// // ...mix L/R into your stereo bus at the strike time. +// +// Materials: bronze brass steel aluminum silver glass gold (pitch + decay) +// Geometries: church handbell tubular bowl glass (timbre) +// +// Pitch is set by `note` (name like "C#5", or a bare Hz number via `freq`); +// `material`/`geometry` shape the timbre and ring; `vel` (0..1) the strike. + +import { spawnSync } from "node:child_process"; +import { existsSync, statSync, mkdtempSync, readFileSync, rmSync } from "node:fs"; +import { tmpdir } from "node:os"; +import { dirname, resolve, join } from "node:path"; +import { fileURLToPath } from "node:url"; + +const HERE = dirname(fileURLToPath(import.meta.url)); +const ENGINE = resolve(HERE, "..", "bell", "c", "bell"); +const CSRC = resolve(HERE, "..", "bell", "c", "bell.c"); +const BUILD = resolve(HERE, "..", "bell", "c", "build.sh"); + +export const BELL_MATERIALS = [ + "bronze", "brass", "steel", "aluminum", "silver", "glass", "gold", +]; +export const BELL_GEOMETRIES = ["church", "handbell", "tubular", "bowl", "glass"]; +export const BELL_SR = 48_000; + +// Build the engine if it is missing or older than its source. +function ensureEngine() { + const stale = + !existsSync(ENGINE) || + (existsSync(CSRC) && statSync(CSRC).mtimeMs > statSync(ENGINE).mtimeMs); + if (stale) { + const b = spawnSync("bash", [BUILD], { stdio: "inherit" }); + if (b.status !== 0) throw new Error("bell: engine build failed"); + } +} + +// Read a stereo f32 WAV into { L, R, sampleRate }. +function readWavStereoF32(path) { + const buf = readFileSync(path); + let ch = 2, sr = BELL_SR, bits = 32, fmt = 3, i = 12; + while (i < buf.length - 8) { + const id = buf.toString("ascii", i, i + 4); + const size = buf.readUInt32LE(i + 4); + const body = i + 8; + if (id === "fmt ") { + fmt = buf.readUInt16LE(body); + ch = buf.readUInt16LE(body + 2); + sr = buf.readUInt32LE(body + 4); + bits = buf.readUInt16LE(body + 14); + } else if (id === "data") { + const bps = bits / 8; + const frames = Math.floor(size / (bps * ch)); + const L = new Float32Array(frames); + const R = new Float32Array(frames); + for (let f = 0; f < frames; f++) { + const o = body + f * ch * bps; + const rd = (k) => + fmt === 3 ? buf.readFloatLE(o + k * bps) : buf.readInt16LE(o + k * bps) / 32768; + L[f] = rd(0); + R[f] = ch > 1 ? rd(1) : L[f]; + } + return { L, R, sampleRate: sr }; + } + i = body + size + (size & 1); + } + throw new Error(`bell: no data chunk in ${path}`); +} + +function engineArgs({ note, freq, material, geometry, dur, vel, sr, maxm }) { + const pitch = freq != null ? String(freq) : note ?? "A4"; + const a = ["--note", pitch]; + if (material) a.push("--material", material); + if (geometry) a.push("--geometry", geometry); + if (dur != null) a.push("--dur", String(dur)); + if (vel != null) a.push("--vel", String(vel)); + if (sr != null) a.push("--sr", String(sr)); + if (maxm != null) a.push("--maxm", String(maxm)); + return a; +} + +// Render one strike -> { L, R, sampleRate } stereo Float32Arrays. +export function renderBell(opts = {}) { + ensureEngine(); + const tmp = mkdtempSync(join(tmpdir(), "bell-")); + const wav = join(tmp, "out.wav"); + try { + const r = spawnSync(ENGINE, [...engineArgs(opts), "--out", wav], { + stdio: ["ignore", "ignore", "inherit"], + }); + if (r.status !== 0) throw new Error("bell: render failed"); + return readWavStereoF32(wav); + } finally { + rmSync(tmp, { recursive: true, force: true }); + } +} + +// Convenience mono downmix (many track buses are mono per voice). +export function renderBellMono(opts = {}) { + const { L, R, sampleRate } = renderBell(opts); + const m = new Float32Array(L.length); + for (let i = 0; i < L.length; i++) m[i] = 0.5 * (L[i] + R[i]); + return { samples: m, sampleRate }; +} + +// Solve + return the modes (and geometry/material) without rendering audio — +// useful for inspection or feeding pop/bell/bin/viz.mjs. +export function bellModes(opts = {}) { + ensureEngine(); + const tmp = mkdtempSync(join(tmpdir(), "bell-")); + const json = join(tmp, "modes.json"); + try { + const r = spawnSync(ENGINE, [...engineArgs(opts), "--modes", json], { + stdio: ["ignore", "ignore", "inherit"], + }); + if (r.status !== 0) throw new Error("bell: solve failed"); + return JSON.parse(readFileSync(json, "utf8")); + } finally { + rmSync(tmp, { recursive: true, force: true }); + } +} + +// Render a bank of pitched one-shots to WAV files (for sample-playback +// pipelines). `notes` is an array of note names or Hz. Returns the paths. +export function bellSampleBank({ notes, material = "bronze", geometry = "church", dur = 6, dir }) { + ensureEngine(); + if (!dir) throw new Error("bell: bellSampleBank needs a `dir`"); + const paths = []; + for (const n of notes) { + const safe = String(n).replace(/[^a-zA-Z0-9.+-]/g, "_"); + const out = join(dir, `bell-${geometry}-${material}-${safe}.wav`); + const r = spawnSync( + ENGINE, + [...engineArgs({ note: n, material, geometry, dur }), "--out", out], + { stdio: ["ignore", "ignore", "inherit"] }, + ); + if (r.status !== 0) throw new Error(`bell: bank render failed for ${n}`); + paths.push(out); + } + return paths; +}