4dff647Bead-sheet cymbal: 2D FPUT plate with bending and dome, chaos scan, web appJeremy Magland 1"""Where does the struck sheet become chaotic? Pairs of strikes differing by 1e-6 in strength, with the web app's
2damping and time stepping; how different are they after a second?
4 WORKERS=32 python chaos_scan.py # results/chaos_scan.npz and a printed table
6For each setting: bending Kb, dome Hd, strike energy E, ring time of the highest mode (lowest mode 6 s), 1.5 s of
7sound with the lowest mode at 100 Hz. Reported: relative difference of the pickup velocity in 1.0-1.5 s, and the
8RMS difference of the final energy share per log-spaced band (dB, bands holding > 0.1%).
9"""
11from __future__ import annotations
13import itertools
14import os
15from multiprocessing import Pool
17os.environ.setdefault("OMP_NUM_THREADS", "1")
19import numpy as np # noqa: E402
21from plate import Plate, strike # noqa: E402
23HERE = os.path.dirname(os.path.abspath(__file__))
24RES = os.path.join(HERE, "results")
25N, SR, F1, DUR, DELTA = 24, 44100, 100.0, 1.5, 1e-6
28def one(args):
29 Kb, Hd, E, thi = args
30 pl0 = Plate(N=N, Kb=Kb, Hd=Hd)
31 om = pl0.omega()
32 om1, omN = om.min(), om.max()
33 tau = om1 / (2 * np.pi * F1)
34 M = N + 1
35 s1, sN = 4 * np.sin(np.pi / (2 * M)) ** 2, 4 * np.sin(N * np.pi / (2 * M)) ** 2
36 lam1, lamN = 2 * s1, 2 * sN
37 sLo, sHi = np.log(1000) / 6.0 * tau, np.log(1000) / thi * tau
38 c1 = max(0.0, 2 * (sHi - sLo) / (lamN - lam1))
39 c0 = max(0.0, 2 * sLo - c1 * lam1)
40 pl = Plate(N=N, Kb=Kb, Hd=Hd, c0=c0, c1=c1)
41 dtOut = 1 / SR / tau
42 nsub = int(np.ceil(dtOut / (1.2 / omN)))
43 h = dtOut / nsub
44 w, v = strike(pl, E, 9, 10, 1.6, M=2, scale=[1.0, 1 + DELTA])
45 n = int(DUR * SR)
46 sig = np.zeros((n, 2))
47 f = pl.force(w)
48 for i in range(n):
49 for _ in range(nsub):
50 v = v + 0.5 * h * f
51 w = w + h * v
52 f = pl.force(w)
53 v = v + 0.5 * h * f
54 v = v - h * (pl.c0 * v - pl.c1 * pl.lap(v))
55 if not np.isfinite(v).all():
56 return args, np.nan, np.nan, True
57 sig[i] = v[:, 15, 17]
58 a, b = sig[int(1.0 * SR):, 0], sig[int(1.0 * SR):, 1]
59 wave = float(np.linalg.norm(a - b) / np.linalg.norm(a))
60 Ek = pl.mode_energies(w, v).reshape(2, -1)
61 fk = (om / om1 * F1).ravel()
62 edges = np.geomspace(F1, fk.max() * 1.0001, 25)
63 bands = np.stack([Ek[:, (fk >= lo) & (fk < hi)].sum(1) for lo, hi in zip(edges[:-1], edges[1:])], 1)
64 bands /= bands.sum(1, keepdims=True)
65 keep = bands.mean(0) > 1e-3
66 spec = float(np.sqrt(np.mean((10 * np.log10(bands[0, keep] / bands[1, keep])) ** 2)))
67 return args, wave, spec, False
70if __name__ == "__main__":
71 os.makedirs(RES, exist_ok=True)
72 grid = list(itertools.product((8.0, 25.0), (0.0, 10.0), (1e4, 3e4, 1e5), (0.8, 2.0, 5.0)))
73 rows = []
74 with Pool(int(os.environ.get("WORKERS", "32"))) as pool:
75 for args, wave, spec, blew in pool.imap_unordered(one, grid):
76 rows.append((*args, wave, spec, blew))
77 Kb, Hd, E, thi = args
78 print(f"Kb {Kb:4g} dome {Hd:4g} E {E:7g} high ring {thi:3g} s: "
79 + ("BLEW UP" if blew else f"waveform diff 1-1.5 s {wave:.1e}, final band spectra differ {spec:.2f} dB"), flush=True)
80 np.savez(os.path.join(RES, "chaos_scan.npz"), rows=np.array(rows, dtype=float))