1"""Independent numpy implementation of the player's modal engine, for checking the worklet.
3Reads a modes file, builds the contact patches exactly as the app does (linear interpolation on
4the mesh, best-fit face plane), and integrates the modal resonators with the same scheme: complex
5state z per mode with a = Im z, exact pole per (sub-)step, zero-order-hold forcing, symplectic
6Euler for the exciter, implicit (Newton) contact forces. Output is the unbaffled listener pressure at fs.
8 python modal_reference.py modes.bin hits.json out.f32 [--fs 44100] [--seconds 0.5]
9"""
10from __future__ import annotations
12import argparse
13import json
14import struct
16import numpy as np
19def read_bin(path):
20 data = open(path, "rb").read()
21 assert data[:4] == b"CYMB"
22 (hlen,) = struct.unpack("<I", data[4:8])
23 meta = json.loads(data[8:8 + hlen].decode())
24 base = 8 + hlen
25 arrays = {}
26 for name, e in meta.pop("arrays").items():
27 arrays[name] = np.frombuffer(data, dtype=np.dtype(e["dtype"]).newbyteorder("<"),
28 count=int(np.prod(e["shape"])), offset=base + e["offset"]).reshape(e["shape"])
29 return meta, arrays
32def locate(xy, tri, px, py):
33 x, y = xy[tri, 0], xy[tri, 1]
34 det = (x[:, 1] - x[:, 0]) * (y[:, 2] - y[:, 0]) - (x[:, 2] - x[:, 0]) * (y[:, 1] - y[:, 0])
35 l1 = ((y[:, 1] - y[:, 2]) * (px - x[:, 2]) + (x[:, 2] - x[:, 1]) * (py - y[:, 2])) / det
36 l2 = ((y[:, 2] - y[:, 0]) * (px - x[:, 2]) + (x[:, 0] - x[:, 2]) * (py - y[:, 2])) / det
37 l3 = 1 - l1 - l2
38 inside = (l1 >= -1e-9) & (l2 >= -1e-9) & (l3 >= -1e-9)
39 e = int(np.argmax(inside))
40 assert inside[e]
41 return tri[e], np.array([l1[e], l2[e], l3[e]])
44STICK = dict(mass=0.03, k=5e8, alpha=1.5, patch_length=0.03, patch_points=7, patch_angle=0.6, gap=1e-4)
45MALLET = dict(mass=0.3, alpha=2.5, compression=12e-3, radius=0.012, gap=1e-4)
48def build_patch(arrays, kind, r, phi_angle, v, mallet_points=7):
49 xy, tri, z0 = arrays["xy"].astype(np.float64), arrays["tri"], arrays["z0"]
50 phi = (arrays["phi_q"].astype(np.float32) * arrays["phi_scale"][None, :].astype(np.float32))
51 xs, ys = r * np.cos(phi_angle), r * np.sin(phi_angle)
52 if kind == "stick":
53 ex = STICK
54 ang = phi_angle + ex["patch_angle"]
55 xi = np.linspace(-0.5, 0.5, ex["patch_points"]) * ex["patch_length"]
56 px, py = xs + xi * np.cos(ang), ys + xi * np.sin(ang)
57 k_total, alpha = ex["k"], ex["alpha"]
58 else:
59 ex = MALLET
60 alpha = ex["alpha"]
61 k_total = (alpha + 1) * 0.5 * ex["mass"] * v ** 2 / ex["compression"] ** (alpha + 1)
62 th = np.arange(6 if mallet_points > 1 else 0) * np.pi / 3
63 px = np.concatenate([[xs], xs + ex["radius"] * np.cos(th)])
64 py = np.concatenate([[ys], ys + ex["radius"] * np.sin(th)])
65 P = len(px)
66 phip = np.zeros((P, phi.shape[1]), np.float32)
67 z0p = np.zeros(P)
68 for p in range(P):
69 nodes, w = locate(xy, tri, px[p], py[p])
70 acc = np.zeros(phi.shape[1], np.float32)
71 for i in range(3):
72 acc += np.float32(w[i]) * arrays["phi_q"][nodes[i]].astype(np.float32)
73 phip[p] = acc * arrays["phi_scale"].astype(np.float32)
74 z0p[p] = float(w @ z0[nodes])
75 A = np.column_stack([px - xs, py - ys, np.ones(P)])
76 offsets = z0p - A @ np.linalg.lstsq(A, z0p, rcond=None)[0]
77 return dict(phi=phip.astype(np.float64), offsets=offsets.astype(np.float32).astype(np.float64),
78 k_point=k_total / P, alpha=alpha, mass=ex["mass"], gap=ex["gap"], velocity=v)
81def main():
82 ap = argparse.ArgumentParser()
83 ap.add_argument("modes"); ap.add_argument("hits"); ap.add_argument("out")
84 ap.add_argument("--fs", type=int, default=44100); ap.add_argument("--seconds", type=float, default=0.5)
85 ap.add_argument("--signal", choices=["unbaffled", "front"], default="unbaffled")
86 ap.add_argument("--no-direct", action="store_true", help="drop the direct term Re(G) f")
87 ap.add_argument("--mallet-points", type=int, default=7)
88 args = ap.parse_args()
89 meta, arrays = read_bin(args.modes)
90 fs, B = args.fs, 128
91 freq = arrays["freq"]
92 mat = meta["config"]["material"]
93 sigma = mat["sigma0"] + mat["sigma1"] * arrays["damp_q"]
94 freq = freq * np.sqrt((mat["E"] / meta["E0"]) / (mat["rho"] / meta["rho0"]))
95 gb = arrays["gain_back"] if args.signal == "unbaffled" else np.zeros_like(arrays["gain_back"])
96 G = (arrays["gain_front"][:, 0] - gb[:, 0]) + 1j * (arrays["gain_front"][:, 1] - gb[:, 1])
97 G = G.astype(np.complex128)
98 w0 = 2 * np.pi * freq
99 wd = np.sqrt(np.maximum(w0 ** 2 - sigma ** 2, 1e-6))
100 lam = -sigma + 1j * wd
101 C = G * lam ** 2
102 D = G.real * (0.0 if args.no_direct else 1.0)
103 dt = 1.0 / fs
104 pole1 = np.exp(lam * dt)
106 hits = json.load(open(args.hits))
107 b1 = (pole1 - 1) / (lam * wd)
108 n_blocks = int(np.ceil(args.seconds * fs / B))
109 z = np.zeros(len(freq), np.complex128)
110 out = np.zeros(n_blocks * B)
111 active = []
112 n = 0
113 for b in range(n_blocks):
114 for h in hits:
115 if round(h["t"] * fs / B) == b:
116 p = build_patch(arrays, h["kind"], h["r"], h["phi"], h["v"], args.mallet_points)
117 # the worklet computes the compliance from float32 phi in float64 arithmetic
118 p["c"] = (p["phi"] * b1.imag[None, :]) @ p["phi"].T
119 p["F"] = np.zeros(len(p["offsets"]))
120 p["released"] = False; p["touched"] = False; p["peak"] = 0.0; p["samples"] = 0; p["kind"] = h["kind"]
121 active.append(p)
122 for ex in active:
123 if not ex["released"]:
124 wp = ex["phi"] @ z.imag + ex["offsets"]
125 ex["zs"] = wp.max() + ex["gap"]; ex["vs"] = -ex["velocity"]; ex["released"] = True
126 if not active:
127 for k in range(B):
128 z *= pole1
129 out[n + k] = (C * z).imag.sum()
130 else:
131 for k in range(B):
132 z = pole1 * z
133 f = np.zeros(len(freq))
134 any_force = False
135 for ex in active:
136 wp = ex["phi"] @ z.imag + ex["offsets"]
137 base = ex["zs"] + dt * ex["vs"]
138 em = dt * dt / ex["mass"]
139 F = ex["F"]
140 if (wp > base).any():
141 kp, a, c = ex["k_point"], ex["alpha"], ex["c"]
142 P = len(F)
143 for it in range(12):
144 d = wp - base - em * F.sum() - c @ F
145 pos = d > 0
146 g = np.where(pos, F - kp * np.maximum(d, 0) ** a, F)
147 J = np.eye(P)
148 s = np.where(pos, kp * a * np.maximum(d, 1e-300) ** (a - 1), 0.0)
149 J += s[:, None] * (c + em)
150 J[~pos, :] = np.eye(P)[~pos, :]
151 dF = np.linalg.solve(J, g)
152 F = np.maximum(0, F - dF)
153 if np.abs(dF).max() < 1e-7 * max(1.0, F.max()):
154 break
155 else:
156 F = np.zeros_like(F)
157 ex["F"] = F
158 Ft = F.sum()
159 if Ft > 0:
160 any_force = True; ex["touched"] = True; ex["samples"] += 1; ex["peak"] = max(ex["peak"], Ft)
161 f -= F @ ex["phi"]
162 ex["vs"] += dt * Ft / ex["mass"]; ex["zs"] += dt * ex["vs"]
163 if any_force:
164 z = z + b1 * f
165 out[n + k] = (C * z).imag.sum() + (D * f).sum()
166 keep = []
167 for ex in active:
168 top = (ex["phi"] @ z.imag + ex["offsets"]).max()
169 done = ex["touched"] and ex["vs"] > 0 and ex["zs"] - top > 2e-3
170 if not done:
171 keep.append(ex)
172 else:
173 print(f" {ex['kind']}: peak force {ex['peak']:.0f} N, contact {ex['samples'] / fs * 1e3:.2f} ms, rebound {ex['vs']:+.2f} m/s")
174 active = keep
175 n += B
176 out.astype(np.float32).tofile(args.out)
177 print(f"wrote {args.out}: peak {np.abs(out).max():.3f} Pa")
180if __name__ == "__main__":
181 main()