#!/usr/bin/env python3
# =============================================================================
# growth_regime_adjudication.py
#
# DFD growth-regime adjudication (Appendix GR companion).
#
# Resolves the AC (Hubble-EFE, x_bar=5.85, G_eff=1.02-1.17 G) vs
# AE (FRW deep-MOND, a_ext=0, x_bar=0) regime question and produces the
# SINGLE DFD sigma8/S8/fsigma8 prediction, honest to the independently verified tally:
#
#   - contradiction = scope-distinction with a real localized residual
#   - geff_derived  = FALSE  (the AQUAL G_eff is NOT the cosmological growth carrier)
#   - the cosmological linear growth of delta_m is carried by chi-matter
#     (Omega_chi h^2 = 0.12, Q ~ 1, GR-like growth), NOT by AC's 1.12 G
#     and NOT by AE's runaway deep-MOND self-field.
#   - fsigma8 fit is ACCEPTABLE, NOT superior to LCDM (chi2/N ~ 1.6 vs ~1.1).
#   - S8 = 0.755 lands on the lensing-LOW side of the S8 tension.
#
# Outputs:
#   * console table of x = g/a0(z) on R8 (proves deep-MOND on all linear scales)
#   * the AC envelope numbers (x_bar=5.85, G_perp, G_par, <Geff>)
#   * D(a), f(a), sigma8(z), S8 for the chi-matter (Q=1) growth law
#   * the fsigma8(z) table and chi2/N vs DESI/eBOSS, vs LCDM
#   * the integrated-sigma8 over-growth flag if AC's Q=1.118 is (wrongly)
#     applied over the whole history  (honesty-guard: shows why AC is NOT cosmological)
#   * figure  dfd_class/fig_growth_regime.png  (panel A + panel B)
#
# No fitted parameters. All cosmology numbers come from the frozen DFD corpus.
# =============================================================================

import numpy as np
from scipy.integrate import solve_ivp, quad
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

# -----------------------------------------------------------------------------
# 0. Frozen DFD / cosmology constants (corpus values, no refit)
# -----------------------------------------------------------------------------
alpha     = 1.0 / 137.036          # fine-structure constant
sqrt_a    = np.sqrt(alpha)
H0_kms    = 72.09                  # km/s/Mpc  (DFD-CLASS derived)
Om0       = 0.274                  # total matter (baryon + chi); DFD-CLASS background
Ob0       = 0.04304               # baryon
OL0       = 1.0 - Om0             # optical / cosmological-constant-like late-time DE (w=-1)
sigma8_0  = 0.790                 # corpus CMB-normalized central value
c_kms     = 299792.458            # km/s
R8_Mpc    = 8.0 / (H0_kms/100.0)   # 8 h^-1 Mpc in Mpc  (h = H0/100)

# a0 today (deep-MOND crossover) : a0 = 2 sqrt(alpha) c H0
# H0 in SI:
Mpc_m     = 3.085677581e22
H0_SI     = H0_kms * 1.0e3 / Mpc_m          # 1/s
c_SI      = 2.99792458e8
a0_today  = 2.0 * sqrt_a * c_SI * H0_SI     # m/s^2
print("="*78)
print("DFD GROWTH-REGIME ADJUDICATION")
print("="*78)
print(f"a0(today) = 2 sqrt(alpha) c H0 = {a0_today:.4e} m/s^2   (corpus 1.197e-10)")

def E(a):
    """H(a)/H0 for flat matter + DE (w=-1)."""
    return np.sqrt(Om0 * a**-3 + OL0)

def Om_a(a):
    return Om0 * a**-3 / E(a)**2

# -----------------------------------------------------------------------------
# 1. AC ENVELOPE NUMBERS  (Hubble-EFE hypothesis; redshift-independent)
# -----------------------------------------------------------------------------
xbar_AC = 1.0 / (2.0 * sqrt_a)              # = 5.85, the AC value
mu0     = xbar_AC / (1.0 + xbar_AC)
muln0   = xbar_AC / (1.0 + xbar_AC)**2
q_AC    = 1.0 / (1.0 + xbar_AC)             # = muln0/mu0
G_perp  = 1.0 / mu0                          # k perp g
G_par   = 1.0 / (mu0 + muln0)                # k || g
# angular mean  <Geff>/G = (1/mu0) * arctan(sqrt q)/sqrt q
G_mean  = (1.0 / mu0) * np.arctan(np.sqrt(q_AC)) / np.sqrt(q_AC)
print("\n--- AC Hubble-EFE envelope (admitted hypothesis, NOT action-derived) ---")
print(f"  x_bar = 1/(2 sqrt alpha)       = {xbar_AC:.4f}   (redshift-INDEPENDENT: cH cancels)")
print(f"  mu0   = {mu0:.4f}   mu_ln,0 = {muln0:.4f}")
print(f"  G_perp/G = {G_perp:.4f}   G_par/G = {G_par:.4f}   <Geff>/G = {G_mean:.4f}")

# -----------------------------------------------------------------------------
# 2. x = g/a0(z) ON R8  (proves linear cosmological scales are DEEP-MOND)
# -----------------------------------------------------------------------------
# Spherical top-hat on R8: g_N = (1/2) Om(a) H(a)^2 R8 delta   (physical accel)
# a0(z) = 2 sqrt(alpha) c H(z)
def x_on_R8(z, delta):
    a   = 1.0/(1.0+z)
    H   = H0_SI * E(a)                       # 1/s
    R8m = R8_Mpc * Mpc_m / a                 # physical radius (proper) at scale factor a
    # use comoving R8 in proper units at epoch a:
    gN  = 0.5 * Om0 * a**-3 * H0_SI**2 * (R8_Mpc*Mpc_m) * delta  # SI accel (comoving R8)
    a0z = 2.0 * sqrt_a * c_SI * H
    return gN / a0z

print("\n--- x = g_N/a0(z) on R8 (linear scales) ---")
print("   z      delta     x=g/a0        regime")
for z, d in [(0.0,0.01),(0.0,1.0),(0.5,1.0),(1.0,1.0),(2.0,1.0),(1090.0,1e-5)]:
    xx = x_on_R8(z,d)
    reg = "deep-MOND" if xx < 1 else "Newtonian"
    print(f"  {z:6.1f}  {d:8.2e}  {xx:10.3e}   {reg}")

# collapsed-halo and cluster checks (Newtonian / deep-MOND boundary)
G_grav = 6.674e-11
Msun   = 1.989e30
def x_point(M_Msun, r_kpc):
    r  = r_kpc * 3.0857e19
    gN = G_grav * M_Msun*Msun / r**2
    return gN / a0_today
print(f"  collapsed halo 1e12 Msun @ 30 kpc : x = {x_point(1e12,30):.3f}  (Newtonian-ish, transition)")
print(f"  cluster        1e14 Msun @ 1 Mpc  : x = {x_point(1e14,1000):.3f}  (still deep-MOND)")

# -----------------------------------------------------------------------------
# 3. GROWTH ODE  -- D(a) in scale factor.  delta'' form.
#    d2D/da2 + (3/a + E'/E) dD/da - (3/2) Om(a)/a^2 * Q * D = 0
# -----------------------------------------------------------------------------
def growth_D(Qfun, a_ini=1e-3):
    def Eprime_over_E(a):
        # dlnE/da
        num = -1.5*Om0*a**-4
        return num / E(a)**2
    def rhs(a, y):
        D, Dp = y
        coef1 = (3.0/a + Eprime_over_E(a))
        coef2 = 1.5 * Om_a(a) / a**2 * Qfun(a)
        return [Dp, -coef1*Dp + coef2*D]
    sol = solve_ivp(rhs, [a_ini, 1.0], [a_ini, 1.0],
                    t_eval=np.linspace(a_ini, 1.0, 4000),
                    rtol=1e-9, atol=1e-12, method="RK45")
    return sol.t, sol.y[0], sol.y[1]

# (a) chi-matter / GR-like growth : Q = 1  -- the DFD COSMOLOGICAL carrier
a_arr, D_chi, Dp_chi = growth_D(lambda a: 1.0)
D_chi_norm = D_chi / D_chi[-1]                     # D(a)/D(1)
f_chi = a_arr * Dp_chi / D_chi                     # f = dlnD/dlna

# (b) LCDM reference (Q=1 too, Om0=0.30) for f(0) sanity
def growth_D_om(Om_ref):
    OLr = 1.0-Om_ref
    Eref = lambda a: np.sqrt(Om_ref*a**-3+OLr)
    Omr  = lambda a: Om_ref*a**-3/Eref(a)**2
    def EpE(a): return (-1.5*Om_ref*a**-4)/Eref(a)**2
    def rhs(a,y):
        D,Dp=y
        return [Dp, -(3.0/a+EpE(a))*Dp + 1.5*Omr(a)/a**2*D]
    sol=solve_ivp(rhs,[1e-3,1.0],[1e-3,1.0],t_eval=np.linspace(1e-3,1,4000),
                  rtol=1e-9,atol=1e-12)
    return sol.t,sol.y[0],sol.y[1]
a_l,D_l,Dp_l=growth_D_om(0.30)
f_lcdm0 = (a_l*Dp_l/D_l)[-1]
print(f"\n--- growth check ---")
print(f"  f_LCDM(0) (Om=0.30) = {f_lcdm0:.3f}   (Om^0.55 = {0.30**0.55:.3f})")
print(f"  f_DFD(0)  (Om=0.274, chi Q=1) = {f_chi[-1]:.3f}")

# -----------------------------------------------------------------------------
# 4. sigma8(z), S8, fsigma8(z)   for the chi-matter (Q=1) carrier
# -----------------------------------------------------------------------------
from numpy import interp
def Dz(z):  return interp(1.0/(1.0+z), a_arr, D_chi_norm)
def fz(z):  return interp(1.0/(1.0+z), a_arr, f_chi)
def s8z(z): return sigma8_0 * Dz(z)
def fs8z(z):return fz(z) * s8z(z)

S8 = sigma8_0 * np.sqrt(Om0/0.3)
print(f"\n--- DFD growth prediction (chi-matter carrier, Q=1) ---")
print(f"  sigma8(0) = {sigma8_0:.3f}   f(0) = {f_chi[-1]:.3f}")
print(f"  S8 = sigma8 sqrt(Om/0.3) = {S8:.3f}")

# DESI/eBOSS fsigma8 compilation (z, fs8_obs, sigma)
desi = [
    (0.15, 0.53, 0.16),    # 6dF/SDSS MGS-ish low-z
    (0.38, 0.497,0.045),   # BOSS
    (0.51, 0.459,0.038),   # BOSS / DESI
    (0.70, 0.473,0.041),   # eBOSS LRG / DESI
    (0.85, 0.52, 0.10),    # eBOSS
    (1.48, 0.462,0.045),   # eBOSS QSO
]
print("\n   z     fs8_DFD   fs8_obs   sigma    pull(DFD)")
chi2_dfd=0.0; chi2_lcdm=0.0
# LCDM fs8 with sigma8=0.811, Om=0.31
s8_l0=0.811
a_l31,D_l31,Dp_l31=growth_D_om(0.31)
D_l31n=D_l31/D_l31[-1]; f_l31=a_l31*Dp_l31/D_l31
def fs8_lcdm(z):
    a=1.0/(1.0+z)
    return interp(a,a_l31,f_l31)*s8_l0*interp(a,a_l31,D_l31n)
fs8_dfd_tab=[]
for z,obs,sig in desi:
    m=fs8z(z); fs8_dfd_tab.append(m)
    ml=fs8_lcdm(z)
    pull=(m-obs)/sig
    chi2_dfd += pull**2
    chi2_lcdm += ((ml-obs)/sig)**2
    print(f"  {z:4.2f}   {m:6.3f}   {obs:6.3f}   {sig:5.3f}   {pull:+5.2f}")
N=len(desi)
print(f"\n  chi2/N  DFD  = {chi2_dfd/N:.3f}   (chi2={chi2_dfd:.2f}, N={N})")
print(f"  chi2/N  LCDM = {chi2_lcdm/N:.3f}   (chi2={chi2_lcdm:.2f}, N={N})")
print("  verdict: DFD fsigma8 fit ACCEPTABLE, NOT superior to LCDM.")

# -----------------------------------------------------------------------------
# 5. HONESTY-GUARD FLAG : if AC's Q=1.118 were applied over the WHOLE history,
#    integrated sigma8 over-grows (NOT AC's advertised +1-9% instantaneous f).
# -----------------------------------------------------------------------------
a_ac, D_ac, _ = growth_D(lambda a: G_mean)   # constant Q=1.118 from z~1000
sigma8_AC = sigma8_0 * (D_ac[-1]/D_ac[0]) / (D_chi[-1]/D_chi[0])
ratio = (D_ac[-1]/D_ac[0]) / (D_chi[-1]/D_chi[0])
print(f"\n--- honesty-guard: AC's constant Q={G_mean:.3f} integrated from z~1000 ---")
print(f"  D-ratio(AC/chi) = {ratio:.3f}  ->  sigma8 would be {sigma8_0*ratio:.3f}")
print(f"  i.e. +{100*(ratio-1):.0f}% integrated sigma8 (NOT the +1-9% instantaneous f AC advertises)")
print("  => AC's 1.118 G compounds; cannot be the cosmological growth law. Confirms chi carries growth.")

# -----------------------------------------------------------------------------
# 6. FIGURE
# -----------------------------------------------------------------------------
fig, (axA, axB) = plt.subplots(1, 2, figsize=(13, 5.2))

# Panel A : G_eff/G across regimes vs x = g/a0, with transition x=1
xs = np.logspace(-5, 3, 400)
mu_x = xs/(1.0+xs)
muln_x = xs/(1.0+xs)**2
Geff_perp = 1.0/mu_x
Geff_par  = 1.0/(mu_x+muln_x)
axA.loglog(xs, Geff_perp, color="#1f77b4", lw=2.2, label=r"$G_{\rm eff}/G$  ($k\perp g$)")
axA.loglog(xs, Geff_par,  color="#2ca02c", lw=2.2, ls="--", label=r"$G_{\rm eff}/G$  ($k\parallel g$)")
axA.axvline(1.0, color="0.4", ls=":", lw=1.5)
axA.text(1.15, 50, r"transition $x=g/a_0=1$", rotation=90, va="center", fontsize=9, color="0.3")
# mark AC and AE points
axA.axvline(xbar_AC, color="#d62728", ls="-.", lw=1.5, alpha=0.8)
axA.text(xbar_AC*1.1, 1.4, f"AC hypothesis\n"+r"$\bar x=5.85$ ($\langle G_{\rm eff}\rangle=1.12G$)",
         fontsize=8, color="#d62728", va="bottom")
# linear-scale band (deep-MOND): x ~ 2.5e-5 .. 2.5e-3
axA.axvspan(2.5e-5, 2.5e-3, color="orange", alpha=0.15)
axA.text(1.5e-4, 200, "linear cosmological\nscales (R8): deep-MOND\n"+r"but $\delta_m$ carried by $\chi$ ($Q\!\approx\!1$)",
         fontsize=8, color="#b35900", ha="left")
axA.axhline(1.0, color="k", lw=0.8, alpha=0.5)
axA.set_xlabel(r"$x = g/a_0(z)$")
axA.set_ylabel(r"$G_{\rm eff}/G$  (AQUAL $\psi$-operator)")
axA.set_title(r"A. AQUAL $G_{\rm eff}$ vs regime ($\psi$-operator is NOT the growth carrier)", fontsize=10)
axA.legend(loc="upper right", fontsize=8)
axA.set_ylim(0.8, 1e3)
axA.grid(True, which="both", alpha=0.2)

# Panel B : fsigma8(z) DFD vs DESI vs LCDM ; inset S8
zz = np.linspace(0.0, 1.6, 200)
axB.plot(zz, [fs8z(z) for z in zz], color="#d62728", lw=2.2,
         label=r"DFD ($\chi$-matter, $\sigma_8=0.790$, $\Omega_m=0.274$)")
axB.plot(zz, [fs8_lcdm(z) for z in zz], color="#1f77b4", lw=1.8, ls="--",
         label=r"$\Lambda$CDM ($\sigma_8=0.811$, $\Omega_m=0.31$)")
zo=[d[0] for d in desi]; yo=[d[1] for d in desi]; so=[d[2] for d in desi]
axB.errorbar(zo, yo, yerr=so, fmt="o", color="k", ms=5, capsize=3,
             label="DESI/eBOSS $f\\sigma_8$")
axB.set_xlabel(r"redshift $z$")
axB.set_ylabel(r"$f\sigma_8(z)$")
axB.set_title(r"B. $f\sigma_8$: DFD $\chi^2/N$=%.2f vs $\Lambda$CDM %.2f (acceptable, not superior)" % (chi2_dfd/N, chi2_lcdm/N), fontsize=10)
axB.legend(loc="lower center", fontsize=8)
axB.grid(True, alpha=0.2)
axB.set_ylim(0.30, 0.58)

# S8 inset
axins = axB.inset_axes([0.60, 0.62, 0.37, 0.34])
labels = ["DFD", "Planck\nLCDM", "DES-Y3", "KiDS"]
vals   = [S8, 0.834, 0.776, 0.759]
errs   = [0.0, 0.016, 0.017, 0.024]
cols   = ["#d62728", "#1f77b4", "#7f7f7f", "#7f7f7f"]
xpos=np.arange(len(labels))
axins.errorbar(xpos, vals, yerr=errs, fmt="s", ms=4, capsize=2,
               color="none", ecolor="k")
for xp,v,co in zip(xpos,vals,cols):
    axins.plot(xp,v,"s",ms=6,color=co)
axins.axhspan(0.74,0.78, color="orange", alpha=0.15)
axins.set_xticks(xpos); axins.set_xticklabels(labels, fontsize=6)
axins.set_ylabel(r"$S_8$", fontsize=7)
axins.tick_params(labelsize=6)
axins.set_title("S8: DFD on lensing-low side", fontsize=7)
axins.axhline(S8, color="#d62728", ls=":", lw=1)

plt.tight_layout()
outfig = "./fig_growth_regime.png"
plt.savefig(outfig, dpi=140)
print(f"\nFigure written: {outfig}")

# -----------------------------------------------------------------------------
# 7. MANIFEST SUMMARY
# -----------------------------------------------------------------------------
print("\n" + "="*78)
print("MANIFEST (honest to verified tally)")
print("="*78)
print(f"  HEADLINE      : scope-partition (AE governs FRW linear growth;")
print(f"                  AC's x_bar=5.85 is an admitted hypothesis the action does not produce)")
print(f"  geff_derived  : FALSE  (AQUAL G_eff is not the cosmological growth carrier; chi-matter is)")
print(f"  growth carrier: chi-matter (Omega_chi h^2=0.12, Q~1, GR-like)")
print(f"  sigma8(0)     : {sigma8_0:.3f}")
print(f"  S8            : {S8:.3f}   -> lensing-LOW side: {'YES' if S8<0.79 else 'NO'}")
print(f"  fsigma8 chi2/N: DFD {chi2_dfd/N:.2f}  vs  LCDM {chi2_lcdm/N:.2f}  (acceptable, NOT superior)")
print(f"  S8 vs Planck  : {(S8-0.834)/0.016:+.1f} sigma")
print(f"  S8 vs DES-Y3  : {(S8-0.776)/0.017:+.1f} sigma")
print(f"  S8 vs KiDS    : {(S8-0.759)/0.024:+.1f} sigma")
print("="*78)
