#!/usr/bin/env python3
"""
verify_P138.py — Verification script for Addendum P138
Pion FSI and rho-omega Interference Corrections to Delta_alpha_had

Verified claims:
  1. FSI Omnes correction:        +2.182e-4
  2. rho-omega cross-term:        -2.316e-4
  3. Net correction:              -1.34e-5
  4. P138 total Delta_alpha_had:  3.718e-3
  5. Residual gap:                -9.55%
  6. W-boson mass:                80.383 GeV

Author: verification script for León Fernando Vlegels, P138
Copyright: León Fernando Vlegels. License: MIT.
"""

import math
import sys

# ---------------------------------------------------------------------------
# TOE constants
# ---------------------------------------------------------------------------
ALPHA_INV = 4 * math.pi**3 + math.pi**2 + math.pi
ALPHA     = 1.0 / ALPHA_INV
PHI       = (1.0 + math.sqrt(5.0)) / 2.0

# ---------------------------------------------------------------------------
# Physical constants (PDG 2022)
# ---------------------------------------------------------------------------
M_Z    = 91.1876          # GeV
M_PI   = 0.13957          # GeV, charged pion mass
M_RHO  = 0.77526          # GeV, rho(770) mass
G_RHO  = 0.1478           # GeV, rho total width
M_OME  = 0.78265          # GeV, omega mass
G_OME  = 0.00849          # GeV, omega total width
G_F    = 1.1663787e-5     # GeV^-2, Fermi constant
ALPHA_S_MZ = 0.1186       # strong coupling at M_Z  (P106)

# Derived
S_THRESH = 4.0 * M_PI**2  # 4m_pi^2  (GeV^2)
M_RHO2   = M_RHO**2
M_OME2   = M_OME**2
M_Z2     = M_Z**2

# ---------------------------------------------------------------------------
# P118 rho-omega mixing amplitude
# ---------------------------------------------------------------------------
EPS = -ALPHA_S_MZ / 3.0   # = -0.03953...

# ---------------------------------------------------------------------------
# P137 baseline and targets (from P138 text)
# ---------------------------------------------------------------------------
DA_HAD_P137  = 3.731e-3   # GS baseline
DA_HAD_FERMI = 4.110e-3   # Fermi target (P135)
DA_LEP       = 0.031422   # leptonic running (P134)
PDG_MW       = 80.377     # GeV

# ---------------------------------------------------------------------------
# Helper functions
# ---------------------------------------------------------------------------

def beta(s):
    """Pion velocity beta(s) = sqrt(1 - 4m_pi^2/s)."""
    arg = 1.0 - S_THRESH / s
    if arg <= 0.0:
        return 0.0
    return math.sqrt(arg)


def bw_ff_sq(s):
    """|F_pi|^2_BW = M_rho^4 / [(s-M_rho^2)^2 + M_rho^2 Gamma_rho^2]."""
    denom = (s - M_RHO2)**2 + M_RHO2 * G_RHO**2
    return M_RHO2**2 / denom


def omnes_factor(s):
    """
    Omnes ratio |F_pi|^2_Om / |F_pi|^2_BW = s^(-s/M_rho^2)
    using linear phase-shift approximation delta_1^1(s) ~ (s/M_rho^2)*(pi/2)
    and Lambda = 1 GeV cutoff.
    = exp[(s/M_rho^2) * ln(1/s)]  with s in GeV^2.
    """
    exp_arg = (s / M_RHO2) * math.log(1.0 / s)
    return math.exp(exp_arg)


def re_b_omega(s):
    """Re[B_omega(s)] = M_omega^2 (M_omega^2 - s) / [(M_omega^2-s)^2 + (M_omega Gamma_omega)^2]."""
    num   = M_OME2 * (M_OME2 - s)
    denom = (M_OME2 - s)**2 + (M_OME * G_OME)**2
    return num / denom


def dispersive_kernel(s):
    """Common dispersive weight: M_Z^2 / (M_Z^2 - s)."""
    return M_Z2 / (M_Z2 - s)


# ---------------------------------------------------------------------------
# Gaussian-Legendre quadrature (pure Python, no scipy dependency)
# ---------------------------------------------------------------------------

def gauss_legendre_nodes_weights(n):
    """
    Compute n-point Gauss-Legendre nodes and weights on [-1, 1].
    Uses Newton iteration on Legendre polynomial recurrence.
    """
    nodes   = [0.0] * n
    weights = [0.0] * n
    for i in range((n + 1) // 2):
        # Initial guess
        xi = math.cos(math.pi * (i + 0.75) / (n + 0.5))
        for _ in range(100):
            p0, p1 = 1.0, xi
            for k in range(2, n + 1):
                p0, p1 = p1, ((2*k - 1)*xi*p1 - (k-1)*p0) / k
            dp = n * (p0 - xi * p1) / (1.0 - xi**2)
            dx = p1 / dp
            xi -= dx
            if abs(dx) < 1e-15:
                break
        w = 2.0 / ((1.0 - xi**2) * dp**2)
        nodes[i]         = -xi
        nodes[n - 1 - i] =  xi
        weights[i]         = w
        weights[n - 1 - i] = w
    return nodes, weights


def integrate_gl(f, a, b, n):
    """Integrate f on [a,b] using n-point Gauss-Legendre."""
    nodes, weights = gauss_legendre_nodes_weights(n)
    mid  = 0.5 * (a + b)
    half = 0.5 * (b - a)
    return half * sum(w * f(mid + half * x) for x, w in zip(nodes, weights))


def integrate_subintervals(f, breakpoints, n_per_sub):
    """Integrate f over subintervals defined by a list of breakpoints."""
    total = 0.0
    for k in range(len(breakpoints) - 1):
        total += integrate_gl(f, breakpoints[k], breakpoints[k+1], n_per_sub)
    return total


# ---------------------------------------------------------------------------
# Claim 1: FSI Omnes correction
# ---------------------------------------------------------------------------

def da_fsi_integrand(s):
    b   = beta(s)
    ff2 = bw_ff_sq(s)
    om  = omnes_factor(s)
    dR  = 0.25 * b**3 * ff2 * (om - 1.0)
    return dR / s * dispersive_kernel(s)


# Integration range: 4m_pi^2 to (0.6 GeV)^2 = 0.36 GeV^2
S_HI_FSI = 0.36   # GeV^2

# Five sub-intervals matching the paper's Table 1 (sqrt(s) boundaries in GeV)
SQRT_S_BREAKS_FSI = [math.sqrt(S_THRESH), 0.346, 0.400, 0.458, 0.520, 0.600]
S_BREAKS_FSI = [x**2 for x in SQRT_S_BREAKS_FSI]

DA_FSI = (ALPHA / (3.0 * math.pi)) * integrate_subintervals(
    da_fsi_integrand, S_BREAKS_FSI, n_per_sub=400
)

# ---------------------------------------------------------------------------
# Claim 2: rho-omega cross-term
# ---------------------------------------------------------------------------

def da_rw_integrand(s):
    b   = beta(s)
    ff2 = bw_ff_sq(s)
    reBw = re_b_omega(s)
    dR  = 0.25 * b**3 * ff2 * 2.0 * EPS * reBw
    return dR / s * dispersive_kernel(s)


# Integration range: 4m_pi^2 to 1 GeV^2, with fine grid around omega pole
S_HI_RW = 1.0   # GeV^2
DELTA_OMEGA = 5.0 * M_OME * G_OME   # half-window around omega pole

S_BREAKS_RW = [
    S_THRESH,
    max(S_THRESH + 1e-6, M_OME2 - DELTA_OMEGA),
    M_OME2 - 0.5 * M_OME * G_OME,
    M_OME2,
    M_OME2 + 0.5 * M_OME * G_OME,
    M_OME2 + DELTA_OMEGA,
    S_HI_RW,
]
# Remove any breakpoints below threshold or above upper limit, keep sorted unique
S_BREAKS_RW = sorted(set(
    max(S_THRESH + 1e-9, min(S_HI_RW - 1e-9, bp)) for bp in S_BREAKS_RW
    if S_THRESH < bp < S_HI_RW
))
S_BREAKS_RW = [S_THRESH] + S_BREAKS_RW + [S_HI_RW]

DA_RW = (ALPHA / (3.0 * math.pi)) * integrate_subintervals(
    da_rw_integrand, S_BREAKS_RW, n_per_sub=600
)

# ---------------------------------------------------------------------------
# Derived claims
# ---------------------------------------------------------------------------

DA_NET    = DA_FSI + DA_RW
DA_P138   = DA_HAD_P137 + DA_NET
GAP_ABS   = DA_P138 - DA_HAD_FERMI
GAP_PCT   = GAP_ABS / DA_HAD_FERMI * 100.0

# W-boson mass from Fermi on-shell constraint
DA_TOTAL_P138 = DA_LEP + DA_P138
ALPHA_MZ      = ALPHA / (1.0 - DA_TOTAL_P138)
C = math.pi * ALPHA_MZ / (math.sqrt(2.0) * G_F * M_Z2)
x = 0.5 * (1.0 + math.sqrt(1.0 - 4.0 * C))   # larger root gives M_W > M_Z*sin(theta_W)
MW_P138 = M_Z * math.sqrt(x)
MW_RESIDUAL_PCT = (MW_P138 - PDG_MW) / PDG_MW * 100.0

# ---------------------------------------------------------------------------
# Report
# ---------------------------------------------------------------------------

PASS = FAIL = 0

def check(n, desc, computed, claimed, tol_rel, unit=""):
    global PASS, FAIL
    rel_err = abs(computed - claimed) / abs(claimed) if claimed != 0 else abs(computed)
    ok = rel_err <= tol_rel
    PASS += ok; FAIL += (not ok)
    print(f"  [{'PASS' if ok else 'FAIL'}] {n:>2}. {desc}")
    print(f"        computed={computed:>12.5g}  claimed={claimed:>12.5g}  rel err={rel_err * 100:>7.3f}%")
    return ok


print("  P138 Verification: Pion FSI + rho-omega Corrections to Delta_alpha_had")
print()

# Core constants
print(f"  ALPHA_INV  = {ALPHA_INV:.10f}")
print(f"  ALPHA      = {ALPHA:.10e}")
print(f"  epsilon    = {EPS:.6f}   (claimed -0.03953)")
print(f"  S_thresh   = {S_THRESH:.6f} GeV^2   (4m_pi^2)")
print()

# Sub-interval FSI breakdown
print("  FSI sub-interval contributions (x1e-4):")
fsi_subs = []
for k in range(len(S_BREAKS_FSI) - 1):
    sub = (ALPHA / (3.0 * math.pi)) * integrate_gl(
        da_fsi_integrand, S_BREAKS_FSI[k], S_BREAKS_FSI[k+1], 400
    )
    claimed_str = [0.052, 0.151, 0.304, 0.520, 1.155][k]
    fsi_subs.append(sub)
    lo = math.sqrt(S_BREAKS_FSI[k])
    hi = math.sqrt(S_BREAKS_FSI[k+1])
    print(f"    [{lo:.3f}, {hi:.3f}] GeV: computed={sub*1e4:.3f}  claimed={claimed_str:.3f}  (x1e-4)")

print()

check(1, "epsilon (rho-omega mixing)",      EPS,        -0.03953,   5e-3)
check(2, "FSI Omnes correction (Da_FSI)",   DA_FSI,     +2.182e-4,  2e-2)
check(3, "rho-omega cross-term (Da_rw)",    DA_RW,      -2.316e-4,  2e-2)
check(4, "Net correction (Da_FSI+Da_rw)",   DA_NET,     -1.34e-5,   0.20)   # looser: near-cancellation
check(5, "P138 total (Da_had_P138)",        DA_P138,    3.718e-3,   2e-3)
check(6, "Residual gap fraction (%)",       GAP_PCT,    -9.55,      2e-2)
check(7, "W-boson mass MW (GeV)",           MW_P138,    80.383,     5e-4)

print()
print(f"  DA_FSI          = {DA_FSI:.4e}   (claimed +2.182e-4)")
print(f"  DA_RW           = {DA_RW:.4e}   (claimed -2.316e-4)")
print(f"  Net             = {DA_NET:.4e}   (claimed -1.34e-5)")
print(f"  DA_had(P138)    = {DA_P138:.6e}   (claimed 3.718e-3)")
print(f"  Gap             = {GAP_PCT:.4f}%       (claimed -9.55%)")
print(f"  MW              = {MW_P138:.6f} GeV  (claimed 80.383 GeV)")

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