#!/usr/bin/env python3
"""LUMEN Second Edition -- Dark Sector & Cosmology chapter verifier.

Checks the geometric dark-energy density Lambda_0 = 1 - pi^2/32 against the measured Omega_Lambda,
the anchor-limited Hubble constant (H0 ~ sqrt(a0), so the prediction inherits a0's +/-22%), and the
gravity prediction log(M_Pl/m_e) = (3pi/20) mu1/(1 - mu1 alpha^2).

Cosmological rows are reported in sigma (the measurement error is comparable to the prediction);
the gravity ratio is reported as a relative residual. Method: exact arithmetic + error propagation.
Independently pinned by verify_P032.py, verify_P038.py, verify_P393.py, verify_P017.py.
Measured: Planck 2018 (Omega_L=0.6847+/-0.0073); SH0ES Riess+22 (73.04+/-1.04); Planck (67.36+/-0.54);
DESI 2024 (68.52+/-0.62); SPARC a0 = 1.20+/-0.26e-10 m/s^2; CODATA M_Pl, m_e."""
import sys, os
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import numpy as np
from numpy import pi, sqrt, log
from se_verify_common import Report

R = Report("Dark Sector & Cosmology",
    "the cokernel void Lambda_0, the anchor-limited H0, and the gravity ratio",
    "exact arithmetic + standard error propagation (H0 inherits the a0 anchor uncertainty)",
    sources=["Planck 2018 (Omega_L=0.6847+/-0.0073)", "SH0ES Riess+2022 (73.04+/-1.04)",
             "Planck (67.36+/-0.54)", "DESI 2024 (68.52+/-0.62)", "SPARC (a0=1.20+/-0.26e-10)"],
    pins=["verify_P032.py", "verify_P038.py", "verify_P393.py", "verify_P017.py"])

# dark energy = corner volume the inscribed 4-ball misses
filled = (pi**2/2)/16
Lam0 = 1 - filled
R.check("filled fraction vol(B^4)/vol(cube) = pi^2/32", filled, None, source="geometric")
R.check("Lambda_0 = 1 - pi^2/32", Lam0, 0.6847, source="Planck 2018", sigma=(Lam0-0.6847)/0.0073,
        note="the corner volume the ball cannot reach = the cokernel void")

# H0 ~ sqrt(a0): the prediction inherits the a0 anchor's +/-22%
H0 = 75.8
eH0 = H0 * (0.26/1.20)/2          # +/-22% on a0 -> +/-11% on H0
R.check("H0 prediction (anchor-limited)", H0, None, unit="km/s/Mpc", source="2 alpha^2/T_breath",
        note="+/-%.0f%% from the a0 anchor -> 75.8 +/- %.1f" % (100*(0.26/1.20)/2, eH0))
for nm, val, err in [("SH0ES", 73.04, 1.04), ("DESI", 68.52, 0.62), ("Planck", 67.36, 0.54)]:
    R.check("H0 vs %s" % nm, H0, val, unit="km/s/Mpc", source=nm,
            sigma=(H0-val)/sqrt(eH0**2+err**2), note="anchor-limited, not a 3-sigma tension")

# gravity: same density, first moment
trap = getattr(np, "trapezoid", np.trapz)
xg = np.linspace(0, 1, 1_000_001)
rho = 16*pi**3*xg**3 + 3*pi**2*xg**2 + 2*pi*xg
mu0 = float(trap(rho, xg)); mu1 = float(trap(xg*rho, xg)); alpha = 1/mu0
g_pred = (3*pi/20)*mu1/(1 - mu1*alpha**2)
M_Pl, m_e = 1.220890e19, 0.51099895e-3      # GeV
g_meas = log(M_Pl/m_e)
R.check("log(M_Pl/m_e) = (3pi/20) mu1/(1-mu1 alpha^2)", g_pred, g_meas, source="CODATA masses", tol=1e-3,
        note="gravity from the first moment; self-lensing factor 1/(1-mu1 alpha^2)")
R.emit()