1"""Plots for analysis/chaos-check.ts: python3 analysis/plot_chaos.py chaos-check.json out_dir"""
3import json
4import sys
5from pathlib import Path
7import matplotlib
9matplotlib.use("Agg")
10import matplotlib.pyplot as plt
11import numpy as np
13d = json.load(open(sys.argv[1]))
14out = Path(sys.argv[2] if len(sys.argv) > 2 else ".")
15out.mkdir(parents=True, exist_ok=True)
16t = np.array(d["time"]) * 1000 # ms
17TMAX = 30
20def sep_axes(ax, title):
21 ax.set_yscale("log")
22 ax.set_xlim(0, TMAX)
23 ax.set_ylim(1e-12, 3e-2)
24 ax.set_xlabel("time (ms)")
25 ax.set_ylabel("RMS separation of the beads (m)")
26 ax.set_title(title, fontsize=10)
27 ax.grid(alpha=0.3)
30# 1. step size
31fig, ax = plt.subplots(figsize=(6.5, 4))
32for s in d["stepSize"]:
33 y = np.maximum(np.array(s["d"]), 1e-13)
34 ax.plot(t, y, lw=1, label=f"cfl {s['cfl']} vs {s['against']}")
35sep_axes(ax, "Same pluck (10 mm), time step halved: runs agree for longer")
36ax.legend(fontsize=8)
37fig.tight_layout()
38fig.savefig(out / "step-size.png", dpi=150)
40# 2. perturbations
41fig, axs = plt.subplots(1, 2, figsize=(11, 4), sharey=True)
42for ax, cfl in zip(axs, sorted({p["cfl"] for p in d["perturbation"]}, reverse=True)):
43 for p in [p for p in d["perturbation"] if p["cfl"] == cfl]:
44 y = np.maximum(np.array(p["d"]), 1e-13)
45 g = p["growth"]
46 lab = f"eps {p['eps']:g}" + (f", growth {g['lambda']:.0f}/s" if g else "")
47 ax.plot(t, y, lw=1, label=lab)
48 sep_axes(ax, f"Pluck amplitude times (1 + eps), cfl {cfl}")
49 ax.legend(fontsize=8)
50fig.tight_layout()
51fig.savefig(out / "perturbation.png", dpi=150)
53# 3. energy error vs time step
54fig, ax = plt.subplots(figsize=(5, 4))
55dt = np.array([r["dt"] for r in d["runs"]])
56err = np.array([r["energyError"] for r in d["runs"]])
57ax.loglog(dt, err, "o-", label="largest energy balance error / peak energy")
58ax.loglog(dt, err[0] * (dt / dt[0]) ** 2, "k--", lw=0.8, label="slope 2")
59ax.set_xlabel("time step (s)")
60ax.set_ylabel("relative energy error")
61ax.legend(fontsize=8)
62ax.grid(alpha=0.3, which="both")
63fig.tight_layout()
64fig.savefig(out / "energy.png", dpi=150)
66# 4. sound statistics: octave bands, step sizes vs perturbations
67fig, ax = plt.subplots(figsize=(6.5, 4))
68for r in d["runs"]:
69 ax.plot([b["f"] for b in r["bands"]], [b["db"] for b in r["bands"]], "o-", color="C0", alpha=0.6, lw=1)
70for p in d["perturbation"]:
71 ax.plot([b["f"] for b in p["bands"]], [b["db"] for b in p["bands"]], "s--", color="C1", alpha=0.6, lw=1)
72ax.plot([], [], "o-", color="C0", label="6 time steps")
73ax.plot([], [], "s--", color="C1", label="8 tiny perturbations")
74ax.set_xscale("log")
75ax.set_xlabel("octave band centre (Hz)")
76ax.set_ylabel("level, 1 s to 2 s (dB re loudest band)")
77ax.set_title("Sound statistics: time step vs perturbation", fontsize=10)
78ax.legend(fontsize=8)
79ax.grid(alpha=0.3)
80fig.tight_layout()
81fig.savefig(out / "bands.png", dpi=150)
83# 5. aliasing
84fig, ax = plt.subplots(figsize=(6.5, 4))
85for a in d["aliasing"]:
86 ax.plot([b["f"] for b in a["bands"]], [b["db"] for b in a["bands"]], "o-",
87 label=f"{a['amp'] * 1000:g} mm: {a['fractionAbove24k'] * 100:.1f}% of power above 24 kHz")
88ax.axvline(24000, color="k", lw=0.8, ls="--")
89ax.text(25000, -5, "24 kHz", fontsize=8)
90ax.set_xscale("log")
91ax.set_xlabel("octave band centre (Hz)")
92ax.set_ylabel("pickup level (dB re loudest band)")
93ax.set_title("Pickup at 768 kHz: content above the 48 kHz Nyquist limit", fontsize=10)
94ax.legend(fontsize=8)
95ax.grid(alpha=0.3)
96fig.tight_layout()
97fig.savefig(out / "aliasing.png", dpi=150)
98print("wrote", out)