#!/usr/bin/env python3
"""
verify_P072.py -- Addendum 72: Bottom Yukawa and theta23 PMNS.

This verifier checks the numerical claims in
72_Addendum_BottomYukawa.tex.

The theta23 CKM/PMNS arithmetic using the observed Wolfenstein A is clean.
The flagged issues are in the claimed derivation of A:
  * sin(2x)/sin^2(x) equals 2*cot(x), not 2*cos(x).
  * The paper alternates between incompatible RG factors: 0.423, 0.094, and
    0.024.
  * The final theta23 closure uses A_obs as an input, so it is not a closed
    TOE derivation of A from y_b.
  * The residual is 1.15 sigma, so statements that it is "within 1 sigma" are
    arithmetically false, though "about 1.2 sigma" is correct.
"""

import math
import sys
from pathlib import Path

PASS = FAIL = 0
_N = 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}")
    return ok


class Verifier:
    """Target-style output adapter; check/record numerics are unchanged."""

    def __init__(self, name):
        print(name)

    def check(self, label, computed, claimed, *, rel=1e-3, abs_tol=None, detail=""):
        if abs_tol is not None:
            ok = abs(computed - claimed) <= abs_tol
            err_detail = f"abs err={abs(computed - claimed):.6g}, tol={abs_tol:.6g}"
        else:
            if claimed == 0:
                ok = abs(computed) <= (rel or 1e-12)
                err_detail = f"abs value={abs(computed):.6g}, tol={rel:.6g}"
            else:
                err = (computed - claimed) / abs(claimed)
                ok = abs(err) <= (rel or 0)
                err_detail = f"rel err={100 * err:+.6g}%, tol={100 * (rel or 0):.6g}%"
        return self.record(label, ok, computed, claimed,
                           err_detail + (f"; {detail}" if detail else ""))

    def record(self, label, ok, computed="", claimed="", detail=""):
        global _N
        ok = bool(ok)
        desc, extra = label, detail
        if not ok and detail:
            desc, extra = f"{label} -- {detail}", ""
        _N += 1
        check(_N, desc, ok)
        if computed != "" or claimed != "":
            print(f"       computed: {computed}")
            print(f"       claimed : {claimed}")
        if extra:
            print(f"       {extra}")
        return ok

    def summary(self):
        print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
        return 1 if FAIL else 0


v = Verifier("P072 -- Bottom Yukawa and theta23 PMNS")

DEG = math.pi / 180
LAM = math.sin(math.pi / 14)
Y2 = math.sin(2 * math.pi / 14)
A_OBS = 0.823
PDG_TH23 = 42.20
PDG_SIGMA_TH23 = 0.83


def asind(x):
    return math.asin(x) / DEG


v.check("lambda = sin(pi/14)", LAM, 0.22252, rel=5e-6)
v.check("lambda^2", LAM**2, 0.04952, rel=1e-4)
v.check("sin(2*pi/14)", Y2, 0.43388, rel=1e-4)
v.check("2*cos(pi/14)", 2 * math.cos(math.pi / 14), 1.9499, rel=5e-5)

a_ratio = Y2 / (LAM**2)
v.check(
    "A_GUT from sin(2*pi/14)/sin^2(pi/14)",
    a_ratio,
    1.9499,
    rel=1e-3,
    detail="Expected fail: this ratio is 2*cot(pi/14) ~= 8.7626, not 2*cos(pi/14).",
)
v.check("A_GUT corrected as 2*cot(pi/14)", a_ratio, 8.769, rel=1e-3)

v.check("abstract physical A from 2*cos(pi/14)*r_b with r_b=0.423", 2 * math.cos(math.pi / 14) * 0.423, 0.826, rel=2e-3)
v.check("physical A from corrected A_GUT and r_b=0.094", a_ratio * 0.094, 0.824, rel=2e-3)
v.check("naive physical A from r_b=0.024 and y2/lambda^2", 0.024 * a_ratio, 0.210, rel=2e-3)

v.record(
    "RG reduction factor is internally consistent",
    False,
    "0.423 in abstract; 0.094 in theorem proof; 0.024 in Definition 3.1",
    "one factor",
    "Expected fail: these are different quantities but are all presented as r_b-like reductions.",
)
v.record(
    "theta23 closure uses derived A rather than observed A",
    False,
    "Section 4 explicitly adopts A_obs = 0.823 as input",
    "A derived from y_b",
    "Expected fail: the angle calculation is valid conditional on A_obs, but A is not derived in this paper.",
)

v.check("|V_cb| = A_obs*lambda^2", A_OBS * LAM**2, 0.04075, rel=5e-4)
theta23_ckm = asind(A_OBS * LAM**2)
v.check("theta23 CKM from arcsin(A lambda^2)", theta23_ckm, 2.335, rel=5e-4)

nlo_shift_deg = (LAM**2 / 2) / DEG
theta23_nlo = 45 - nlo_shift_deg
v.check("lambda^2/2 in degrees", nlo_shift_deg, 1.42, rel=2e-3)
v.check("theta23 NLO = 45 - lambda^2/2", theta23_nlo, 43.58, rel=5e-4)

theta23_toe = theta23_nlo - theta23_ckm
residual = PDG_TH23 - theta23_toe
sigma = residual / PDG_SIGMA_TH23
v.check("theta23 PMNS TOE prediction", theta23_toe, 41.245, rel=5e-4)
v.check("theta23 residual in degrees", residual, 0.955, rel=2e-3)
v.check("theta23 residual percent", 100 * residual / PDG_TH23, 2.3, rel=2e-2)
v.check("theta23 residual sigma", sigma, 1.15, rel=5e-3)
v.record(
    "theta23 residual is within the 1 sigma PDG band",
    abs(residual) <= PDG_SIGMA_TH23,
    f"{residual:.4f} deg = {sigma:.3f} sigma",
    "<= 0.83 deg",
    "Expected fail: 0.955 deg is about 1.15 sigma, not within 1 sigma.",
)

third_order = (LAM**3 / math.sqrt(2)) / DEG
v.check("suggested -lambda^3/sqrt(2) correction in degrees", third_order, 0.45, rel=1e-2)

sys.exit(v.summary())