1"""Figures from results/*.npz (written by experiments.py). Run: python plots.py"""
3from __future__ import annotations
5import os
7import matplotlib
9matplotlib.use("Agg")
10import matplotlib.pyplot as plt # noqa: E402
11import numpy as np # noqa: E402
13HERE = os.path.dirname(os.path.abspath(__file__))
14RES = os.path.join(HERE, "results")
15CHAOS = 0.02 # finite-time Lyapunov exponent above which an orbit counts as chaotic
18def sections():
19 d = np.load(os.path.join(RES, "sections.npz"))
20 energies, w2s = d["energies"], d["w2s"]
21 fig, ax = plt.subplots(len(w2s), len(energies), figsize=(3 * len(energies), 3.1 * len(w2s)))
22 for i, w2 in enumerate(w2s):
23 for j, E in enumerate(energies):
24 k = f"w2_{w2:g}_E_{E:g}"
25 lam, owner, pts = d[k + "_lyap"], d[k + "_owner"], d[k + "_pts"]
26 a = ax[i, j]
27 chaotic = lam[owner] > CHAOS
28 reg = ~chaotic
29 a.scatter(pts[reg, 0], pts[reg, 1], s=0.15, c=owner[reg] % 10, cmap="tab10", vmin=0, vmax=9, rasterized=True)
30 a.scatter(pts[chaotic, 0], pts[chaotic, 1], s=0.15, c="k", alpha=0.5, rasterized=True)
31 a.set_title(f"w2 = {w2:g}, E = {E:g}\n{np.mean(lam > CHAOS):.0%} of orbits chaotic", fontsize=9)
32 a.set_xticks([])
33 a.set_yticks([])
34 if j == 0:
35 a.set_ylabel("p_x (section y = 0, p_y > 0)")
36 if i == len(w2s) - 1:
37 a.set_xlabel("x")
38 fig.suptitle("Poincare sections: colored curves are regular orbits (tori), black dots chaotic orbits", y=1.0)
39 fig.tight_layout()
40 fig.savefig(os.path.join(RES, "sections.png"), dpi=110)
43def lyapunov():
44 d = np.load(os.path.join(RES, "lyapunov.npz"))
45 energies, w2s = d["energies"], d["w2s"]
46 fig, ax = plt.subplots(1, 2, figsize=(12, 4.5))
47 lines = []
48 for w2 in w2s:
49 lams = [d[f"w2_{w2:g}_E{i}"] for i in range(len(energies))]
50 frac = [np.mean(l > CHAOS) for l in lams]
51 med = [np.median(l[l > CHAOS]) if np.any(l > CHAOS) else np.nan for l in lams]
52 ax[0].semilogx(energies, frac, "o-", label=f"w2 = {w2:g}")
53 ax[1].loglog(energies, med, "o-", label=f"w2 = {w2:g}")
54 lines.append((w2, energies, frac, med))
55 ax[0].set(xlabel="energy E (beta = 1)", ylabel="fraction of orbits chaotic", title="Onset of chaos")
56 ax[1].set(xlabel="energy E", ylabel="median Lyapunov exponent of chaotic orbits",
57 title="Rate of separation in the chaotic regime")
58 for a in ax:
59 a.grid(True, which="both", alpha=0.3)
60 a.legend()
61 fig.tight_layout()
62 fig.savefig(os.path.join(RES, "lyapunov.png"), dpi=110)
63 return lines
66def damped():
67 d = np.load(os.path.join(RES, "damped.npz"))
68 E0s = sorted({float(k.split("_")[1]) for k in d.files if k.startswith("E0_")})
69 fig, ax = plt.subplots(3, len(E0s), figsize=(5.5 * len(E0s), 10), sharex=True)
70 summary = []
71 for j, E0 in enumerate(E0s):
72 k = f"E0_{E0:g}"
73 t, Ex, Ey, E, logsum = d[k + "_t"], d[k + "_Ex"], d[k + "_Ey"], d[k + "_E"], d[k + "_logsum"]
74 frac = Ex / (Ex + Ey)
75 ax[0, j].semilogy(t, E[:, 0], "k")
76 ax[0, j].set(title=f"strike energy E0 = {E0:g}", ylabel="energy (member 0)")
77 ax[1, j].plot(t, frac, lw=0.6)
78 ax[1, j].set(ylabel="share of energy in the x oscillator\n(16 strikes differing by 1e-10)", ylim=(0, 1))
79 # local Lyapunov exponent: slope of the log-growth sum over windows of 20 time units
80 w = max(1, int(round(20 / (t[1] - t[0]))))
81 loc = (logsum[w:, 0] - logsum[:-w, 0]) / (t[w:] - t[:-w])
82 ax[2, j].plot(t[w:], loc, "k", lw=0.8)
83 ax[2, j].axhline(0, color="0.6", lw=0.5)
84 ax[2, j].set(xlabel="time", ylabel="local Lyapunov exponent (20-unit windows)")
85 late = t > t[-1] * 0.8
86 summary.append((E0, float(frac[late].mean(0).std()), float(frac[late].mean(0).min()), float(frac[late].mean(0).max())))
87 fig.tight_layout()
88 fig.savefig(os.path.join(RES, "damped.png"), dpi=110)
89 return summary
92if __name__ == "__main__":
93 for name, f in (("sections", sections), ("lyapunov", lyapunov), ("damped", damped)):
94 if os.path.exists(os.path.join(RES, name + ".npz")):
95 out = f()
96 if name == "lyapunov":
97 for w2, energies, frac, med in out:
98 print(f"w2 = {w2:g}:", " ".join(f"E {E:.2g}: {fr:.2f}" for E, fr in zip(energies, frac)))
99 if name == "damped":
100 for E0, sd, lo, hi in out:
101 print(f"damped E0 {E0:g}: late share of energy in x across strikes: std {sd:.3f}, range {lo:.3f} to {hi:.3f}")
102 print("wrote", os.path.join(RES, name + ".png"))