#!/usr/bin/env python3
"""
verify_P135.py — Computational verification of Addendum P135
"W-Boson Mass via the On-Shell Fermi Consistency Equation"

Re-derives every numerical claim in P135 from TOE first principles,
using only the inputs stated in the addendum: α (exact TOE), GF (P123),
MZ (PDG), Δα_lep (P134), Δα_had,SM (PDG 5-quark).

Each claim is compared against the value stated in the paper.
Tolerances are ≤ 0.2% for physical predictions (generous for radiative
corrections) and tight absolutes for intermediate quantities.

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

import math
import sys

# ══════════════════════════════════════════════════════════════════════════════
# TOE / stated inputs (as given in P135 §2 "Computation")
# ══════════════════════════════════════════════════════════════════════════════
ALPHA_INV  = 4*math.pi**3 + math.pi**2 + math.pi  # exact TOE expression ≈ 137.0363
ALPHA      = 1 / ALPHA_INV                          # ≈ 0.0072973363
GF         = 1.16638e-5   # GeV^-2  (P123; matches PDG to 6 s.f.)
MZ         = 91.1876      # GeV     (PDG Z-pole mass)
PHI        = (1 + math.sqrt(5)) / 2                 # golden ratio
SIN2_W     = 3 / (8 * PHI)                          # exact Weinberg angle (P112)
COS2_W     = 1 - SIN2_W
DALEP      = 0.031422     # leptonic running (P134)
DAHAD_SM   = 0.02750      # SM LO 5-quark hadronic running (PDG)
MW_PDG     = 80.377       # GeV  (PDG W mass used throughout P135)
MW_P122    = 80.313       # GeV  (P122 geometric route result; taken as given)
DRHO_P122  = 0.009718     # Veltman ρ-correction from P122

# ══════════════════════════════════════════════════════════════════════════════
# Core formula: on-shell Fermi consistency equation
#   MW²(1 − MW²/MZ²) = π α(MZ) / (√2 GF)
#
# Physical (upper) root of the equivalent quadratic in MW²:
#   MW² = MZ²/2 · [1 + √(1 − 4πα(MZ)/(√2 GF MZ²))]
# ══════════════════════════════════════════════════════════════════════════════
def fermi_mw(dalpha):
    """Return (MW, alpha_mz, discriminant) for given total running Δα."""
    alpha_mz     = ALPHA / (1 - dalpha)
    discriminant = 1 - 4 * math.pi * alpha_mz / (math.sqrt(2) * GF * MZ**2)
    mw           = math.sqrt(MZ**2 / 2 * (1 + math.sqrt(discriminant)))
    return mw, alpha_mz, discriminant


# ══════════════════════════════════════════════════════════════════════════════
# Result collector
# ══════════════════════════════════════════════════════════════════════════════
_results = []

def check(label, computed, claimed, *, tol_pct=None, tol_abs=None):
    """
    Compare computed vs claimed.  Supply exactly one of tol_pct (relative %)
    or tol_abs (absolute).  Accumulate pass/fail.
    """
    if tol_abs is not None:
        err  = abs(computed - claimed)
        ok   = err <= tol_abs
        info = f"|Δ| = {err:.4g}   tol = {tol_abs}"
    else:
        if claimed == 0:
            err_pct = 0.0 if computed == 0 else 1e9
        else:
            err_pct = abs(computed - claimed) / abs(claimed) * 100
        ok   = err_pct <= tol_pct
        info = f"|Δ%| = {err_pct:.5f}%   tol = {tol_pct}%"

    _results.append((label, ok))
    n = len(_results)
    print(f"  [{'PASS' if ok else 'FAIL'}] {n:>2}. {label}")
    print(f"       computed = {computed:.8g}   claimed = {claimed:.8g}   {info}")


# ══════════════════════════════════════════════════════════════════════════════
print("=" * 72)
print("Addendum P135 — Independent Computational Verification")
print("W-Boson Mass via the On-Shell Fermi Consistency Equation")
print("=" * 72)


# ──────────────────────────────────────────────────────────────────────────────
# §1  TOE constants
# ──────────────────────────────────────────────────────────────────────────────
print("\n[1]  TOE constants ───────────────────────────────────────────────────")

check(
    "ALPHA_INV = 4π³ + π² + π  (≈ 137.0363)",
    ALPHA_INV, 137.0363,
    tol_pct=0.001,
)
check(
    "ALPHA = 1/(4π³+π²+π)  (≈ 0.0072973363)",
    ALPHA, 0.0072973363,
    tol_pct=0.001,
)
# SIN2_W = 3/(8φ); paper uses it numerically in §3 via cos²θW/sin²θW = 3.3148
cos_sin_ratio_claimed = 3.3148
cos_sin_ratio_derived = COS2_W / SIN2_W
check(
    "cos²θW / sin²θW = (1 − 3/8φ) / (3/8φ)  (≈ 3.3148)",
    cos_sin_ratio_derived, cos_sin_ratio_claimed,
    tol_pct=0.05,
)


# ──────────────────────────────────────────────────────────────────────────────
# §2  Leptonic-only prediction  (Δα = Δα_lep = 0.031422)
# ──────────────────────────────────────────────────────────────────────────────
print("\n[2]  Leptonic-only prediction  (Δα = Δα_lep = 0.031422) ─────────────")

MW_lep, alpha_lep, delta_lep = fermi_mw(DALEP)

check(
    "α_lep(MZ) = α / (1 − Δα_lep)  (≈ 0.0075341)",
    alpha_lep, 0.0075341,
    tol_pct=0.01,
)
check(
    "1/α_lep  (≈ 132.730)",
    1 / alpha_lep, 132.730,
    tol_pct=0.01,
)
check(
    "δ_lep = 1 − 4πα_lep/(√2·GF·MZ²)  (≈ 0.30974)",
    delta_lep, 0.30974,
    tol_pct=0.05,
)
check(
    "MW_lep = MZ·√[½(1+√δ_lep)]  (≈ 80.4454 GeV)",
    MW_lep, 80.4454,
    tol_pct=0.02,
)
res_lep_pct = 100 * (MW_lep - MW_PDG) / MW_PDG
check(
    "MW_lep residual vs PDG  (claimed +0.085%)",
    res_lep_pct, 0.085,
    tol_abs=0.030,           # absolute on the %-value
)


# ──────────────────────────────────────────────────────────────────────────────
# §3  Baseline: Δα = 0  (bare coupling, no running)
# ──────────────────────────────────────────────────────────────────────────────
print("\n[3]  Baseline: Δα = 0 (no running) ─────────────────────────────────")

MW_0, _alpha_0, _delta_0 = fermi_mw(0.0)

check(
    "MW_0 (bare, no running)  (≈ 80.939 GeV)",
    MW_0, 80.939,
    tol_pct=0.02,
)
res_0_pct = 100 * (MW_0 - MW_PDG) / MW_PDG
check(
    "MW_0 residual vs PDG  (claimed +0.699%)",
    res_0_pct, 0.699,
    tol_abs=0.050,
)
lep_shift = MW_0 - MW_lep
check(
    "Leptonic running shift  MW_0 − MW_lep  (claimed 0.494 GeV)",
    lep_shift, 0.494,
    tol_abs=0.005,
)


# ──────────────────────────────────────────────────────────────────────────────
# §4  Full SM hadronic prediction  (Δα = Δα_lep + Δα_had,SM)
# ──────────────────────────────────────────────────────────────────────────────
print("\n[4]  Full SM hadronic  (Δα = Δα_lep + 0.02750) ──────────────────────")

DA_full = DALEP + DAHAD_SM   # exact: 0.031422 + 0.027500 = 0.058922
# Paper writes "0.05892" — 5-significant-figure rounding of 0.058922
check(
    "Δα_total = Δα_lep + Δα_had,SM  (paper: 0.05892)",
    DA_full, 0.05892,
    tol_abs=1e-4,
)

MW_full, alpha_full, delta_full = fermi_mw(DA_full)

check(
    "α_full(MZ) = α/(1−Δα_total)  (≈ 0.0077542)",
    alpha_full, 0.0077542,
    tol_pct=0.01,
)
check(
    "1/α_full  (≈ 128.962)",
    1 / alpha_full, 128.962,
    tol_pct=0.01,
)
check(
    "δ_full discriminant  (≈ 0.28957)",
    delta_full, 0.28957,
    tol_pct=0.05,
)
check(
    "MW_full  (≈ 79.9678 GeV)",
    MW_full, 79.9678,
    tol_pct=0.02,
)
res_full_pct = 100 * (MW_full - MW_PDG) / MW_PDG
check(
    "MW_full residual vs PDG  (claimed −0.509%)",
    res_full_pct, -0.509,
    tol_abs=0.030,
)
had_shift = MW_lep - MW_full
check(
    "Hadronic running shift  MW_lep − MW_full  (claimed 0.478 GeV)",
    had_shift, 0.478,
    tol_abs=0.005,
)


# ──────────────────────────────────────────────────────────────────────────────
# §5  Hadronic correction needed to reproduce PDG MW exactly
# ──────────────────────────────────────────────────────────────────────────────
print("\n[5]  Hadronic gap / numerical target ─────────────────────────────────")

# RHS of the Fermi equation evaluated at MW = MW_PDG
fermi_rhs_pdg = MW_PDG**2 * (1 - MW_PDG**2 / MZ**2)
check(
    "MW²(1−MW²/MZ²)|_{PDG}  (claimed 1441.018 GeV²)",
    fermi_rhs_pdg, 1441.018,
    tol_pct=0.01,
)

# Running coupling needed to match that RHS
alpha_needed = fermi_rhs_pdg * math.sqrt(2) * GF / math.pi
check(
    "α_needed(MZ) = RHS·√2·GF/π  (≈ 0.0075661)",
    alpha_needed, 0.0075661,
    tol_pct=0.01,
)
check(
    "1/α_needed  (≈ 132.168)",
    1 / alpha_needed, 132.168,
    tol_pct=0.01,
)

# Total running needed  (from α_needed = α/(1−Δα_needed))
DA_needed = 1 - ALPHA / alpha_needed
check(
    "Δα_needed = 1 − α/α_needed  (≈ 0.035528)",
    DA_needed, 0.035528,
    tol_pct=0.5,
)

# Hadronic target: subtract the already-known leptonic piece
DAHAD_target = DA_needed - DALEP
check(
    "Δα_had,target = Δα_needed − Δα_lep  (claimed 0.00411)",
    DAHAD_target, 0.00411,
    tol_pct=2.0,
)

# How far below SM value is the target?
factor = DAHAD_SM / DAHAD_target
check(
    "Factor  Δα_had,SM / Δα_had,target  (claimed 6.7)",
    factor, 6.7,
    tol_pct=5.0,
)


# ──────────────────────────────────────────────────────────────────────────────
# §6  P122 / P135 bracket and average
# ──────────────────────────────────────────────────────────────────────────────
print("\n[6]  Bracket: P122 (−0.080%) and P135 leptonic (+0.085%) ───────────")

MW_avg = (MW_lep + MW_P122) / 2
check(
    "(MW_lep + MW_P122)/2  (claimed 80.379 GeV)",
    MW_avg, 80.379,
    tol_abs=0.002,
)
bracket_gap = MW_avg - MW_PDG
check(
    "Bracket gap  (MW_avg − MW_PDG)  (claimed 0.002 GeV)",
    bracket_gap, 0.002,
    tol_abs=0.001,
)
res_avg_pct = 100 * (MW_avg - MW_PDG) / MW_PDG
check(
    "Bracket residual vs PDG  (claimed +0.003%)",
    res_avg_pct, 0.003,
    tol_abs=0.002,
)

# Verify P122 residual using its stated MW value
res_P122_pct = 100 * (MW_P122 - MW_PDG) / MW_PDG
check(
    "P122 residual vs PDG  (claimed −0.080%)",
    res_P122_pct, -0.080,
    tol_abs=0.005,
)


# ──────────────────────────────────────────────────────────────────────────────
# §7  Near-cancellation in the leading-order ΔR formula
#     ΔR⁽⁰⁾ = Δα_lep − Δρ·(cos²θW/sin²θW)
# ──────────────────────────────────────────────────────────────────────────────
print("\n[7]  Near-cancellation in ΔR⁽⁰⁾ formula ─────────────────────────────")

check(
    "cos²θW / sin²θW  (≈ 3.3148)  [repeated from §1 for context]",
    cos_sin_ratio_derived, 3.3148,
    tol_pct=0.05,
)

drho_term = DRHO_P122 * cos_sin_ratio_derived
check(
    "Δρ · (cos²θW/sin²θW) = 0.009718 × 3.3148  (≈ 0.032213)",
    drho_term, 0.032213,
    tol_pct=0.05,
)

ratio_lep_drho = DALEP / drho_term
check(
    "Δα_lep / [Δρ·(cos²/sin²)]  (claimed 0.9754)",
    ratio_lep_drho, 0.9754,
    tol_pct=0.10,
)

Dr0 = DALEP - drho_term
check(
    "ΔR⁽⁰⁾ = Δα_lep − Δρ·cos²θW/sin²θW  (claimed −0.000791)",
    Dr0, -0.000791,
    tol_abs=5e-5,
)

# Sanity: the near-cancellation implies |ΔR⁽⁰⁾| << Δα_lep (within ~2.5%)
cancellation_pct = abs(Dr0) / DALEP * 100
print(f"\n       [info]  Near-cancellation: |ΔR⁽⁰⁾| / Δα_lep = {cancellation_pct:.2f}%  "
      f"(paper: within 2.5% → ratio ≈ 0.9754)")


# ──────────────────────────────────────────────────────────────────────────────
# §8  P112 tree-level W mass   MW,tree = MZ·√(1 − sin²θW)
# ──────────────────────────────────────────────────────────────────────────────
print("\n[8]  P112 tree-level mass  MW,tree = MZ·√(1 − 3/(8φ)) ───────────────")

MW_tree = MZ * math.sqrt(1 - SIN2_W)
check(
    "MW,tree (≈ 79.925 GeV)",
    MW_tree, 79.925,
    tol_pct=0.02,
)


# ══════════════════════════════════════════════════════════════════════════════
# Summary
# ══════════════════════════════════════════════════════════════════════════════
n_pass  = sum(1 for _, ok in _results if ok)
n_fail  = sum(1 for _, ok in _results if not ok)

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