"""
verify_P129.py
Standalone verification of Addendum P129: Muon Lifetime from the TOE Fermi Sector.

Paper claims:
  - LO  (TOE GF + PDG m_mu):  tau = 2.1873e-6 s,  residual -0.44%
  - NLO (+ leading QED corr): tau = 2.1945e-6 s,  residual -0.11%  ← primary claim
  - Full-TOE LO  (both TOE):  tau = 2.1221e-6 s,  residual -3.41%
  - Full-TOE NLO:             tau = 2.1290e-6 s,  residual -3.09%

Copyright: Léon Fernando Vlegels. License: MIT.
"""

import math
import sys

# ──────────────────────────────────────────────
# 1. Constants
# ──────────────────────────────────────────────

# TOE fine-structure constant (P01/P36)
ALPHA_INV = 4 * math.pi**3 + math.pi**2 + math.pi   # ≈ 137.036
ALPHA     = 1.0 / ALPHA_INV

# Fermi coupling (GeV^-2) — TOE value from P123
GF_TOE = 1.16638e-5   # GeV^-2
GF_PDG = 1.1663788e-5 # GeV^-2

# Muon mass (GeV)
M_MU_PDG = 105.658e-3  # GeV
M_MU_TOE = 106.300e-3  # GeV  (A100 spectral ladder)

# Electron mass (GeV)
M_E = 0.51100e-3  # GeV

# Reduced Planck constant
HBAR = 6.582119e-25  # GeV·s

# PDG measured muon lifetime
TAU_PDG = 2.1969811e-6  # s

# ──────────────────────────────────────────────
# 2. LO formula   tau = 192π³ ħ / (GF² mμ⁵)
# ──────────────────────────────────────────────

def tau_LO(GF, m_mu):
    """Leading-order muon lifetime (seconds)."""
    return (192.0 * math.pi**3 * HBAR) / (GF**2 * m_mu**5)


# ──────────────────────────────────────────────
# 3. Leading QED radiative correction
# ──────────────────────────────────────────────

def kinoshita_sirlin(m_e, m_mu):
    """
    Kinoshita-Sirlin phase-space factor f(x), x = me/mmu.
    f(x) = 1 - 8x² + 8x⁶ - x⁸ - 24x⁴ ln(x)
    Corrects for finite electron mass in 3-body phase space.
    A factor <1 reduces decay rate, i.e. increases lifetime.
    """
    x = m_e / m_mu
    return 1.0 - 8*x**2 + 8*x**6 - x**8 - 24*x**4 * math.log(x)


def vertex_correction(alpha):
    """
    Leading O(alpha) vertex correction.
    delta_vert = (25/4 - pi²/2) * alpha/pi
    Reduces decay rate (positive delta_vert → rate multiplied by (1 - delta_vert)).
    """
    coeff = 25.0/4.0 - math.pi**2/2.0   # ≈ 1.3152
    return coeff * alpha / math.pi


def F_QED(m_e, m_mu, alpha):
    """
    Combined leading QED factor applied to the decay rate:
      Gamma_NLO = Gamma_LO * F_QED
    So tau_NLO = tau_LO / F_QED.
    """
    f  = kinoshita_sirlin(m_e, m_mu)
    dv = vertex_correction(alpha)
    return f * (1.0 - dv)


# ──────────────────────────────────────────────
# 4. Compute all branches
# ──────────────────────────────────────────────

def pct(computed, reference):
    return 100.0 * (computed - reference) / reference


def run():
    failures = []

    # ── Main branch: TOE GF, PDG mmu ──────────
    t_lo_main  = tau_LO(GF_TOE, M_MU_PDG)
    fqed_main  = F_QED(M_E, M_MU_PDG, ALPHA)
    t_nlo_main = t_lo_main / fqed_main

    res_lo_main  = pct(t_lo_main,  TAU_PDG)
    res_nlo_main = pct(t_nlo_main, TAU_PDG)

    # ── Full-TOE branch: TOE GF, TOE mmu ──────
    t_lo_toe  = tau_LO(GF_TOE, M_MU_TOE)
    fqed_toe  = F_QED(M_E, M_MU_TOE, ALPHA)
    t_nlo_toe = t_lo_toe / fqed_toe

    res_lo_toe  = pct(t_lo_toe,  TAU_PDG)
    res_nlo_toe = pct(t_nlo_toe, TAU_PDG)

    # ── Intermediate diagnostics ───────────────
    x     = M_E / M_MU_PDG
    f_val = kinoshita_sirlin(M_E, M_MU_PDG)
    dv    = vertex_correction(ALPHA)

    # ──────────────────────────────────────────
    # 5. Print results
    # ──────────────────────────────────────────
    print("=" * 60)
    print("P129 Muon Lifetime Verification")
    print("=" * 60)
    print()
    print(f"  ALPHA_INV               = {ALPHA_INV:.6f}  (paper: ~137.036)")
    print(f"  GF_TOE                  = {GF_TOE:.5e} GeV^-2")
    print(f"  m_mu PDG                = {M_MU_PDG*1e3:.3f} MeV")
    print(f"  m_mu TOE                = {M_MU_TOE*1e3:.3f} MeV")
    print(f"  hbar                    = {HBAR:.6e} GeV·s")
    print()
    print("── Kinoshita-Sirlin factor ─────────────────")
    print(f"  x = me/mmu              = {x:.4e}  (paper: 4.836e-3)")
    print(f"  f(x)                    = {f_val:.6f}  (paper: 0.999813)")
    print()
    print("── Vertex correction ───────────────────────")
    print(f"  25/4 - pi²/2            = {25/4 - math.pi**2/2:.4f}  (paper: 1.315)")
    print(f"  delta_vert              = {dv:.6e}  (paper: 3.056e-3)")
    print()
    print(f"  F_QED (main branch)     = {fqed_main:.6f}  (paper: 0.996759)")
    print()
    print("── Main branch (TOE GF + PDG m_mu) ─────────")
    print(f"  tau_LO                  = {t_lo_main:.4e} s  (paper: 2.1873e-6 s)")
    print(f"  residual LO             = {res_lo_main:+.3f}%         (paper: -0.44%)")
    print(f"  tau_NLO                 = {t_nlo_main:.4e} s  (paper: 2.1945e-6 s)")
    print(f"  residual NLO            = {res_nlo_main:+.3f}%         (paper: -0.11%)")
    print()
    print("── Full-TOE branch (TOE GF + TOE m_mu) ─────")
    print(f"  tau_LO                  = {t_lo_toe:.4e} s  (paper: 2.1221e-6 s)")
    print(f"  residual LO             = {res_lo_toe:+.3f}%         (paper: -3.41%)")
    print(f"  tau_NLO                 = {t_nlo_toe:.4e} s  (paper: 2.1290e-6 s)")
    print(f"  residual NLO            = {res_nlo_toe:+.3f}%         (paper: -3.09%)")
    print()

    # ──────────────────────────────────────────
    # 6. Assertions  (tolerance ±0.02 percentage-points)
    # ──────────────────────────────────────────
    TOL = 0.02  # percentage-point tolerance

    checks = [
        ("f(x) Kinoshita-Sirlin",        f_val,        0.999813, 1e-5),
        ("delta_vert",                    dv,           3.056e-3, 1e-5),
        ("F_QED main branch",             fqed_main,    0.996759, 1e-5),
        ("tau_LO main [1e-6 s]",          t_lo_main*1e6,  2.1873, 5e-4),
        ("residual LO main [%]",          res_lo_main,   -0.44,   TOL),
        ("tau_NLO main [1e-6 s]",         t_nlo_main*1e6, 2.1945, 5e-4),
        ("residual NLO main [%]",         res_nlo_main,  -0.11,   TOL),   # PRIMARY CLAIM
        ("tau_LO full-TOE [1e-6 s]",      t_lo_toe*1e6,   2.1221, 5e-4),
        ("residual LO full-TOE [%]",      res_lo_toe,    -3.41,   TOL),
        ("tau_NLO full-TOE [1e-6 s]",     t_nlo_toe*1e6,  2.1290, 5e-4),
        ("residual NLO full-TOE [%]",     res_nlo_toe,   -3.09,   TOL),
    ]

    print("── Assertion results ────────────────────────")
    PASS = FAIL = 0
    for n, (name, got, expected, tol) in enumerate(checks, 1):
        ok = abs(got - expected) <= tol
        PASS += ok
        FAIL += (not ok)
        status = "PASS" if ok else "FAIL"
        print(f"  [{status}] {n:>2}. {name}")
        print(f"         got={got:.6g}  expected={expected:.6g}  tol=±{tol:.1g}")
        if not ok:
            failures.append((name, got, expected, tol))

    print()
    if not failures:
        print("Primary claim: tau_NLO (main branch) residual vs PDG = "
              f"{res_nlo_main:+.3f}%  (paper: -0.11%)")
    print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
    if failures:
        sys.exit(1)


if __name__ == "__main__":
    run()
