#!/usr/bin/env python3
"""
verify_P127.py — Verification script for Addendum 127
Strong Coupling α_s — Geometric Matching at the Non-Perturbative Boundary

Verifies:
  1. μ_match = m_conf · exp(2π(6−√3)/23) ≈ 705.9 MeV
  2. LO   αs(MZ) ≈ 0.1261
  3. NLO  αs(MZ) ≈ 0.1183  (+0.4σ from PDG)
  4. N²LO αs(MZ) ≈ 0.1176  (−0.3σ from PDG)
  5. NLO and N²LO bracket the PDG central value 0.1179
  6. Matched series is monotone (LO > NLO > N²LO)

Note on exponent: the paper's eq.(T) writes ln(μ²_match/m²_conf) = 2π(6−√3)/23,
but the correct reading of the algebra — and confirmed by P132 — is that the
exponent in the boxed formula is 2π, i.e. ln(μ_match/m_conf) = 2π(6−√3)/23.
This script uses the boxed (correct) formula.

Copyright: Léon Fernando Vlegels. License: MIT.
Addendum 127 to the LumenOS TOE Corpus. Date: 2026-05-14.
"""

import math
import sys

# ─── TOE constants ────────────────────────────────────────────────────────────
ALPHA_INV     = 4 * math.pi**3 + math.pi**2 + math.pi   # ≈ 137.036304
ALPHA         = 1 / ALPHA_INV
BREATH_PERIOD = math.pi * ALPHA_INV                       # ≈ 432

# ─── Physical constants ───────────────────────────────────────────────────────
M_ELECTRON_MEV = 0.51099895000   # MeV  (PDG 2024)
MZ_MEV         = 91187.6         # MeV  (Z boson mass)
PDG_AS_MZ      = 0.1179          # PDG 2024 central value
PDG_AS_MZ_ERR  = 0.0009          # PDG 2024 uncertainty

# ─── TOE confinement scale (Addendum 106) ─────────────────────────────────────
# m_conf = π · α⁻¹ · m_e
M_CONF = math.pi * ALPHA_INV * M_ELECTRON_MEV

# ─── Matching scale (this addendum, eq. boxed) ────────────────────────────────
# G₂ root geometry: αs(μ_match) = sin(π/6) = 1/2
# Inverting one-loop running: ln(μ_match/m_conf) = 2π(6−√3)/23
EXPONENT = 2 * math.pi * (6 - math.sqrt(3)) / 23
MU_MATCH = M_CONF * math.exp(EXPONENT)
AS_MATCH = math.sin(math.pi / 6)   # = 0.5 exactly

# ─── MSbar β-function coefficients (nf=5, as used throughout the paper) ──────
Nc = 3
CF = 4 / 3
nf = 5
b0 = (11 * Nc - 2 * nf) / 3                                  # = 23/3
b1 = (34 * Nc**2 - 10 * Nc * nf - 6 * CF * nf) / 3          # = 116/3
b2 = 9769 / 54                                                # from J₃(𝕆) Casimir triple (P125)


# ─── LO running (analytical closed form) ──────────────────────────────────────
def run_lo_analytic(as_start: float, mu_start: float, mu_end: float) -> float:
    """
    LO: 1/αs(μ2) = 1/αs(μ1) + (b0/2π)·ln(μ2/μ1)
    ↔  αs(μ2) = αs(μ1) / (1 + αs(μ1)·(b0/2π)·ln(μ2/μ1))
    """
    return as_start / (1 + as_start * (b0 / (2 * math.pi)) * math.log(mu_end / mu_start))


# ─── NLO / N²LO RGE right-hand side ──────────────────────────────────────────
def rge_rhs(a: float, order: int) -> float:
    """
    dαs / d(lnμ) = −(b0/2π)·αs² · [1 + (b1/4π·b0)·αs + (b2/(4π)²·b0)·αs²]
    order=1: LO bracket (just 1)
    order=2: NLO bracket (adds b1 term)
    order=3: N²LO bracket (adds b2 term)
    """
    bracket = 1.0
    if order >= 2:
        bracket += (b1 / (4 * math.pi * b0)) * a
    if order >= 3:
        bracket += (b2 / ((4 * math.pi)**2 * b0)) * a**2
    return -(b0 / (2 * math.pi)) * a**2 * bracket


def run_rk4(as_start: float, mu_start: float, mu_end: float, order: int,
            n_steps: int = 200_000) -> float:
    """
    4th-order Runge–Kutta integration of dαs/d(lnμ).
    Uses lnμ as the independent variable for numerical stability.
    """
    t_start = math.log(mu_start)
    t_end   = math.log(mu_end)
    h       = (t_end - t_start) / n_steps
    a       = as_start
    for _ in range(n_steps):
        k1 = rge_rhs(a,            order)
        k2 = rge_rhs(a + h*k1/2,  order)
        k3 = rge_rhs(a + h*k2/2,  order)
        k4 = rge_rhs(a + h*k3,    order)
        a += h * (k1 + 2*k2 + 2*k3 + k4) / 6
    return a


# ─── Run ──────────────────────────────────────────────────────────────────────
as_lo   = run_lo_analytic(AS_MATCH, MU_MATCH, MZ_MEV)
as_nlo  = run_rk4(AS_MATCH, MU_MATCH, MZ_MEV, order=2)
as_n2lo = run_rk4(AS_MATCH, MU_MATCH, MZ_MEV, order=3)

# ─── Paper claims ─────────────────────────────────────────────────────────────
CLAIM_MU_MATCH = 705.9    # MeV
CLAIM_LO       = 0.1261
CLAIM_NLO      = 0.1183
CLAIM_N2LO     = 0.1176

# ─── Tolerances ───────────────────────────────────────────────────────────────
TOL_MU   = 1.0    # MeV  (paper rounds to 705.9)
TOL_LO   = 5e-4   # 4th decimal place
TOL_NLO  = 2e-4
TOL_N2LO = 2e-4

# ─── Helpers ──────────────────────────────────────────────────────────────────
PASS = FAIL = 0

def check(n, desc, cond):
    global PASS, FAIL
    ok = bool(cond); PASS += ok; FAIL += (not ok)
    print(f"  [{'PASS' if ok else 'FAIL'}] {n:>2}. {desc}")


# ─── Report ───────────────────────────────────────────────────────────────────
print("Addendum P127 — Strong Coupling Boundary")

print()
print("TOE constants")
print(f"  ALPHA_INV     = {ALPHA_INV:.8f}   (exact: 4π³+π²+π)")
print(f"  BREATH_PERIOD = {BREATH_PERIOD:.6f}   (π · α⁻¹ ≈ 432)")

print()
print("Matching scale derivation")
print(f"  m_conf        = {M_CONF:.4f} MeV   (π · α⁻¹ · m_e, Addendum 106)")
print(f"  exponent T    = 2π(6−√3)/23 = {EXPONENT:.6f}")
print(f"  μ_match       = {MU_MATCH:.4f} MeV   (claim: {CLAIM_MU_MATCH} MeV)")
print(f"  αs(μ_match)   = {AS_MATCH:.4f}        (sin(π/6) = 1/2, exact)")

print()
print("β coefficients (nf=5)")
print(f"  b0 = {b0:.6f}   (= 23/3  ✓)")
print(f"  b1 = {b1:.6f}   (= 116/3 ✓, Addendum 121)")
print(f"  b2 = {b2:.6f}   (= 9769/54, Addendum 125)")

print()
delta_t = 2 * math.log(MZ_MEV / MU_MATCH)
print(f"  Δt = ln(MZ²/μ_match²) = {delta_t:.4f}   (paper: 9.722)")

print()
print("PDG consistency")
print(f"  PDG αs(MZ) = {PDG_AS_MZ} ± {PDG_AS_MZ_ERR}")
for label, computed in [("LO", as_lo), ("NLO", as_nlo), ("N²LO", as_n2lo)]:
    sigma = (computed - PDG_AS_MZ) / PDG_AS_MZ_ERR
    print(f"  {label:<8}  αs = {computed:.4f}  →  {sigma:+.1f}σ from PDG")

print()
print("S1  Running to MZ")
for i, (label, computed, claimed, tol) in enumerate([
    ("LO",   as_lo,   CLAIM_LO,   TOL_LO),
    ("NLO",  as_nlo,  CLAIM_NLO,  TOL_NLO),
    ("N²LO", as_n2lo, CLAIM_N2LO, TOL_N2LO),
], start=1):
    diff = computed - claimed
    check(i, f"{label} αs(MZ) = {computed:.4f} vs claimed {claimed:.4f} "
             f"(diff {diff:+.5f})", abs(computed - claimed) <= tol)

print("S2  Structural checks")
check(4, f"bracketing (NLO > PDG > N²LO): {as_nlo:.4f} > {PDG_AS_MZ} > {as_n2lo:.4f}",
      as_nlo > PDG_AS_MZ > as_n2lo)
check(5, f"monotone convergence (LO > NLO > N²LO): {as_lo:.4f} > {as_nlo:.4f} > {as_n2lo:.4f}",
      as_lo > as_nlo > as_n2lo)
check(6, f"matching scale {MU_MATCH:.2f} MeV vs claimed {CLAIM_MU_MATCH} MeV "
         f"(diff {MU_MATCH - CLAIM_MU_MATCH:+.2f} MeV)",
      abs(MU_MATCH - CLAIM_MU_MATCH) <= TOL_MU)

print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
sys.exit(0 if FAIL == 0 else 1)