"""Independent numpy implementation of the player's modal engine, for checking the worklet. Reads a modes file, builds the contact patches exactly as the app does (linear interpolation on the mesh, best-fit face plane), and integrates the modal resonators with the same scheme: complex state z per mode with a = Im z, exact pole per (sub-)step, zero-order-hold forcing, symplectic Euler for the exciter, implicit (Newton) contact forces. Output is the unbaffled listener pressure at fs. python modal_reference.py modes.bin hits.json out.f32 [--fs 44100] [--seconds 0.5] """ from __future__ import annotations import argparse import json import struct import numpy as np def read_bin(path): data = open(path, "rb").read() assert data[:4] == b"CYMB" (hlen,) = struct.unpack("= -1e-9) & (l2 >= -1e-9) & (l3 >= -1e-9) e = int(np.argmax(inside)) assert inside[e] return tri[e], np.array([l1[e], l2[e], l3[e]]) STICK = dict(mass=0.03, k=5e8, alpha=1.5, patch_length=0.03, patch_points=7, patch_angle=0.6, gap=1e-4) MALLET = dict(mass=0.3, alpha=2.5, compression=12e-3, radius=0.012, gap=1e-4) def build_patch(arrays, kind, r, phi_angle, v, mallet_points=7): xy, tri, z0 = arrays["xy"].astype(np.float64), arrays["tri"], arrays["z0"] phi = (arrays["phi_q"].astype(np.float32) * arrays["phi_scale"][None, :].astype(np.float32)) xs, ys = r * np.cos(phi_angle), r * np.sin(phi_angle) if kind == "stick": ex = STICK ang = phi_angle + ex["patch_angle"] xi = np.linspace(-0.5, 0.5, ex["patch_points"]) * ex["patch_length"] px, py = xs + xi * np.cos(ang), ys + xi * np.sin(ang) k_total, alpha = ex["k"], ex["alpha"] else: ex = MALLET alpha = ex["alpha"] k_total = (alpha + 1) * 0.5 * ex["mass"] * v ** 2 / ex["compression"] ** (alpha + 1) th = np.arange(6 if mallet_points > 1 else 0) * np.pi / 3 px = np.concatenate([[xs], xs + ex["radius"] * np.cos(th)]) py = np.concatenate([[ys], ys + ex["radius"] * np.sin(th)]) P = len(px) phip = np.zeros((P, phi.shape[1]), np.float32) z0p = np.zeros(P) for p in range(P): nodes, w = locate(xy, tri, px[p], py[p]) acc = np.zeros(phi.shape[1], np.float32) for i in range(3): acc += np.float32(w[i]) * arrays["phi_q"][nodes[i]].astype(np.float32) phip[p] = acc * arrays["phi_scale"].astype(np.float32) z0p[p] = float(w @ z0[nodes]) A = np.column_stack([px - xs, py - ys, np.ones(P)]) offsets = z0p - A @ np.linalg.lstsq(A, z0p, rcond=None)[0] return dict(phi=phip.astype(np.float64), offsets=offsets.astype(np.float32).astype(np.float64), k_point=k_total / P, alpha=alpha, mass=ex["mass"], gap=ex["gap"], velocity=v) def main(): ap = argparse.ArgumentParser() ap.add_argument("modes"); ap.add_argument("hits"); ap.add_argument("out") ap.add_argument("--fs", type=int, default=44100); ap.add_argument("--seconds", type=float, default=0.5) ap.add_argument("--signal", choices=["unbaffled", "front"], default="unbaffled") ap.add_argument("--no-direct", action="store_true", help="drop the direct term Re(G) f") ap.add_argument("--mallet-points", type=int, default=7) args = ap.parse_args() meta, arrays = read_bin(args.modes) fs, B = args.fs, 128 freq = arrays["freq"] mat = meta["config"]["material"] sigma = mat["sigma0"] + mat["sigma1"] * arrays["damp_q"] freq = freq * np.sqrt((mat["E"] / meta["E0"]) / (mat["rho"] / meta["rho0"])) gb = arrays["gain_back"] if args.signal == "unbaffled" else np.zeros_like(arrays["gain_back"]) G = (arrays["gain_front"][:, 0] - gb[:, 0]) + 1j * (arrays["gain_front"][:, 1] - gb[:, 1]) G = G.astype(np.complex128) w0 = 2 * np.pi * freq wd = np.sqrt(np.maximum(w0 ** 2 - sigma ** 2, 1e-6)) lam = -sigma + 1j * wd C = G * lam ** 2 D = G.real * (0.0 if args.no_direct else 1.0) dt = 1.0 / fs pole1 = np.exp(lam * dt) hits = json.load(open(args.hits)) b1 = (pole1 - 1) / (lam * wd) n_blocks = int(np.ceil(args.seconds * fs / B)) z = np.zeros(len(freq), np.complex128) out = np.zeros(n_blocks * B) active = [] n = 0 for b in range(n_blocks): for h in hits: if round(h["t"] * fs / B) == b: p = build_patch(arrays, h["kind"], h["r"], h["phi"], h["v"], args.mallet_points) # the worklet computes the compliance from float32 phi in float64 arithmetic p["c"] = (p["phi"] * b1.imag[None, :]) @ p["phi"].T p["F"] = np.zeros(len(p["offsets"])) p["released"] = False; p["touched"] = False; p["peak"] = 0.0; p["samples"] = 0; p["kind"] = h["kind"] active.append(p) for ex in active: if not ex["released"]: wp = ex["phi"] @ z.imag + ex["offsets"] ex["zs"] = wp.max() + ex["gap"]; ex["vs"] = -ex["velocity"]; ex["released"] = True if not active: for k in range(B): z *= pole1 out[n + k] = (C * z).imag.sum() else: for k in range(B): z = pole1 * z f = np.zeros(len(freq)) any_force = False for ex in active: wp = ex["phi"] @ z.imag + ex["offsets"] base = ex["zs"] + dt * ex["vs"] em = dt * dt / ex["mass"] F = ex["F"] if (wp > base).any(): kp, a, c = ex["k_point"], ex["alpha"], ex["c"] P = len(F) for it in range(12): d = wp - base - em * F.sum() - c @ F pos = d > 0 g = np.where(pos, F - kp * np.maximum(d, 0) ** a, F) J = np.eye(P) s = np.where(pos, kp * a * np.maximum(d, 1e-300) ** (a - 1), 0.0) J += s[:, None] * (c + em) J[~pos, :] = np.eye(P)[~pos, :] dF = np.linalg.solve(J, g) F = np.maximum(0, F - dF) if np.abs(dF).max() < 1e-7 * max(1.0, F.max()): break else: F = np.zeros_like(F) ex["F"] = F Ft = F.sum() if Ft > 0: any_force = True; ex["touched"] = True; ex["samples"] += 1; ex["peak"] = max(ex["peak"], Ft) f -= F @ ex["phi"] ex["vs"] += dt * Ft / ex["mass"]; ex["zs"] += dt * ex["vs"] if any_force: z = z + b1 * f out[n + k] = (C * z).imag.sum() + (D * f).sum() keep = [] for ex in active: top = (ex["phi"] @ z.imag + ex["offsets"]).max() done = ex["touched"] and ex["vs"] > 0 and ex["zs"] - top > 2e-3 if not done: keep.append(ex) else: print(f" {ex['kind']}: peak force {ex['peak']:.0f} N, contact {ex['samples'] / fs * 1e3:.2f} ms, rebound {ex['vs']:+.2f} m/s") active = keep n += B out.astype(np.float32).tofile(args.out) print(f"wrote {args.out}: peak {np.abs(out).max():.3f} Pa") if __name__ == "__main__": main()