"""Our analysis for https://factchecked.science/galactic-bending-waves/ : re-derive Fig. 4 of
Griv, E., and Yahalom, A. (2026), Astrophysics and Space Science 371, 118, https://doi.org/10.1007/s10509-026-04651-8
(CC BY 4.0), from its eqs. (48) and (49), and compare with our reading of the figure (fig4_readings.csv, hosted beside
this script). Standard library only.

The paper's eqs. (48)-(49), outer cyclotron resonance l = +1, with G rho0 h_z^2 = c_z^2/(4 pi) (h_z^2 ~ c_z^2/mu^2 and
mu^2 ~ 4 pi G rho0, p. 11):
  Re w*  = kappa * sqrt(1 - (8/pi^2) (c_z/c_r)^2 e^-x I1(x))
  Im w*  = (c_z^2 kappa^3 / (pi m Re w* |D| c_r^4)) * J1(chi')^2 * exp(-a/2),  a = (kappa - Re w*)/(m |D| c_r^2) > 0
with x = k^2 c_r^2 / kappa^2. The paper prints chi' = (k c_r/kappa) sqrt(a) = sqrt(x) sqrt(a) (after its eq. 41);
we also try chi' = x sqrt(a), and both J1 exactly and as the paper's approximation J1^2 ~ chi'^2/4 (after eq. 47).
Drift: the caption does not say which value Fig. 4 uses. The text gives r_epic <~ 1 kpc and |L_Omega| ~ 3 kpc
(p. 8), so |D| c_r^2 ~ r_epic^2/(r0 L_Omega) Omega = 0.0407 Omega (2D c_r^2 = -0.081 Omega; p. 11 rounds it to
"~ -0.1 Omega"); eq. (24) for the caption's flat rotation curve at r0 = 8.2 kpc gives |D| c_r^2 = (c_r/(kappa r0))^2
Omega = 0.0114 Omega instead.
Usage: python3 -I fig4_rederive.py [fig4_readings.csv]
"""
import csv
import math
import os
import sys

OMEGA, KAPPA, C_R, C_Z, R0 = 28.29, 40.01, 35.0, 20.0, 8.2      # caption of Fig. 4 (km/s, km/s/kpc, kpc)
K_O = KAPPA / OMEGA                                             # kappa / Omega
D_PAPER = 1.0 ** 2 / (R0 * 3.0)                                 # |D| c_r^2 / Omega from r_epic = 1 kpc, |L_Omega| = 3 kpc
D_FLAT = (C_R / (KAPPA * R0)) ** 2                              # |D| c_r^2 / Omega from eq. (24), flat curve


def i1(x, n=60):
    return sum((x / 2) ** (2 * j + 1) / (math.factorial(j) * math.factorial(j + 1)) for j in range(n))


def j1(z, n=80):
    return sum((-1) ** j * (z / 2) ** (2 * j + 1) / (math.factorial(j) * math.factorial(j + 1)) for j in range(n))


def re_w(x):
    """Re(omega*)/Omega, eq. (48)."""
    return K_O * math.sqrt(1 - (8 / math.pi ** 2) * (C_Z / C_R) ** 2 * math.exp(-x) * i1(x))


def im_w(x, m, d, chi="x", bessel="exact"):
    """Im(omega*)/Omega, eq. (49); d = |D| c_r^2 / Omega; chi = 'x' (x sqrt a) or 'printed' (sqrt(x) sqrt a)."""
    re = re_w(x)
    a = (K_O - re) / (m * d)
    c = (x if chi == "x" else math.sqrt(x)) * math.sqrt(a)
    jj = j1(c) ** 2 if bessel == "exact" else c * c / 4
    return (C_Z / C_R) ** 2 * K_O ** 3 * jj * math.exp(-a / 2) / (math.pi * m * re * d)


def readings(path):
    out = {}
    with open(path, encoding="utf-8") as f:
        rows = [ln for ln in f if not ln.startswith("#")]
    for r in csv.DictReader(rows):
        for m in (1, 2, 3, 4):
            v = r[f"im_m{m}"]
            if v != "":
                out[(round(float(r["x"]), 2), m)] = float(v)
    return out


path = sys.argv[1] if len(sys.argv) > 1 else os.path.join(os.path.dirname(os.path.abspath(__file__)), "fig4_readings.csv")
read = readings(path) if os.path.isfile(path) else {}

print(f"Drift |D| c_r^2: {D_PAPER:.4f} Omega from the paper's p. 8 values (2D c_r^2 = {-2 * D_PAPER:+.4f} Omega); "
      f"{D_FLAT:.4f} Omega from eq. (24) for the flat curve (2D c_r^2 = {-2 * D_FLAT:+.4f} Omega)")
print(f"Re(omega*)/Omega, eq. (48): x = 0 {re_w(0):.3f}, 0.5 {re_w(0.5):.3f}, 1 {re_w(1.0):.3f}")

print("\n1. How well each form reproduces our reading of Fig. 4 (m = 1-4, x = 0.20-1.00, every 0.02)")
pts = [(x, m) for (x, m) in read if 0.2 <= x <= 1.0]
best = None
for chi in ("x", "printed"):
    for bessel in ("exact", "approx"):
        for label, d in (("p. 8 drift", D_PAPER), ("eq. (24) drift", D_FLAT)):
            if not pts:
                break
            res = [im_w(x, m, d, chi, bessel) - read[(x, m)] for (x, m) in pts]
            rms = math.sqrt(sum(r * r for r in res) / len(res))
            print(f"   chi' = {'x sqrt(a)' if chi == 'x' else 'sqrt(x) sqrt(a) (printed)':<28} J1 {bessel:<6} {label:<15} rms {rms:.4f} Omega over {len(res)} points")
            if best is None or rms < best[0]:
                best = (rms, chi, bessel, label)
        # best-fitting drift for this form, by a scan
        if pts:
            scan = []
            for i in range(200, 1601):
                d = i / 20000                                  # 0.0100 to 0.0800 Omega
                res = [im_w(x, m, d, chi, bessel) - read[(x, m)] for (x, m) in pts]
                scan.append((math.sqrt(sum(r * r for r in res) / len(res)), d))
            rms, d = min(scan)
            print(f"      best drift for this form: |D| c_r^2 = {d:.4f} Omega (2D c_r^2 = {-2 * d:+.4f} Omega), rms {rms:.4f} Omega")
if best:
    print(f"   closest: chi' {best[1]}, J1 {best[2]}, {best[3]} (rms {best[0]:.4f})")

print("\n2. The re-derived curves (chi' = x sqrt(a), J1 exact, p. 8 drift) beside our reading, Im(omega*)/Omega")
for x in (0.1, 0.2, 0.4, 0.6, 0.8, 1.0):
    cells = []
    for m in (1, 2, 3, 4):
        r = read.get((round(x, 2), m))
        cells.append(f"m={m} {im_w(x, m, D_PAPER):.3f}" + (f" (read {r:.3f})" if r is not None else ""))
    print(f"   x = {x:.1f}: " + "; ".join(cells))

inv_omega_myr = 1 / (OMEGA / 3.085677581e16 * 365.25 * 86400) / 1e6
print("\n3. Growth of one-armed waves (m = 1), Im(omega*)/Omega, at the scales on the page"
      f"\n   (e-folding time in Gyr at the caption's Omega, 1/Omega = {inv_omega_myr:.1f} Myr, in brackets)")
L = C_R / KAPPA
cases = [("1/k = 3 kpc", (L / 3.0) ** 2), ("1/k = 2 kpc", (L / 2.0) ** 2),
         ("2 pi/k = 3 kpc", (2 * math.pi * L / 3.0) ** 2), ("2 pi/k = 2.5 kpc", (2 * math.pi * L / 2.5) ** 2),
         ("2 pi/k = 2 kpc", (2 * math.pi * L / 2.0) ** 2)]
for label, x in cases:
    vals = []
    for dl, d in (("p. 8 drift", D_PAPER), ("eq. (24) drift", D_FLAT)):
        for chi in ("x", "printed"):
            g = im_w(x, 1, d, chi)
            vals.append(f"{dl}, chi' {'x sqrt a' if chi == 'x' else 'printed'}: {g:.3f} [{inv_omega_myr / g / 1000:.2f}]")
    note = "" if x <= 1 else "   [outside the theory's x <~ 1 (p. 10): the formula's approximations no longer hold]"
    print(f"   {label:<17} x = {x:.3f}: " + "; ".join(vals) + note)
for d, dl in ((D_PAPER, "p. 8 drift"), (D_FLAT, "eq. (24) drift")):
    lo, hi = 0.0, 1.0                                                  # where m = 1 reaches 0.1 Omega (chi' = x sqrt a)
    for _ in range(60):
        mid = (lo + hi) / 2
        lo, hi = (mid, hi) if im_w(mid, 1, d) < 0.1 else (lo, mid)
    print(f"   m = 1 reaches 0.1 Omega at x = {lo:.3f} with the {dl}: 1/k = {L / math.sqrt(lo):.2f} kpc, 2 pi/k = {2 * math.pi * L / math.sqrt(lo):.1f} kpc")

print(f"\n4. At the caption's Omega, 1/Omega = {inv_omega_myr:.1f} Myr: one e-fold at 0.1 Omega takes {inv_omega_myr / 0.1 / 1000:.2f} Gyr, "
      f"at 0.02 Omega {inv_omega_myr / 0.02 / 1000:.2f} Gyr (the paper's round figure, t_grow ~ 1 Gyr, takes Omega ~ 1e-8 per yr)")

if len(sys.argv) > 2:                                            # optional: write the re-derived curves for the chart
    with open(sys.argv[2], "w", newline="", encoding="utf-8") as f:
        f.write("# Im(omega*)/Omega from eq. (49) of Griv & Yahalom (2026) with chi' = x sqrt(a), J1 exact, and the p. 8 drift "
                f"(|D| c_r^2 = {D_PAPER:.4f} Omega); m1_flat uses the eq. (24) drift for the caption's flat curve ({D_FLAT:.4f} Omega). Our analysis.\n")
        w = csv.writer(f)
        w.writerow(["x", "m1", "m2", "m3", "m4", "m1_flat"])
        for i in range(0, 101):
            x = i / 100
            if x == 0:
                w.writerow(["0.00", "0", "0", "0", "0", "0"])
                continue
            w.writerow([f"{x:.2f}"] + [f"{im_w(x, m, D_PAPER):.4f}" for m in (1, 2, 3, 4)] + [f"{im_w(x, 1, D_FLAT):.4f}"])
    print("wrote", sys.argv[2])