# SCRIPT: BH-FIGS-NEW-05
# Figures f8 to f12 of v5, computed from the same formulas as BH-LIFT-05, BH-PETROV-06,
# BH-SHADOW-07 (its saved data file) and BH-KNOT-08.
import os as _os
_HERE=_os.path.dirname(_os.path.abspath(__file__))
_FIG=_os.path.join(_HERE,'..','paper','fig') if _os.path.isdir(_os.path.join(_HERE,'..','paper','fig')) else _os.path.join(_HERE,'..','fig')
import numpy as np, json, matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from matplotlib.patches import FancyArrowPatch, Circle
OUT=_FIG
BG='#0f1115'; INK='#e6e8ec'; MUT='#9aa3ae'; RULE='#2a313c'
C1='#3987e5'; C2='#d95926'; C3='#199e70'; C4='#c98500'; C5='#d55181'
plt.rcParams.update({'figure.facecolor':BG,'axes.facecolor':BG,'savefig.facecolor':BG,'text.color':INK,
  'axes.labelcolor':INK,'xtick.color':MUT,'ytick.color':MUT,'axes.edgecolor':RULE,'font.family':'serif',
  'font.serif':['TeX Gyre Pagella','URW Palladio L','Palatino','DejaVu Serif'],'mathtext.fontset':'cm','font.size':11})
def clean(ax):
    for s in ['top','right']: ax.spines[s].set_visible(False)
def save(fig,name): fig.savefig(f'{OUT}/{name}.png',dpi=260,bbox_inches='tight',pad_inches=0.12); plt.close(fig)
d0=0.4
def shapeA(w3,d,t):
    t1,t2,t3=t; w1=(1-w3+d)/2; w2=(1-w3-d)/2
    return np.array([[(t1+t2)-(t1-t2)*(w1-w2),((t1-t2)-(t1+t2-2*t3)*(w1-w2))/np.sqrt(3)],
                     [np.sqrt(3)*(t1-t2)*w3, 2*t3+w3*(t1+t2-2*t3)]],dtype=complex)
def Breg(s,d,t):
    t1,t2,t3=t; c=t1+t2-2*t3
    return np.array([[-d*(t1-t2)+t1+t2,2*s*(d*c-t1+t2)],[s*(t2-t1)/2,s*s*c+2*t3]],dtype=complex)
def gshape(w3,d):
    w1=(1-w3+d)/2; w2=(1-w3-d)/2; e1=np.array([1,-1,0])/np.sqrt(2); e2=np.array([1,1,-2])/np.sqrt(6); wv=np.array([w1,w2,w3])
    f=lambda a,b: np.sum(a*b/(4*wv)); return np.array([[f(e1,e1),f(e1,e2)],[f(e2,e1),f(e2,e2)]])

# ---------------- f8: two loops, two charts ----------------
def phase_track(Mf,path):
    M=Mf(path[0]); ev,V=np.linalg.eig(M); k=0; lam=ev[k]; v=V[:,k]
    evl,U=np.linalg.eig(M.T); u=U[:,np.argmin(abs(evl-lam))]; u=u/(u@v); v0=v.copy(); ph=[0.0]
    for p in path[1:]:
        M=Mf(p); ev,V=np.linalg.eig(M); k=np.argmin(abs(ev-lam)); lam=ev[k]; vn=V[:,k]; vn=vn*(u@v)/(u@vn)
        evl,U=np.linalg.eig(M.T); un=U[:,np.argmin(abs(evl-lam))]; un=un/(un@vn); v,u=vn,un
        r=(v0.conj()@v)/(v0.conj()@v0); ph.append(np.angle(r))
    return np.unwrap(np.array(ph))
tt=(2.0,1.0,1.3); th=np.linspace(0,2*np.pi,6000); path=0.02*np.exp(1j*th)
pg=phase_track(lambda s: shapeA(s*s,d0,tt),path); pr=phase_track(lambda s: Breg(s,d0,tt),path)
fig=plt.figure(figsize=(11.2,4.4))
ax=fig.add_axes([0.05,0.14,0.27,0.74]); clean(ax)
ax.plot(th/np.pi,pg/np.pi,color=C2,lw=2.4); ax.plot(th/np.pi,pr/np.pi,color=C1,lw=2.4)
ax.set_xlabel(r'angle around the wall point, units of $\pi$  (complex $\sigma_3$)'); ax.set_ylabel(r'transported phase / $\pi$')
ax.text(0.05,-0.93,'Gram frame: the phase runs to $-\\pi$;\nthe eigenvector returns as $-v$',color=INK,fontsize=9.5)
ax.text(0.95,-0.09,'regular frame on the cover:\nreturns as $+v$',color=INK,fontsize=9.5)
ax.set_title('a  complex loop, fixed joint data',color=INK,fontsize=11,loc='left')
# (b) real loop in (s, t3) on the cover, regular frame: eigenline arrows
bx=fig.add_axes([0.38,0.14,0.27,0.74]); clean(bx)
R=0.05; n=26; phis=np.linspace(0,2*np.pi,n,endpoint=False)
loop=np.linspace(0,2*np.pi,400)
bx.plot(R*np.cos(loop),1.3+R*np.sin(loop),color=MUT,lw=1)
bx.axvline(0,color=RULE,lw=1.2,ls=(0,(4,3)))
vprev=None
for k,ph in enumerate(np.linspace(0,2*np.pi,2000)):
    s=R*np.cos(ph); t3=1.3+R*np.sin(ph); M=Breg(s,d0,(2,1,t3)).real
    ev,V=np.linalg.eig(M); v=V[:,np.argmax(ev)]
    g=gshape(s*s,d0); J=np.diag([1/np.sqrt(2),-np.sqrt(6)*s]); gr=J.T@g@J
    L=np.linalg.cholesky(gr); y=L.T@v; y=y/np.linalg.norm(y)
    if vprev is not None and y@vprev<0: y=-y
    vprev=y
    if k==0: ystart=y.copy(); pstart=(s,t3)
    if k%(2000//n)==0 and k>0:
        x0,y0=s,t3; sc=0.016
        bx.add_patch(FancyArrowPatch((x0-sc*y[0],y0-sc*y[1]),(x0+sc*y[0],y0+sc*y[1]),arrowstyle='-|>',mutation_scale=9,color=C3 if s>=0 else C5,lw=1.4))
x0,y0=pstart; sc=0.022
bx.add_patch(FancyArrowPatch((x0,y0),(x0+sc*ystart[0],y0+sc*ystart[1]),arrowstyle='-|>',mutation_scale=12,color=INK,lw=2))
bx.add_patch(FancyArrowPatch((x0,y0),(x0+sc*vprev[0],y0+sc*vprev[1]),arrowstyle='-|>',mutation_scale=12,color=C4,lw=2))
bx.text(x0+0.012,y0-0.022,'start',color=INK,fontsize=9); bx.text(x0-0.03,y0+0.02,'end',color=INK,fontsize=9)
bx.plot([0],[1.3],'o',color=INK,ms=4)
bx.set_xlim(-0.085,0.085); bx.set_ylim(1.22,1.38); bx.set_aspect('equal')
bx.set_xlabel(r'$\sigma_3$  (left: mirror side)'); bx.set_ylabel(r'$t_3$')
bx.set_title('b  real loop through the mirror',color=INK,fontsize=11,loc='left')
bx.text(-0.082,1.225,'the eigenline comes home reversed',color=MUT,fontsize=9)
# (c) the Gram observer sees the same loop folded
cx=fig.add_axes([0.72,0.14,0.26,0.74]); clean(cx)
ph=np.linspace(0,2*np.pi,400); s=R*np.cos(ph); t3=1.3+R*np.sin(ph)
cx.plot(s**2,t3,color=C4,lw=2.4); cx.axvline(0,color=RULE,lw=1.2,ls=(0,(4,3)))
cx.plot([0],[1.3],'o',color=INK,ms=4)
cx.set_xlabel(r'$w_3=\sigma_3^2$'); cx.set_ylabel(r'$t_3$'); cx.set_xlim(-0.0004,0.003)
cx.set_title('c  the Gram observer',color=INK,fontsize=11,loc='left')
cx.text(0.00015,1.292,'the path retraces itself;\nno loop in $w_3\\geq0$ can\nencircle the wall point',color=MUT,fontsize=9)
save(fig,'f8_signs')

# ---------------- f9: Majorana stars and Petrov types ----------------
import sympy as sp
z=sp.symbols('z'); nv=sp.Matrix([1-z**2,sp.I*(1+z**2),-2*z]); Bs=[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; Bs.append(M)
cols=[[sp.Poly(sp.expand((nv.T*M*nv)[0]),z).coeff_monomial(z**k) for k in range(5)] for M in Bs]
Mapn=np.array(sp.Matrix(cols).T.evalf(),dtype=complex); Bn=[np.array(M,dtype=complex) for M in Bs]
def Qr(roots):
    p=np.poly1d([1.0])
    for r in roots: 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 mingap(Q):
    ev=np.linalg.eigvals(Q); return min(abs(ev[a]-ev[b]) for a in range(3) for b in range(a+1,3))
fig=plt.figure(figsize=(11.2,4.6))
types=[('I','four distinct stars','distinct, diagonal',[(1,40),(1,140),(1,220),(1,310)]),
       ('II','one double star','Jordan block',[(2,60),(1,180),(1,290)]),
       ('D','two double stars (Kerr)','double, diagonal',[(2,90),(2,270)]),
       ('III','one triple star','one 3x3 block',[(3,70),(1,250)]),
       ('N','one quadruple star','nilpotent',[(4,90)])]
for k,(nm,desc,jd,stars) in enumerate(types):
    a=fig.add_axes([0.02+k*0.105,0.52,0.095,0.36]); a.set_aspect('equal'); a.axis('off')
    a.add_patch(Circle((0,0),1,fill=False,ec=MUT,lw=1)); a.plot([-1,1],[0,0],color=RULE,lw=0.6)
    for m,ang in stars:
        x,y=np.cos(np.radians(ang)),np.sin(np.radians(ang))
        a.plot(x,y,'o',ms=5+3*m,color=C4 if m==1 else (C5 if nm in('II','III','N') else C3),mec=BG,mew=1.5)
        if m>1: a.text(x*1.42,y*1.42,f'{m}',color=INK,ha='center',va='center',fontsize=9)
    a.set_xlim(-1.6,1.6); a.set_ylim(-1.6,1.6)
    a.text(0,1.95,nm,ha='center',color=INK,fontsize=13); a.text(0,-2.0,desc,ha='center',color=MUT,fontsize=8.2)
    a.text(0,-2.55,jd,ha='center',color=INK,fontsize=8.5)
fig.text(0.02,0.94,'a  four Majorana stars fix the Jordan form of the Weyl operator',color=INK,fontsize=11)
fig.text(0.03,0.30,'Kerr: only $\\Psi_2=-M/(r-ia\\cos\\theta)^3$ is nonzero, so the stars sit in two\ndoubled pairs at every point (type D), and the Weyl operator has\neigenvalues $(\\lambda,\\lambda,-2\\lambda)$ with no Jordan block.',color=INK,fontsize=10,va='top')
fig.text(0.03,0.10,'Small circles: the celestial sphere seen edge-on; dot size grows with multiplicity.',color=MUT,fontsize=9)
ax=fig.add_axes([0.6,0.14,0.37,0.74]); clean(ax)
ds=np.logspace(-4,-1,30)
gII=[mingap(Qr([0.5-d/2,0.5+d/2,-1.0,2.2])) for d in ds]
gD=[mingap(Qr([1-d/2,1+d/2,-1-d/2,-1+d/2])) for d in ds]
gW=[abs(np.diff(np.linalg.eigvals(Breg(d,d0,tt)))[0]) for d in ds]
ax.loglog(ds,gII,color=C5,lw=2.4); ax.loglog(ds,gD,color=C3,lw=2.4); ax.loglog(ds,gW,color=C1,lw=2.4,ls=(0,(5,3)))
ax.text(1.4e-3,1.0e-4,'toward II: gap $\\propto\\delta$ (Jordan; the stars unfold it)',color=INK,fontsize=9.5)
ax.text(1.2e-4,3e-6,'toward D: gap $\\propto\\delta^2$ (diagonalisable)',color=INK,fontsize=9.5)
ax.text(1.2e-4,2.5e-2,'our wall on the cover:\ngap $\\propto\\sigma_3$, diagonalisable',color=INK,fontsize=9.5)
ax.set_xlabel(r'star separation $\delta$  (for the wall: $\sigma_3$)'); ax.set_ylabel('smallest eigenvalue gap')
ax.set_title('b  how each degeneracy unfolds',color=INK,fontsize=11,loc='left')
save(fig,'f9_petrov')

# ---------------- f10: Ramanujan's shadow ----------------
D=json.load(open('bh_shadow_07.json'))
fig=plt.figure(figsize=(11.2,4.4)); ax=fig.add_axes([0.06,0.14,0.4,0.74]); bx=fig.add_axes([0.57,0.14,0.4,0.74])
for a in (ax,bx): clean(a)
cols={'zeta = -1':C1,'zeta = i':C2,'zeta = e^(i pi/3)':C3}
labs={'zeta = -1':r'$q\to-1$','zeta = i':r'$q\to i$','zeta = e^(i pi/3)':r'$q\to e^{i\pi/3}$'}
for k,v in D.items():
    r=np.array(v['rows']); e=r[:,0]
    ax.semilogx(e,r[:,1],'o-',color=cols[k],lw=2,ms=6)
    ax.text(e[-1]*0.85,r[-1,1],labs[k],color=INK,fontsize=10,va='center')
    err=np.hypot(r[:,3]-v['L'][0],r[:,4]-v['L'][1])
    bx.loglog(e,err,'o-',color=cols[k],lw=2,ms=6)
    bx.text(e[0]*0.97,err[0]*1.18,labs[k]+'  limit '+{'zeta = -1':'$4$','zeta = i':'$4i$','zeta = e^(i pi/3)':'$4e^{-2\\pi i/3}$'}[k],color=INK,fontsize=9.5,va='bottom')
ax.set_xlim(1.3e-2,1.3e-4); bx.set_xlim(1.3e-2,2.2e-4); bx.set_ylim(8e-4,1.0)
ax.set_xlabel(r'$1-|q|$ (approach along the radius)'); ax.set_ylabel(r'$\log_{10}|f(q)|$  (and $|b(q)|$, the same to plotting accuracy)')
ax.set_title('a  each one blows up: $10^{597}$ at $1-|q|=3\\times10^{-4}$',color=INK,fontsize=11,loc='left')
bx.set_xlabel(r'$1-|q|$'); bx.set_ylabel(r'$|\,f-(-1)^k b-L(\zeta)\,|$')
bx.set_title('b  the pair stays finite and lands on $L(\\zeta)$',color=INK,fontsize=11,loc='left')
save(fig,'f10_shadow')

# ---------------- f11: knots at roots of unity ----------------
import mpmath as mp
VOL=2.029883212819307
Ns=np.array([25,50,100,200,400,800,1600]); raw=[]; cor=[]
for N in Ns:
    q=np.exp(2j*np.pi/N); lp=0.0; terms=[0.0]
    for k in range(1,N): lp+=2*np.log(abs(1-q**k)); terms.append(lp)
    t=np.array(terms); m=t.max(); lk=m+np.log(np.sum(np.exp(t-m)))
    raw.append(2*np.pi/N*lk); cor.append(2*np.pi/N*(lk-1.5*np.log(N)+0.25*np.log(3)))
Fr=[]
for N in Ns:
    mp.mp.dps=int(40+0.1*N); q=mp.exp(2j*mp.pi/N); s=0; p=mp.mpc(1)
    for n in range(N): s+=p; p*=(1-q**(n+1))
    Fr.append(float(abs(s))/N**1.5)
fig=plt.figure(figsize=(11.2,4.4)); ax=fig.add_axes([0.06,0.14,0.4,0.74]); bx=fig.add_axes([0.57,0.14,0.4,0.74])
for a in (ax,bx): clean(a)
ax.semilogx(Ns,raw,'o-',color=C4,lw=2,ms=6); ax.semilogx(Ns,cor,'o-',color=C1,lw=2,ms=6); ax.axhline(VOL,color=MUT,lw=1,ls=(0,(4,3)))
ax.text(30,2.6,r'$\frac{2\pi}{N}\log\langle 4_1\rangle_N$',color=INK,fontsize=11)
ax.text(30,2.09,'with the $N^{3/2}$ and $3^{-1/4}$ factors removed',color=INK,fontsize=9.5)
ax.set_ylim(1.9,3.3); ax.text(300,1.93,'Vol$(4_1)=2.029883$',color=MUT,fontsize=10)
ax.set_xlabel('$N$'); ax.set_title('a  figure-eight (amphichiral): the invariant is real and grows like its volume',color=INK,fontsize=10.5,loc='left')
bx.semilogx(Ns,Fr,'o-',color=C2,lw=2,ms=6); bx.axhline(1,color=MUT,lw=1,ls=(0,(4,3))); bx.set_ylim(0.9,1.1)
bx.set_xlabel('$N$'); bx.set_ylabel(r'$|F(e^{2\pi i/N})|/N^{3/2}$')
bx.set_title('b  trefoil (chiral): no volume, power growth; the mirror conjugates',color=INK,fontsize=10.5,loc='left')
save(fig,'f11_knots')

# ---------------- f12: three walls ----------------
fig=plt.figure(figsize=(11.2,4.8))
a1=fig.add_axes([0.05,0.16,0.27,0.66]); a2=fig.add_axes([0.385,0.16,0.27,0.66]); a3=fig.add_axes([0.72,0.16,0.27,0.66])
for a in (a1,a2,a3): clean(a)
# spectral wall: H(gamma) = sigma_x + i gamma sigma_z, eigenvalues +-sqrt(1-gamma^2)
g=np.linspace(0,1.6,800); ev=np.sqrt((1-g**2).astype(complex))
a1.plot(g,ev.real,color=C1,lw=2.4); a1.plot(g,-ev.real,color=C1,lw=2.4)
a1.plot(g,ev.imag,color=C2,lw=2.0,ls=(0,(5,3))); a1.plot(g,-ev.imag,color=C2,lw=2.0,ls=(0,(5,3)))
a1.axvline(1,color=RULE,lw=1.2,ls=(0,(4,3))); a1.set_xlabel('physical parameter',labelpad=2); a1.set_ylabel('eigenvalues')
a1.text(0.05,1.25,'real part',color=INK,fontsize=9.5); a1.text(1.12,0.95,'imaginary part',color=INK,fontsize=9.5)
a1.set_title('Spectral wall',color=INK,fontsize=12,loc='left')
a1.text(0.0,-0.24,'real pair turns complex; a Jordan block\nin every frame; loop swaps the pair.\nThe Sun’s dynamo, an optical threshold.',color=MUT,fontsize=9.5,va='top',transform=a1.transAxes)
# geometric wall: eigenvalues vs sigma3 real both sides, linear crossing
s=np.linspace(-0.6,0.6,600)
for t3 in [1.3,1.32,1.28]:
    E=np.array([np.sort(np.linalg.eigvals(Breg(x,d0,(2,1,t3))).real) for x in s])
    col=C3 if t3==1.3 else MUT; lw=2.4 if t3==1.3 else 1.0
    a2.plot(s,E[:,0],color=col,lw=lw); a2.plot(s,E[:,1],color=col,lw=lw)
a2.axvline(0,color=RULE,lw=1.2,ls=(0,(4,3))); a2.set_xlabel(r'$\sigma_3$  (mirror side on the left)')
a2.set_title('Geometric wall (the coplanar fold)',color=INK,fontsize=12,loc='left')
a2.text(0.0,-0.24,'real on both sides; a plain double point\non the cover, a cone; a real loop\nthrough the mirror reverses the eigenline.\nThe Jordan block is the Gram chart’s.',color=MUT,fontsize=9.5,va='top',transform=a2.transAxes)
# horizon: static base lapse metric, finite and even; Kerr Weyl eigenvalues D
lam=np.linspace(-0.9,0.9,400); M=1.0; hl=16*M**2*(1-lam**2)**-4
a3.plot(lam,hl/16,color=C4,lw=2.4)
a3.axvline(0,color=RULE,lw=1.2,ls=(0,(4,3)))
a3.set_xlabel(r'signed lapse $\lambda=\pm\sqrt{N}$'); a3.set_ylabel(r'$h_{\lambda\lambda}/16M^2$')
a3.set_ylim(0,8); a3.set_title('Horizon (Schwarzschild, Kerr)',color=INK,fontsize=12,loc='left')
a3.text(0.0,-0.24,'static base finite and even in the lapse;\nWeyl operator type D, diagonalisable;\nMOTS principal eigenvalue simple.',color=MUT,fontsize=9.5,va='top',transform=a3.transAxes)
save(fig,'f12_walls')
print('figures f8-f12 written')