"""
FIGURE: mond — three-panel reproduction of win_mond.py.
(a) a0_DFD = 2*sqrt(alpha)*c*H0  vs SPARC 1.20e-10 with error band; annotate -0.28%.
(b) Baryonic Tully-Fisher: M_b = V^4/(G a0), with observed 47-63 Msun/(km/s)^4 band, mark V=100,200.
(c) Interpolation mu(x)=x/(1+x) and g/g_N = 1/mu vs x=g/a0, log axes, Newtonian & deep-MOND limits.
Confirms a0 = 1.197e-10 and slope = 63 Msun/(km/s)^4.
"""
import sys
sys.path.insert(0, ".")

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

# ---- exact constants from win_mond.py / dfd_params.py ----
ALPHA = 1.0 / 137.035999084
C = 2.99792458e8           # m/s
MPC = 3.0856775815e22      # m
MSUN = 1.98892e30          # kg
G = 6.67430e-11            # m^3/kg/s^2
H0_DFD = 72.09             # km/s/Mpc, DERIVED


def a0_dfd(H0_kms_Mpc):
    H0 = H0_kms_Mpc * 1000.0 / MPC          # s^-1
    return 2.0 * np.sqrt(ALPHA) * C * H0    # m/s^2


a0 = a0_dfd(H0_DFD)                          # 1.1966e-10
a0_obs, a0_err = 1.20e-10, 0.24e-10         # SPARC (Lelli+2017)
pct = 100.0 * (a0 - a0_obs) / a0_obs        # -0.28%

A_dfd = 1.0 / (G * a0)                       # kg per (m/s)^4
A_dfd_sun_kms = A_dfd * (1000.0 ** 4) / MSUN # Msun per (km/s)^4  -> 63.0

print(f"a0_DFD = {a0:.4e} m/s^2  ({pct:+.2f}% vs SPARC, {abs(a0-a0_obs)/a0_err:.2f} sigma_sys)")
print(f"BTFR slope 1/(G a0) = {A_dfd_sun_kms:.1f} Msun/(km/s)^4")

# ---- style ----
DFD_C = "#c0392b"      # DFD red
OBS_C = "#2c3e50"      # observation dark slate
BAND_C = "#27ae60"     # observed band green
NEWT_C = "#2980b9"     # Newtonian blue

plt.rcParams.update({"font.size": 10, "axes.grid": True,
                     "grid.alpha": 0.3, "axes.axisbelow": True})

fig, (ax0, ax1, ax2) = plt.subplots(1, 3, figsize=(15, 5))
fig.suptitle(r"DFD subsumes MOND: $a_0 = 2\sqrt{\alpha}\,c\,H_0$ derived (no fit)",
             fontsize=14, fontweight="bold", y=1.00)

# ============================================================
# (a) a0 comparison: DFD derived vs SPARC observed
# ============================================================
labels = ["SPARC\n(observed)", "DFD\n(derived)"]
vals = [a0_obs * 1e10, a0 * 1e10]
errs = [a0_err * 1e10, 0.0]
colors = [OBS_C, DFD_C]
x = np.arange(2)

# observed systematic error band across whole panel
ax0.axhspan((a0_obs - a0_err) * 1e10, (a0_obs + a0_err) * 1e10,
            color=BAND_C, alpha=0.15, zorder=0,
            label=r"SPARC $\pm1\sigma_{\rm sys}$")
ax0.axhline(a0_obs * 1e10, color=BAND_C, ls="--", lw=1.2, zorder=1)

ax0.bar(x, vals, yerr=errs, width=0.5, color=colors, alpha=0.85,
        edgecolor="black", linewidth=1.0, capsize=8, zorder=3)
for xi, v in zip(x, vals):
    ax0.text(xi, v + 0.02, f"{v:.3f}", ha="center", va="bottom",
             fontsize=11, fontweight="bold")

ax0.set_ylim(0, 1.55)
ax0.set_xticks(x)
ax0.set_xticklabels(labels)
ax0.set_ylabel(r"$a_0\ \ [10^{-10}\ \mathrm{m\,s^{-2}}]$")
ax0.set_title("(a)  Universal acceleration scale")
ax0.legend(loc="upper left", fontsize=8)
# annotate the headline agreement
ax0.annotate(f"$a_0^{{\\rm DFD}}={a0*1e10:.3f}\\times10^{{-10}}$\n"
             f"agreement: {pct:+.2f}%\n({abs(a0-a0_obs)/a0_err:.2f}$\\,\\sigma_{{\\rm sys}}$)",
             xy=(1, a0 * 1e10), xytext=(0.40, 0.55),
             textcoords="axes fraction",
             bbox=dict(boxstyle="round,pad=0.4", fc="#fff3cd", ec=DFD_C, lw=1.3),
             fontsize=9.5, ha="left",
             arrowprops=dict(arrowstyle="->", color=DFD_C, lw=1.3))

# ============================================================
# (b) Baryonic Tully-Fisher relation
# ============================================================
V = np.logspace(np.log10(20), np.log10(400), 300)   # km/s
Mb = (V * 1000.0) ** 4 / (G * a0) / MSUN             # Msun, DFD slope

# observed normalization band 47-63 Msun/(km/s)^4
Mb_lo = 47.0 * V ** 4
Mb_hi = 63.0 * V ** 4
ax1.fill_between(V, Mb_lo, Mb_hi, color=BAND_C, alpha=0.18, zorder=0,
                 label=r"observed $47$–$63\ M_\odot/(\mathrm{km/s})^4$")

ax1.loglog(V, Mb, color=DFD_C, lw=2.4, zorder=3,
           label=rf"DFD: $M_b=V^4/(G a_0)$, slope $={A_dfd_sun_kms:.0f}$")

# mark V = 100 and 200 km/s
for Vmark, tag in ((100.0, "dwarf-ish"), (200.0, "MW-class")):
    Mbm = (Vmark * 1000.0) ** 4 / (G * a0) / MSUN
    ax1.scatter([Vmark], [Mbm], s=80, color="black", zorder=5,
                edgecolor="white", linewidth=1.2)
    ax1.annotate(f"$V={Vmark:.0f}$\n{Mbm:.2e} $M_\\odot$\n({tag})",
                 xy=(Vmark, Mbm), xytext=(Vmark * 1.08, Mbm * 0.22),
                 fontsize=8.5,
                 arrowprops=dict(arrowstyle="->", color="black", lw=1.0))

ax1.set_xlabel(r"$V_{\rm flat}\ \ [\mathrm{km/s}]$")
ax1.set_ylabel(r"$M_b\ \ [M_\odot]$")
ax1.set_title("(b)  Baryonic Tully–Fisher relation")
ax1.set_xlim(20, 400)
ax1.legend(loc="upper left", fontsize=8.5)
ax1.text(0.97, 0.06,
         rf"slope $1/(G a_0) = {A_dfd_sun_kms:.0f}\ M_\odot/(\mathrm{{km/s}})^4$"
         "\n[DERIVED, no fit]",
         transform=ax1.transAxes, ha="right", va="bottom", fontsize=8.5,
         bbox=dict(boxstyle="round,pad=0.35", fc="#fff3cd", ec=DFD_C, lw=1.1))

# ============================================================
# (c) Interpolation function mu(x) and g/g_N = 1/mu
# ============================================================
xv = np.logspace(-3, 3, 400)        # x = g/a0
mu = xv / (1.0 + xv)                 # DFD optical-Poisson mu-law
g_over_gN = 1.0 / mu                 # boost = 1/mu

ax2.loglog(xv, mu, color=DFD_C, lw=2.4, label=r"$\mu(x)=x/(1+x)$")
ax2.loglog(xv, g_over_gN, color=NEWT_C, lw=2.4, ls="-",
           label=r"$g/g_N = 1/\mu$")

# deep-MOND limit: 1/mu -> 1/x (large boost); Newtonian: 1/mu -> 1
ax2.loglog(xv[xv < 1], 1.0 / xv[xv < 1], color=NEWT_C, lw=1.0, ls=":",
           alpha=0.8, label=r"deep-MOND $g/g_N\!\to\!1/x$")
ax2.axhline(1.0, color="gray", lw=1.0, ls="--", alpha=0.7)

# regime shading
ax2.axvspan(1e-3, 0.1, color=BAND_C, alpha=0.08)
ax2.axvspan(10, 1e3, color=NEWT_C, alpha=0.08)

# mark the three sample points from win_mond.py
for xs, lab in ((0.01, "deep-MOND\n(outskirts)"),
                (1.0, "transition"),
                (100.0, "Newtonian\n(solar sys.)")):
    mus = xs / (1.0 + xs)
    ax2.scatter([xs], [1.0 / mus], s=55, color="black", zorder=5,
                edgecolor="white", linewidth=1.0)

ax2.text(0.03, 0.90, "deep-MOND\n$g/g_N\\simeq1/x$", transform=ax2.transAxes,
         fontsize=8.5, color=BAND_C, fontweight="bold", va="top")
ax2.text(0.78, 0.18, "Newtonian\n$g/g_N\\to1$", transform=ax2.transAxes,
         fontsize=8.5, color=NEWT_C, fontweight="bold", ha="center")

ax2.set_xlabel(r"$x = g/a_0$")
ax2.set_ylabel(r"$\mu(x)$  and  $g/g_N = 1/\mu(x)$")
ax2.set_title("(c)  Interpolation: both MOND regimes")
ax2.set_xlim(1e-3, 1e3)
ax2.set_ylim(5e-3, 2e2)
ax2.legend(loc="upper center", fontsize=8.0, framealpha=0.9)

fig.tight_layout(rect=[0, 0, 1, 0.97])
out = "./figures/mond.png"
fig.savefig(out, dpi=150, bbox_inches="tight")
print("saved:", out)
