# SCRIPT: BH-PETROV-06
# Majorana stars and Petrov types. A Weyl spinor is a spin-2 object: a binary quartic psi(z) whose
# four roots on the Riemann sphere are its principal null directions (its Majorana stars), or a
# complex traceless symmetric 3x3 matrix Q (the Weyl operator on self-dual bivectors). The link used
# here is psi_Q(z) = n(z)^T Q n(z) with the null vector n(z) = (1 - z^2, i(1 + z^2), -2z), n.n = 0,
# so the trace part of Q drops out and Q -> psi_Q is linear on the five traceless directions.
# Kill conditions, written before the run:
#  P1  the map Q -> coefficients of psi_Q is invertible on traceless symmetric Q (nonzero determinant).
#  P2  star multiplicities fix the Jordan structure of Q: (1111) distinct eigenvalues, diagonalisable;
#      (211) type II, one 2x2 Jordan block; (22) type D, diagonalisable with eigenvalues (l,l,-2l);
#      (31) type III, one 3x3 nilpotent block; (4) type N, nilpotent with Q^2 = 0, Q != 0.
#      Kill if any type gives a different Jordan structure.
#  P3  unfolding: split one double star of a type II quartic by delta; the gap of the colliding pair
#      is LINEAR in delta (star separation) and goes like the SQUARE ROOT of the quartic coefficient
#      eps = delta^2/4 (the fold again: stars are the signed coordinate). Kill if gap/delta not constant.
#  P4  split both double stars of a type D quartic symmetrically by delta: the double eigenvalue splits
#      QUADRATICALLY, gap/delta^2 constant, because D is diagonalisable. Kill if not.
#  P5  Kerr: with only Psi_2 = -M/(r - i a cos th)^3 nonzero, Q is type D at every point, eigenvalues
#      in ratio (1,1,-2). Kill if not.
import numpy as np, sympy as sp
ok=True
def report(name,cond,val=''):
    global ok
    print(('PASS ' if cond else 'FAIL ')+name, val); ok=ok and cond
z=sp.symbols('z')
n=sp.Matrix([1-z**2, sp.I*(1+z**2), -2*z])
assert sp.expand((n.T*n)[0])==0
# traceless symmetric basis
B=[sp.diag(1,-1,0), sp.diag(1,1,-2)]
for (i,j) in [(0,1),(0,2),(1,2)]:
    M=sp.zeros(3); M[i,j]=M[j,i]=1; B.append(M)
cols=[]
for M in B:
    p=sp.Poly(sp.expand((n.T*M*n)[0]),z); c=[p.coeff_monomial(z**k) for k in range(5)]; cols.append(c)
Map=sp.Matrix(cols).T
det=sp.simplify(Map.det())
report('P1 map Q -> quartic invertible', det!=0, f'det = {det}')
Mapn=np.array(Map.evalf(),dtype=complex); Bn=[np.array(M,dtype=complex) for M in B]
def Q_from_roots(roots, lead=1.0):
    # quartic with given finite roots (a root at infinity is entered as None: degree drops)
    p=np.poly1d([lead])
    for r in roots:
        if r is not None: p=p*np.poly1d([1,-r])
    c=np.zeros(5,dtype=complex); co=p.coeffs[::-1]; c[:len(co)]=co
    x=np.linalg.solve(Mapn,c); return sum(xi*Mi for xi,Mi in zip(x,Bn))
def jordan_sizes(Q,tol=1e-7):
    ev=np.linalg.eigvals(Q); out=[]
    used=np.zeros(3,bool)
    for k in range(3):
        if used[k]: continue
        grp=[j for j in range(3) if abs(ev[j]-ev[k])<1e-5*max(1,abs(ev[k]))]
        for j in grp: used[j]=True
        lam=np.mean(ev[grp]); m=len(grp); N=Q-lam*np.eye(3)
        ranks=[3]+[np.linalg.matrix_rank(np.linalg.matrix_power(N,p),tol=tol) for p in range(1,m+1)]
        # sizes from rank drops
        d=[ranks[p-1]-ranks[p] for p in range(1,m+1)]+[0]
        sizes=[]
        for p in range(1,m+1):
            cnt=d[p-1]-d[p]; sizes+= [p]*cnt
        out.append((np.round(lam,6),sorted(sizes,reverse=True)))
    return out
cases={'I (1111)':[0.3,-1.1,2.0,0.7],'II (211)':[0.5,0.5,-1.0,2.2],'D (22)':[1.0,1.0,-1.0,-1.0],
       'III (31)':[0.4,0.4,0.4,-1.3],'N (4)':[0.6,0.6,0.6,0.6]}
expect={'I (1111)':[[1],[1],[1]],'II (211)':[[2],[1]],'D (22)':[[1,1],[1]],'III (31)':[[3]],'N (4)':[[2,1]]}
for name,roots in cases.items():
    Q=Q_from_roots(roots); js=jordan_sizes(Q); sizes=sorted([s for _,s in js])
    good=sizes==sorted(expect[name])
    extra=''
    if name.startswith('D'):
        ev=sorted(np.linalg.eigvals(Q),key=lambda x:abs(x)); r=ev[2]/ev[0]
        good=good and abs(r+2)<1e-9; extra=f' eigenvalue ratio {np.round(r,9)}'
    if name.startswith('N'):
        good=good and np.abs(Q@Q).max()<1e-9 and np.abs(Q).max()>1e-6
    report(f'P2 type {name}: Jordan structure', good, str(js)+extra)
# P3 unfolding toward II
print('P3  split a double star of type II by delta')
r3=[]
for dlt in [0.1,0.01,0.001,0.0001]:
    Q=Q_from_roots([0.5-dlt/2,0.5+dlt/2,-1.0,2.2]); ev=np.linalg.eigvals(Q)
    i,j=min(((a,b) for a in range(3) for b in range(a+1,3)),key=lambda ab:abs(ev[ab[0]]-ev[ab[1]]))
    gap=abs(ev[i]-ev[j]); print(f'    delta={dlt:g}: gap={gap:.6e}  gap/delta={gap/dlt:.6f}  gap/sqrt(eps)={gap/np.sqrt(dlt**2/4):.6f}'); r3.append(gap/dlt)
report('P3 gap linear in star separation (II)', abs(r3[-1]-r3[-2])<1e-6*abs(r3[-1]) and r3[-1]>1e-3, f'gap/delta -> {r3[-1]:.6f}')
print('P4  split both double stars of type D by delta')
vals=[]
for dlt in [0.1,0.01,0.001]:
    Q=Q_from_roots([1-dlt/2,1+dlt/2,-1-dlt/2,-1+dlt/2]); ev=np.linalg.eigvals(Q)
    i,j=min(((a,b) for a in range(3) for b in range(a+1,3)),key=lambda ab:abs(ev[ab[0]]-ev[ab[1]]))
    gap=abs(ev[i]-ev[j]); vals.append(gap/dlt**2); print(f'    delta={dlt:g}: gap={gap:.6e}  gap/delta^2={gap/dlt**2:.6f}  gap/delta={gap/dlt:.6f}')
report('P4 gap quadratic in star separation (D)', abs(vals[-1]-vals[-2])<1e-3*abs(vals[-1]), f'gap/delta^2 -> {vals[-1]:.6f}')
# P5 Kerr
M_,a_=1.0,0.7
for r_,th in [(2.5,0.4),(1.9,1.2),(6.0,2.8)]:
    psi2=-M_/(r_-1j*a_*np.cos(th))**3
    # psi = 6 Psi2 z^2 in the Kinnersley-type gauge: stars at 0,0,inf,inf
    c=np.zeros(5,dtype=complex); c[2]=6*psi2; x=np.linalg.solve(Mapn,c); Q=sum(xi*Mi for xi,Mi in zip(x,Bn))
    js=jordan_sizes(Q); ev=sorted(np.linalg.eigvals(Q),key=abs)
    report(f'P5 Kerr r={r_}, th={th}: type D, ratio (1,1,-2)', sorted([s for _,s in js])==[[1],[1,1]] and abs(ev[2]/ev[0]+2)<1e-9, str(js))
print('ALL PASS' if ok else 'SOME FAILED')