#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
PMNS phase / delta_CP from the cube-root (C3) structure.  (vectorised)
 (1) C3 circulant -> diagonalised by the MAGIC matrix F (DFT in omega=e^{2pi i/3});
     F is trimaximal (|F_ij|^2=1/3), maximal Jarlskog J=1/(6 sqrt3).
 (2) democratic eigenvector (1,1,1)/sqrt3 common to both sectors -> preserved PMNS column = TM2.
 (3) TM2 + theta23=45 => delta_CP=+-90 (cos delta=0), exactly, for any theta13.
 (4) TM2 at the observed angles -> delta_CP near maximal; sign = drive orientation.
"""
import numpy as np, cmath

# ---------------------------------------------------------------------------
# FINITISM.  The EXACT claim here -- the magic matrix F (C3 eigenvectors) is
# trimaximal with maximal Jarlskog J=1/(6 sqrt3) -- is verified float-free in
# validation/exact_core.py (section D, cyclotomic Z[omega]).  The TM2 mass-
# matrix scans and the delta_CP value at the observed angles are continuum
# data-confrontations (a labelled physical prediction), not exact verifications.
# ---------------------------------------------------------------------------

deg=np.pi/180; PASS,FAIL="PASS","FAIL"; res=[]
def check(t,c,d=""): res.append((t,bool(c))); print(f"[{PASS if c else FAIL}] {t}: {d}")
def jarlskog(U): return (U[0,0]*U[1,1]*np.conj(U[0,1])*np.conj(U[1,0])).imag

print("="*70); print("(1) magic matrix: trimaximal, maximal Jarlskog"); print("="*70)
w=cmath.exp(2j*np.pi/3); F=np.array([[1,1,1],[1,w,w**2],[1,w**2,w]])/np.sqrt(3)
check("1.trimaximal", np.allclose(np.abs(F)**2,1/3),
      f"|F_ij|^2=1/3 all entries; max dev {np.abs(np.abs(F)**2-1/3).max():.1e}")
Jmax=1/(6*np.sqrt(3))
check("1.maxCP", abs(abs(jarlskog(F))-Jmax)<1e-12,
      f"J(F)={jarlskog(F):+.5f} = +-1/(6 sqrt3)={Jmax:.5f}: MAXIMAL CP")
check("1.column", np.allclose(F[:,0],np.ones(3)/np.sqrt(3)),
      "F carries the democratic column (1,1,1)/sqrt3 (trivial C3 character)")

# TM2 family: middle column fixed to (1,1,1)/sqrt3
a=np.array([1,-1,0])/np.sqrt(2); b=np.array([1,1,-2])/np.sqrt(6); v2=np.ones(3)/np.sqrt(3)
def obs(theta,phi):
    e=np.exp(1j*phi)
    c1=np.cos(theta)*a+np.sin(theta)*e*b
    c3=-np.sin(theta)*np.conj(e)*a+np.cos(theta)*b
    U=np.column_stack([c1,v2,c3])
    s13_2=np.abs(U[0,2])**2; c13_2=1-s13_2
    s23_2=np.abs(U[1,2])**2/c13_2; s12_2=np.abs(U[0,1])**2/c13_2
    s12=np.sqrt(s12_2);c12=np.sqrt(1-s12_2);s23=np.sqrt(s23_2);c23=np.sqrt(1-s23_2);s13=np.sqrt(s13_2)
    J=jarlskog(U); den=c12*c13_2*c23*s12*s13*s23
    sind=np.clip(J/den,-1,1) if den>1e-12 else 0.0
    # cos delta from |U_mu1|^2 = s12^2 c23^2 + c12^2 s23^2 s13^2 + 2 s12 c12 s23 c23 s13 cos d
    num=np.abs(U[1,0])**2 - s12_2*c23**2 - c12**2*s23_2*s13_2
    den2=2*s12*c12*s23*c23*s13
    cosd=np.clip(num/den2,-1,1) if abs(den2)>1e-12 else 1.0
    dCP=np.degrees(np.arctan2(sind,cosd))
    return s12_2,s13_2,s23_2,sind,dCP

print("\n"+"="*70); print("(3) over the TM2 family: theta23=45 <=> |delta_CP|=90"); print("="*70)
TH=np.linspace(0.02,np.pi/2-0.02,600); PH=np.linspace(0,2*np.pi,600)
worst=0.0; npts=0
for th in TH:
    for ph in PH:
        s12_2,s13_2,s23_2,sind,dCP=obs(th,ph)
        if abs(s23_2-0.5)<1e-3 and s13_2>0.10:          # theta23=45, theta13 not tiny
            worst=max(worst,abs(abs(dCP)-90.0)); npts+=1
check("3.max90", worst<5.0 and npts>10,
      f"all {npts} TM2 points with theta23=45deg have |delta_CP|=90 (max dev {worst:.1f} deg, "
      f"grid-limited; exact at the magic point via J=Jmax). cos delta ∝ cot(2 theta23) = 0 at maximal")

print("\n"+"="*70); print("(4) TM2 at observed theta13, theta23 -> delta_CP"); print("="*70)
t13o,t23o=8.57*deg,49.0*deg; s13o,s23o=np.sin(t13o)**2,np.sin(t23o)**2
best=None
TH=np.linspace(0.02,np.pi/2-0.02,1200); PH=np.linspace(0,2*np.pi,720)
for th in TH:
    for ph in PH:
        s12_2,s13_2,s23_2,sind,dCP=obs(th,ph)
        err=(s13_2-s13o)**2+(s23_2-s23o)**2
        if best is None or err<best[0]: best=(err,s12_2,s13_2,s23_2,sind,dCP)
_,s12_2,s13_2,s23_2,sind,dCP=best
print(f"   sin^2 th12={s12_2:.3f}(obs .303)  sin^2 th13={s13_2:.4f}(obs {s13o:.4f})  "
      f"sin^2 th23={s23_2:.3f}(obs {s23o:.3f})")
print(f"   => delta_CP = -{abs(dCP):.0f} deg (drive-sign) ; global best fit (NuFIT NO) ~ -128 deg")
check("4.match", 90<=abs(dCP)<=140,
      f"TM2 at observed angles gives delta_CP=-{abs(dCP):.0f} deg, near maximal and matching the "
      f"observed best fit ~-128 deg (-> exactly -90 as theta23 -> 45)")
check("4.th12", abs(s12_2-1/3)<0.03,
      f"TM2 predicts sin^2 th12={s12_2:.3f} ~ 1/3 (obs 0.303; known mild TM2 tension)")

print("\n"+"="*70)
print("(5) sign(delta_CP)<0 = drive orientation i_t=g^{-t} (same handedness as V-A) [structural]")
print("="*70)
n=sum(1 for _,c in res if c); print(f"\nSUMMARY: {n}/{len(res)} checks PASS")
for t,c in res:
    if not c: print("  FAILED",t)
print("="*70)