1"""Exploration of the struck wire (wire.py). Writes results/<name>.npz; plots.py draws them.
3 python experiments.py explore # undamped: bending x dome x strike energy; high-band share over time
5Time is reported in seconds of sound with the lowest mode at F_LOW Hz, as in the web app.
6"""
8from __future__ import annotations
10import os
11import sys
12from multiprocessing import Pool
14os.environ.setdefault("OMP_NUM_THREADS", "1")
16import numpy as np # noqa: E402
18from wire import Wire, run, strike # noqa: E402
20HERE = os.path.dirname(os.path.abspath(__file__))
21RES = os.path.join(HERE, "results")
22N, CENTER, WIDTH, F_LOW = 64, 14, 3.0, 100.0
25def dt_for(w: Wire, E):
26 om, _ = w.modes()
27 return 0.4 / np.sqrt(om[-1] ** 2 + 12 * E)
30def tau_for(w: Wire):
31 om, _ = w.modes()
32 return om[0] / (2 * np.pi * F_LOW) # seconds of sound per time unit
35def _explore_one(args):
36 Kb, Hd, E = args
37 w = Wire(N=N, Kb=Kb, Hd=Hd)
38 tau = tau_for(w)
39 dt = dt_for(w, E)
40 t_end = 0.3 / tau
41 rec = max(1, int(round(0.002 / tau / dt)))
42 y0, v0 = strike(w, E, CENTER, WIDTH, M=2, scale=[1.0, 1 + 1e-10])
43 r = run(w, y0, v0, t_end, dt, record_every=rec)
44 f = r["omega"] / (2 * np.pi * tau) # mode frequencies in Hz of sound
45 return Kb, Hd, E, r["t"] * tau, r["Ek"][:, 0, :], f, r["E"][:, 0]
48def explore(pool):
49 jobs = [(Kb, Hd, E) for Kb in (0.0, 4.0) for Hd in (0.0, 3.0, 8.0) for E in (2.0, 10.0, 50.0)]
50 out = {}
51 for Kb, Hd, E, t, Ek, f, Etot in pool.imap_unordered(_explore_one, jobs):
52 k = f"Kb{Kb:g}_Hd{Hd:g}_E{E:g}"
53 out[k + "_t"], out[k + "_Ek"], out[k + "_f"], out[k + "_E"] = t, Ek, f, Etot
54 share = Ek / Ek.sum(1, keepdims=True)
55 hi1, hi4 = share[:, f > 1000].sum(1), share[:, f > 4000].sum(1)
56 print(f" Kb {Kb:g} Hd {Hd:g} E {E:g}: top {f[-1]:.0f} Hz; share above 1 kHz at 10/50/300 ms "
57 f"{hi1[4]:.2f} {hi1[24]:.2f} {hi1[-1]:.2f}; above 4 kHz {hi4[4]:.3f} {hi4[24]:.3f} {hi4[-1]:.3f}; "
58 f"energy drift {abs(Etot[-1] / Etot[0] - 1):.1e}", flush=True)
59 np.savez(os.path.join(RES, "explore.npz"), **out)
62if __name__ == "__main__":
63 os.makedirs(RES, exist_ok=True)
64 which = sys.argv[1:] or ["explore"]
65 with Pool(int(os.environ.get("WORKERS", "18"))) as pool:
66 for w in which:
67 print(w, flush=True)
68 dict(explore=explore)[w](pool)