"""
FIGURE gr_tests — visual reproduction of win_gr_classical.py.

Main panel: grouped bars for the four classic GR tests (Mercury perihelion,
solar-limb deflection, Shapiro factor, GW speed ratio) showing DFD = GR = data
overlapping (PPN gamma=beta=1).

Inset panel: 2PN light-bending coefficient 4pi (DFD) vs 15pi/4 (GR) and the
+0.73 uas difference -- the falsifiable DFD edge.

All numbers are computed exactly as in win_gr_classical.py.
"""
import sys
sys.path.insert(0, ".")

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

# ----------------------------------------------------------------------
# Constants and PPN parameters (verbatim from win_gr_classical.py)
# ----------------------------------------------------------------------
ARCSEC = 206264.806        # arcsec per radian
m_sun = 1.476625           # GM_sun/c^2 in km
R_sun = 695700.0           # km

gGR, bGR = 1.0, 1.0
gDFD, bDFD = 1.0, 1.0      # DFD optical metric n=e^psi -> gamma=beta=1 (corpus theorem)

mb = m_sun / R_sun

# (1) Mercury perihelion: (2+2gamma-beta)/3 x 42.98''/cy
base_perih = 42.98
f = lambda g, b: (2 + 2*g - b)/3
mercury_GR = f(gGR, bGR) * base_perih
mercury_DFD = f(gDFD, bDFD) * base_perih
mercury_obs = 42.98        # observed (Mercury perihelion residual)

# (2) Light deflection at solar limb: (1+gamma)/2 x 4m/b
defl1 = lambda g: (1+g)/2 * 4*mb*ARCSEC
defl_GR = defl1(gGR)
defl_DFD = defl1(gDFD)
defl_obs = defl1(1.0)      # VLBI gamma-1 = (-0.8+/-1.2)e-4 -> consistent with 1.7512''

# (3) Shapiro time-delay factor: (1+gamma)/2
shapiro_GR = (1+gGR)/2
shapiro_DFD = (1+gDFD)/2
shapiro_obs = 1.0          # Cassini gamma-1 = (2.1+/-2.3)e-5

# (4) GW speed ratio c_T/c
gw_GR = 1.0
gw_DFD = 1.0               # no extra tensor d.o.f.
gw_obs = 1.0               # GW170817: |c_T/c - 1| < 1e-15

# 2PN light bending
c2_GR, c2_DFD = 15*np.pi/4, 4*np.pi
d2 = lambda c2: c2*mb**2*ARCSEC*1e6     # micro-arcsec
uas_GR = d2(c2_GR)
uas_DFD = d2(c2_DFD)
uas_diff = uas_DFD - uas_GR
ratio = c2_DFD / c2_GR

# ----------------------------------------------------------------------
# Colors
# ----------------------------------------------------------------------
C_DATA = "#444444"
C_GR   = "#1f77b4"
C_DFD  = "#d62728"

# ----------------------------------------------------------------------
# Figure layout: main grouped-bar panel + 2PN inset panel
# ----------------------------------------------------------------------
fig = plt.figure(figsize=(13, 6.4))
gs = fig.add_gridspec(1, 3, width_ratios=[2.05, 0.04, 1.0], wspace=0.28)
ax = fig.add_subplot(gs[0, 0])
ax2 = fig.add_subplot(gs[0, 2])

# --- Main panel: 4 tests, each normalized to its own observed value so they
#     share one axis. Bars sit at ratio = value/observed (all = 1.0). ---
tests = [
    ("Mercury\nperihelion", mercury_obs, mercury_GR, mercury_DFD, "42.98 ''/cy"),
    ("Solar-limb\ndeflection", defl_obs, defl_GR, defl_DFD, "1.7512 ''"),
    ("Shapiro\nfactor", shapiro_obs, shapiro_GR, shapiro_DFD, "1.0000"),
    ("GW speed\n$c_T/c$", gw_obs, gw_GR, gw_DFD, "1.0 ($<10^{-15}$)"),
]

x = np.arange(len(tests))
w = 0.26

data_vals = [t[1]/t[1] for t in tests]   # all 1.0 (observed normalization)
gr_vals   = [t[2]/t[1] for t in tests]
dfd_vals  = [t[3]/t[1] for t in tests]

b1 = ax.bar(x - w, data_vals, w, color=C_DATA, label="Observation", zorder=3,
            edgecolor="white", linewidth=0.7)
b2 = ax.bar(x,     gr_vals,   w, color=C_GR,   label="GR  ($\\gamma=\\beta=1$)", zorder=3,
            edgecolor="white", linewidth=0.7)
b3 = ax.bar(x + w, dfd_vals,  w, color=C_DFD,  label="DFD ($\\gamma=\\beta=1$)", zorder=3,
            edgecolor="white", linewidth=0.7, hatch="///")

ax.axhline(1.0, color="black", lw=1.0, ls="--", zorder=2, alpha=0.6)

# Annotate the actual physical value on each test group
for i, (name, obs, gr, dfd, label) in enumerate(tests):
    ax.text(i, 1.085, label, ha="center", va="bottom", fontsize=10,
            fontweight="bold", color="#222222")
    ax.text(i, 0.50, "DFD = GR\n= data", ha="center", va="center", fontsize=8.5,
            color="white", fontweight="bold", zorder=4)

ax.set_xticks(x)
ax.set_xticklabels([t[0] for t in tests], fontsize=10.5)
ax.set_ylabel("value / observed", fontsize=11.5)
ax.set_ylim(0, 1.25)
ax.set_title("Four classic GR tests — DFD reproduces all four (1PN: $\\gamma=\\beta=1$)",
             fontsize=12.5, fontweight="bold", pad=10)
ax.legend(loc="lower center", ncol=3, fontsize=9.5, framealpha=0.95,
          bbox_to_anchor=(0.5, -0.205))
ax.grid(axis="y", alpha=0.25, zorder=0)
ax.set_axisbelow(True)

# --- 2PN inset panel: light-bending coefficient DFD vs GR (in uas) ---
bx = [0, 1]
bvals = [uas_GR, uas_DFD]
bcols = [C_GR, C_DFD]
bars = ax2.bar(bx, bvals, width=0.6, color=bcols, edgecolor="white", linewidth=0.8,
               zorder=3)
bars[1].set_hatch("///")

ax2.text(0, uas_GR + 0.18, f"{uas_GR:.2f} uas", ha="center", va="bottom",
         fontsize=10.5, fontweight="bold", color=C_GR)
ax2.text(1, uas_DFD + 0.18, f"{uas_DFD:.2f} uas", ha="center", va="bottom",
         fontsize=10.5, fontweight="bold", color=C_DFD)

ax2.set_xticks(bx)
ax2.set_xticklabels(["GR\n$c_2=15\\pi/4$\n$=11.781$",
                     "DFD\n$c_2=4\\pi$\n$=12.566$"], fontsize=10)
ax2.set_ylabel("2PN solar-limb bending  [$\\mu$as]", fontsize=11)
ax2.set_ylim(0, max(bvals) * 1.28)
ax2.set_title("2PN light bending\n(DFD's sharper prediction)",
              fontsize=12, fontweight="bold", pad=10)
ax2.grid(axis="y", alpha=0.25, zorder=0)
ax2.set_axisbelow(True)

# Difference annotation bracket between the two bars
ytop = max(bvals) * 1.10
ax2.annotate("", xy=(1, ytop), xytext=(0, ytop),
             arrowprops=dict(arrowstyle="<->", color="black", lw=1.3))
ax2.text(0.5, ytop + 0.05,
         f"$\\Delta = +{uas_diff:.2f}$ $\\mu$as\n(ratio {ratio:.4f} $= 16/15$, +6.7%)",
         ha="center", va="bottom", fontsize=9.8, fontweight="bold",
         bbox=dict(boxstyle="round,pad=0.3", fc="#fff7e0", ec="#caa84a", lw=0.8))

fig.suptitle("DFD vs GR — classical gravity tests (PPN $\\gamma=\\beta=1$ reproduces all four; "
             "2PN gives a falsifiable +0.73 $\\mu$as edge)",
             fontsize=12.0, y=0.995)

fig.text(0.012, 0.012,
         "Source: win_gr_classical.py  |  $m_\\odot=GM_\\odot/c^2=1.4766$ km, "
         "$R_\\odot=695700$ km  |  DFD optical metric $n=e^\\psi$ (corpus theorem)",
         fontsize=7.8, color="#666666", ha="left")

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

# Echo the confirmed numbers
print("\nCONFIRMED NUMBERS:")
print(f"  Mercury:   GR={mercury_GR:.3f}  DFD={mercury_DFD:.3f}  ''/cy")
print(f"  Deflection:GR={defl_GR:.4f}  DFD={defl_DFD:.4f}  ''")
print(f"  Shapiro:   GR={shapiro_GR:.4f}  DFD={shapiro_DFD:.4f}")
print(f"  GW speed:  GR={gw_GR:.1f}  DFD={gw_DFD:.1f}")
print(f"  2PN: GR={c2_GR:.3f}->{uas_GR:.2f}uas  DFD={c2_DFD:.3f}->{uas_DFD:.2f}uas  "
      f"ratio={ratio:.4f}  diff={uas_diff:.2f}uas")
