#!/usr/bin/env python3
"""verify_P192.py: Spectrum of Δ_{S³} + c·sin(χ)·∂/∂φ on S³, m₁=1, m₂=0, ℓ=0..15.

Key finding: ME[ℓ,ℓ±1]=0 for all ℓ (Parity Null Theorem — integrand odd under
χ→π-χ when ℓ+ℓ' is odd). Leading physical coupling is ℓ↔ℓ±2.
"""
import numpy as np
from math import factorial, sqrt
import warnings; warnings.filterwarnings('ignore')

M1, M2, N = 1, 0, 16   # m₁=1, m₂=0, ℓ=0..15
TARGET = 206.77


def wd(j, m1, m2, b):
    """Wigner small-d element d^j_{m1,m2}(b). Returns 0 for invalid qn."""
    if abs(m1) > j or abs(m2) > j:
        return 0.0
    pf = sqrt(factorial(j+m1)*factorial(j-m1)*factorial(j+m2)*factorial(j-m2))
    cb, sb = np.cos(b/2.0), np.sin(b/2.0)
    total = 0.0
    for s in range(j+1):
        a, c_, d_ = j+m2-s, m1-m2+s, j-m1-s
        if a < 0 or c_ < 0 or d_ < 0:
            continue
        total += ((-1)**c_ * cb**(2*j+m2-m1-2*s) * sb**(m1-m2+2*s)
                  / (factorial(s)*factorial(a)*factorial(c_)*factorial(d_)))
    return pf * total


def integrate(f, n=2000):
    """Composite Simpson on [0,π]."""
    x = np.linspace(0.0, np.pi, n+1)
    y = np.array([f(xi) for xi in x])
    h = np.pi / n
    return h/3 * (y[0] + 4*np.sum(y[1::2]) + 2*np.sum(y[2:-1:2]) + y[-1])


def me(l1, l2):
    return integrate(lambda x: wd(l1,M1,M2,x)*wd(l2,M1,M2,x)*np.sin(x)**3)


# --- Parity Null Theorem check: ME[ℓ,ℓ+1]=0 for all ℓ ---
adj_max = max(abs(me(l, l+1)) for l in range(1, N-1))
assert adj_max < 1e-10, f"Parity Null Theorem violated: max|ME[l,l+1]|={adj_max}"

# --- Corrected matrix: leading coupling is ℓ↔ℓ+2 (ℓ+ℓ' even) ---
ME2 = {}
for l in range(N):
    for lp in range(l+2, min(l+3, N)):   # only |ℓ'-ℓ|=2 (leading term)
        v = me(l, lp)
        ME2[(l, lp)] = ME2[(lp, l)] = v


def build_H(c):
    H = np.diag([complex(l*(l+2)) for l in range(N)])
    for (l, lp), v in ME2.items():
        H[l, lp] += 1j * M1 * c * v
    return H


def ratio_at(c):
    ev = np.linalg.eigvals(build_H(c))
    re = np.sort(np.abs(ev.real))
    nz = re[re > 0.5]
    if len(nz) < 2:
        return None
    return nz[1] / nz[0]


# --- Scan c = 0..500 ---
best_c, best_r, best_d = 0, 0.0, np.inf
for c in range(501):
    r = ratio_at(c)
    if r is None:
        continue
    d = abs(r - TARGET)
    if d < best_d:
        best_c, best_r, best_d = c, r, d

# --- Coarse large scan c = 0..5000 (step 10) ---
best_c5, best_r5, best_d5 = best_c, best_r, best_d
for c in range(0, 5001, 10):
    r = ratio_at(c)
    if r is None:
        continue
    d = abs(r - TARGET)
    if d < best_d5:
        best_c5, best_r5, best_d5 = c, r, d

# --- Min real gap at c_opt ---
ev_opt = np.linalg.eigvals(build_H(best_c))
re_sorted = np.sort(ev_opt.real[ev_opt.real > 0.5])
min_gap = float(np.min(np.diff(re_sorted))) if len(re_sorted) > 1 else 0.0

verdict = ("CLOSED"    if best_d < 0.05 else
           "PROMISING" if best_d < 2.0  else
           "OPEN")

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

check(1, f"Parity Null Theorem: max|ME[l,l+1]| = {adj_max:.2e}", adj_max < 1e-10)
print(f"  c_opt={best_c}  ratio={best_r:.4f}  target={TARGET}  diff={best_d:.4f}")
print(f"  large_scan: c={best_c5}  ratio={best_r5:.4f}  diff={best_d5:.4f}")
print(f"  min_gap={min_gap:.4f}")
print(f"  verdict={verdict}")

print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")