#!/usr/bin/env python3
"""
verify_P068.py — Verification script for Addendum 68 (WolfensteinPMNS)

Checks that the two-loop PMNS corrections are internally consistent:
  - θ₁₂⁽²⁾ computed from the body formula matches the abstract claim ≈ 32.59°
  - δθ₂₃ ≈ -λ²/2 ≈ -1.42° (not -0.81° as in stale abstract)
  - θ₂₃⁽²⁾ ≈ 43.58°
  - θ₁₃⁽²⁾ ≈ 8.93°

All computed at 50 d.p. using mpmath.

Constants:
    λ = sin(π/14)    (Cabibbo angle)
    θ_C = π/14
"""

import mpmath
import sys

mpmath.mp.dps = 50

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}")

lam = mpmath.sin(mpmath.pi / 14)
lam2 = lam ** 2
theta_C = mpmath.pi / 14

print(f"Wolfenstein parameter:")
print(f"  lambda = sin(pi/14) = {float(lam):.8f}")
print(f"  lambda^2 = {float(lam2):.8f}")
print(f"  theta_C = pi/14 = {float(theta_C * 180 / mpmath.pi):.6f} deg")
print()

# ── theta12^(2) (Theorem thm:theta12-corr) ─────────────────────────────────
# sin(θ₁₂^(2)) = sin(π/4 - θ_C) * (1 + λ²/4)
sin_theta12_LO = mpmath.sin(mpmath.pi / 4 - theta_C)
sin_theta12_NLO = sin_theta12_LO * (1 + lam2 / 4)
theta12_NLO = mpmath.asin(sin_theta12_NLO)
theta12_NLO_deg = float(theta12_NLO * 180 / mpmath.pi)
print(f"theta12^(2):")
print(f"  sin(theta12_LO) = sin(pi/4 - pi/14) = {float(sin_theta12_LO):.8f}")
print(f"  correction factor (1 + lam^2/4)     = {float(1 + lam2/4):.8f}")
print(f"  sin(theta12_NLO)                     = {float(sin_theta12_NLO):.8f}")
print(f"  theta12^(2)                          = {theta12_NLO_deg:.4f} deg")

# PDG value
theta12_PDG = 33.44
residual_theta12 = theta12_PDG - theta12_NLO_deg
print(f"  PDG theta12                          = {theta12_PDG:.2f} deg")
print(f"  Residual = PDG - NLO                 = {residual_theta12:.4f} deg ({residual_theta12/theta12_PDG*100:.2f}%)")
print()

# ── delta theta23 (Proposition prop:theta23-lam2) ─────────────────────────
# δθ₂₃ = -λ²/2 / sqrt(1-λ²) ≈ -λ²/2 ≈ -1.42°
dtheta23_exact = -lam2 / 2 / mpmath.sqrt(1 - lam2)
dtheta23_approx = -lam2 / 2
dtheta23_exact_deg = float(dtheta23_exact * 180 / mpmath.pi)
dtheta23_approx_deg = float(dtheta23_approx * 180 / mpmath.pi)
theta23_NLO_deg = 45.0 + dtheta23_exact_deg
print(f"delta_theta23:")
print(f"  exact: -lam^2/(2*sqrt(1-lam^2))     = {float(dtheta23_exact):.8f} rad = {dtheta23_exact_deg:.4f} deg")
print(f"  approx: -lam^2/2                     = {float(dtheta23_approx):.8f} rad = {dtheta23_approx_deg:.4f} deg")
print(f"  theta23^(2) = 45 + delta             = {theta23_NLO_deg:.4f} deg")

theta23_PDG = 42.20
residual_theta23 = theta23_PDG - theta23_NLO_deg
print(f"  PDG theta23                          = {theta23_PDG:.2f} deg")
print(f"  Residual = PDG - NLO                 = {residual_theta23:.4f} deg ({residual_theta23/theta23_PDG*100:.2f}%)")
print()

# ── theta13^(2) (Theorem thm:theta13-corr) ────────────────────────────────
# sin(θ₁₃^(2)) = λ/√2 * (1 - λ²/4)
sin_theta13_NLO = (lam / mpmath.sqrt(2)) * (1 - lam2 / 4)
theta13_NLO = mpmath.asin(sin_theta13_NLO)
theta13_NLO_deg = float(theta13_NLO * 180 / mpmath.pi)
print(f"theta13^(2):")
print(f"  lambda/sqrt(2)                       = {float(lam/mpmath.sqrt(2)):.8f}")
print(f"  correction factor (1 - lam^2/4)      = {float(1 - lam2/4):.8f}")
print(f"  sin(theta13_NLO)                     = {float(sin_theta13_NLO):.8f}")
print(f"  theta13^(2)                          = {theta13_NLO_deg:.4f} deg")

theta13_PDG = 8.620
residual_theta13 = theta13_PDG - theta13_NLO_deg
print(f"  PDG theta13                          = {theta13_PDG:.3f} deg")
print(f"  Residual = PDG - NLO                 = {residual_theta13:.4f} deg ({residual_theta13/theta13_PDG*100:.2f}%)")
print()

# ═══════════════════════════════════════════════════════════════════════════
# ASSERTIONS
# ═══════════════════════════════════════════════════════════════════════════

# ── theta12 ──────────────────────────────────────────────────────────────
check(1, f"theta12^(2) = {theta12_NLO_deg:.4f} deg ≈ 32.59 deg",
      abs(theta12_NLO_deg - 32.59) < 0.05)

# Stale abstract value (32.78) must be rejected
check(2, f"stale value 32.78 rejected (|diff| = {abs(theta12_NLO_deg-32.78):.4f} > 0.10)",
      abs(theta12_NLO_deg - 32.78) > 0.10)

# ── delta theta23 ─────────────────────────────────────────────────────────
# approx value should be close to -1.42 deg
check(3, f"-lam^2/2 = {dtheta23_approx_deg:.4f} deg ≈ -1.42 deg",
      abs(dtheta23_approx_deg - (-1.42)) < 0.10)

# Stale abstract value (-0.81) must be rejected
check(4, f"stale delta_theta23 -0.81 rejected (|diff| = {abs(dtheta23_approx_deg-(-0.81)):.4f} > 0.40)",
      abs(dtheta23_approx_deg - (-0.81)) > 0.40)

# theta23^(2) should be close to 43.58 deg
check(5, f"theta23^(2) = {theta23_NLO_deg:.4f} deg ≈ 43.58 deg",
      abs(theta23_NLO_deg - 43.58) < 0.15)

# Stale abstract value (44.19) must be rejected
check(6, f"stale theta23^(2) 44.19 rejected (|diff| = {abs(theta23_NLO_deg-44.19):.4f} > 0.40)",
      abs(theta23_NLO_deg - 44.19) > 0.40)

# ── theta13 ──────────────────────────────────────────────────────────────
check(7, f"theta13^(2) = {theta13_NLO_deg:.4f} deg ≈ 8.93 deg",
      abs(theta13_NLO_deg - 8.93) < 0.05)

# Stale abstract value (8.78) must be rejected
check(8, f"stale theta13^(2) 8.78 rejected (|diff| = {abs(theta13_NLO_deg-8.78):.4f} > 0.10)",
      abs(theta13_NLO_deg - 8.78) > 0.10)

# ── Consistency: residuals are finite and reasonable ──────────────────────
check(9, "all residuals within reasonable bounds",
      abs(residual_theta12) < 2.0 and abs(residual_theta23) < 5.0
      and abs(residual_theta13) < 2.0)

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