# SCRIPT: BH-FIGS-V6
# Figures f13 (the ISCO and Levin's separatrix) and f14 (Bell, both routes), from the formulas of
# BH-ISCO-10 and BH-BELL-11.
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, matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy.integrate import quad
from scipy.optimize import brentq
from scipy.linalg import expm
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)
# ---------- f13 ----------
W=lambda r:(r-6)/(r**3*(r-3))
r=np.linspace(3.6,12,800); w=W(r)
grow=np.where(w<0,np.sqrt(np.abs(w)),0); osc=np.where(w>0,np.sqrt(np.abs(w)),0)
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.plot(r,grow,color=C2,lw=2.4); ax.plot(r,-grow,color=C2,lw=2.4)
ax.plot(r,osc,color=C1,lw=2.4,ls=(0,(5,3))); ax.plot(r,-osc,color=C1,lw=2.4,ls=(0,(5,3)))
ax.axvline(6,color=RULE,lw=1.2,ls=(0,(4,3)))
ax.text(6.15,0.16,'ISCO, r = 6M:\na Jordan block',color=INK,fontsize=9.5)
ax.text(3.7,0.2,'growth or decay\n(real pair)',color=INK,fontsize=9.5)
ax.text(8.6,0.05,'oscillation\n(imaginary pair,\ndashed)',color=INK,fontsize=9.5)
ax.set_xlabel(r'circular-orbit radius $r_0/M$'); ax.set_ylabel('wobble eigenvalues (real or imaginary part)')
ax.set_title('a  the ISCO is a spectral wall',color=INK,fontsize=11,loc='left')
# separatrix whirls
L2=13.0; V=lambda x:(1-2/x)*(1+L2/x**2)
ru=(L2-np.sqrt(L2**2-12*L2))/2; Es2=V(ru)
Vpp=lambda x:(V(x+1e-5)-2*V(x)+V(x-1e-5))/1e-10
lam=np.sqrt(-Vpp(ru)/2); Om=np.sqrt(L2)/ru**2
xs=[]; ys=[]
for dE in np.logspace(-3,-8,11):
    E2=Es2-dE; rs=np.linspace(ru*1.0001,60,200000); f=E2-V(rs); idx=np.where(np.diff(np.sign(f))!=0)[0]
    rp=brentq(lambda x:E2-V(x),rs[idx[0]],rs[idx[0]+1],xtol=1e-14); ra=brentq(lambda x:E2-V(x),rs[idx[1]],rs[idx[1]+1],xtol=1e-14)
    g=lambda x:(np.sqrt(L2)/x**2)/np.sqrt(max(E2-V(x),1e-300))
    I,_=quad(g,rp,ra,limit=2000); xs.append(np.log(1/dE)); ys.append(2*I/(2*np.pi))
xs=np.array(xs); ys=np.array(ys)
bx.plot(xs,ys,'o',color=C3,ms=7)
fitx=np.linspace(xs.min(),xs.max(),10); b0=np.mean(ys-Om/lam*xs/(2*np.pi))
bx.plot(fitx,Om/lam*fitx/(2*np.pi)+b0,color=C4,lw=2)
bx.set_xlabel(r'$\ln(1/\delta E)$  (closeness to the separatrix)'); bx.set_ylabel('whirls per radial period')
bx.text(xs.min(),ys.max()*0.97,r'slope $\Omega/(2\pi\lambda)$ of the unstable circular orbit',color=INK,fontsize=10)
bx.set_title("b  Levin's separatrix: whirls grow like a logarithm",color=INK,fontsize=11,loc='left')
save(fig,'f13_isco')
# ---------- f14 ----------
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)
x=np.linspace(0.05,12,500); rsq=np.arctanh(np.exp(-x/2)); chsh=2*np.sqrt(1+np.tanh(2*rsq)**2)
ax.plot(x,chsh,color=C5,lw=2.6); ax.axhline(2,color=MUT,lw=1.2,ls=(0,(4,3))); ax.axhline(2*np.sqrt(2),color=RULE,lw=1,ls=(0,(1,3)))
ax.text(7,2.05,'bound for any local hidden record',color=INK,fontsize=9.5)
ax.text(0.3,2.86,r'quantum maximum $2\sqrt{2}$',color=MUT,fontsize=9.5)
ax.text(3.2,2.55,'one Hawking pair, CHSH value',color=INK,fontsize=10)
ax.set_xlabel(r'mode frequency over Hawking temperature, $\omega/T_H$'); ax.set_ylim(1.9,2.92)
ax.set_title('a  entanglement route: the pair beats Bell at every frequency',color=INK,fontsize=10.5,loc='left')
k,a=1.0,1.0; t=np.linspace(0,12,400)
for b,col,lab in [(0.2,C1,'return rate 0.2'),(0.02,C3,'return rate 0.02'),(0.0,C2,'no return (one-way)')]:
    L=np.array([[-(k+a),b],[a,-b]]); esc=[]
    for tt in t:
        # escaped energy = 1 - x - y after an impulse x(0)=1
        v=expm(L*tt)@np.array([1.0,0.0]); esc.append(1-v.sum())
    bx.plot(t,esc,color=col,lw=2.4); bx.text(t[-1]+0.2,esc[-1],lab,color=INK,fontsize=9.5,va='center')
bx.set_xlim(0,17); bx.set_xlabel('time'); bx.set_ylabel('energy that has left through the outside')
bx.set_title('b  anomaly route: what the one-way books must carry',color=INK,fontsize=10.5,loc='left')
save(fig,'f14_bell')
print('f13, f14 written')