# SCRIPT: BH-FIGS-V8
# Figures f16 (the Sun's dynamo ladder: pair births on the real axis, and one loop around the first three
# birth points) from the operator of BH-SUNLADDER-16b, and f17 (the ISCO as the inverse-cube point, and the
# Sun's acoustic cutoff) from the formulas of BH-ISCO-CUBE-17.
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, scipy.sparse as sp, scipy.sparse.linalg as sla
matplotlib.use('Agg')
import matplotlib.pyplot as plt
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)
N=120; h=np.pi/(N+1); x=np.arange(1,N+1)*h
D2=sp.diags([np.ones(N-1),-2*np.ones(N),np.ones(N-1)],[-1,0,1])/h**2
D1=sp.diags([-np.ones(N-1),np.ones(N-1)],[-1,1])/(2*h); Cc=sp.diags(np.cos(x)); K=12
def spec(a):
    A=sp.bmat([[D2,a*Cc],[D1,D2]]).tocsc().astype(complex)
    w=sla.eigs(A,k=K+4,sigma=0,which='LM',return_eigenvectors=False); return w[np.argsort(-w.real)][:K]
births=[3.8717,5.2802,7.1576,9.0850,11.0266]
# (a) real axis
al=np.linspace(2.5,12.0,700); W=np.array([spec(a) for a in al])
# (b) one loop around the first three points, adaptive tracking as in BH-SUNLADDER-16b
base=3.0; a1=(births[2]+births[3])/2; ht=0.2
L1=a1-base; Lc=np.pi*ht; T=2*ht+2*L1+Lc
def P(s):
    if s<=ht: return base-1j*s
    s-=ht
    if s<=L1: return base+s-1j*ht
    s-=L1
    if s<=Lc: return a1+ht*np.exp(1j*(-np.pi/2+s/ht))
    s-=Lc
    if s<=L1: return a1-s+1j*ht
    s-=L1
    return base+1j*(ht-s)
s=0; cur=spec(P(0)); start=cur.copy(); traj=[cur.copy()]; ds=0.01
while s<T-1e-12:
    step=min(ds,T-s)
    while True:
        new=spec(P(s+step)); Dm=np.abs(cur[:,None]-new[None,:]); asg=np.argmin(Dm,axis=1)
        good=len(set(asg))==K and all(np.sort(Dm[j])[0]<=0.25*np.sort(Dm[j])[1] for j in range(K))
        if good or step<1e-6: break
        step/=2
    cur=new[asg]; s+=step; traj.append(cur.copy()); ds=min(0.01,step*1.5)
traj=np.array(traj)
perm=[int(np.argmin(np.abs(start-c))) for c in cur]
print('permutation after one loop around the first three points:', perm)
fig=plt.figure(figsize=(11.6,4.6)); ax=fig.add_axes([0.06,0.14,0.42,0.74]); bx=fig.add_axes([0.57,0.14,0.41,0.74])
for a in (ax,bx): clean(a)
cols=[MUT,C1,C1,C2,C2,C3,C3,C4,C4,C5,C5,MUT]
for j in range(11):
    re=W[:,j].real; real=np.abs(W[:,j].imag)<1e-7*np.maximum(1,np.abs(W[:,j]))
    ax.plot(al[real],re[real],'.',ms=1.6,color=cols[j]); ax.plot(al[~real],re[~real],'.',ms=1.0,color=cols[j],alpha=0.45)
for k,b in enumerate(births):
    ax.axvline(b,color=RULE,lw=1,ls=(0,(3,3)))
    ax.text(b+0.06,-45,f'({2*k+1},{2*k+2})',color=INK,fontsize=9)
ax.set_ylim(-48,1); ax.set_xlabel(r'dynamo number $\alpha_0$'); ax.set_ylabel('real part of the leading eigenvalues')
ax.text(2.6,-38,'labels count the least-damped mode as 0',color=MUT,fontsize=9)
ax.set_title('a  each birth point joins its own pair',color=INK,fontsize=10.5,loc='left')
for j in range(1,7):
    bx.plot(traj[:,j].real,traj[:,j].imag,color=cols[j],lw=1.6)
    bx.plot(start[j].real,start[j].imag,'o',color=cols[j],ms=5)
for j in (0,7,8):
    bx.plot(traj[:,j].real,traj[:,j].imag,color=MUT,lw=1.0); bx.plot(start[j].real,start[j].imag,'o',color=MUT,ms=4)
bx.set_xlabel(r'Re $\lambda$'); bx.set_ylabel(r'Im $\lambda$')
bx.set_title('b  one loop around the first three points: three swaps',color=INK,fontsize=10.5,loc='left')
bx.text(0.02,0.04,'dots: eigenvalues at the base point $\\alpha_0=3$ (all real)\ngrey: modes not enclosed, including the fundamental',transform=bx.transAxes,color=MUT,fontsize=9)
save(fig,'f16_sunladder')
# f17
r=np.linspace(3.2,30,600); p=2+6/r
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,p,color=C1,lw=2.4); ax.axhline(2,color=MUT,lw=1,ls=(0,(4,3))); ax.axhline(3,color=C2,lw=1,ls=(0,(4,3)))
ax.axvline(6,color=RULE,lw=1.2,ls=(0,(4,3)))
ax.text(20,2.06,'Newton, p = 2',color=MUT,fontsize=10); ax.text(20,3.04,'inverse cube, p = 3',color=C2,fontsize=10)
ax.text(6.3,3.6,'ISCO, r = 6M',color=INK,fontsize=10)
ax.text(12,2.72,r'$p=2+6M/r$',color=C1,fontsize=12)
ax.set_xlabel(r'orbit radius $r/M$'); ax.set_ylabel(r'effective power $p$, from $\omega_r^2/\Omega^2=3-p$')
ax.set_title('a  the ISCO is the inverse-cube point',color=INK,fontsize=10.5,loc='left')
om=np.linspace(0.2,2.0,600); kz2=om**2-1
bx.plot(om,np.where(kz2>0,np.sqrt(np.abs(kz2)),0),color=C1,lw=2.4,ls=(0,(5,3)))
bx.plot(om,np.where(kz2<0,np.sqrt(np.abs(kz2)),0),color=C2,lw=2.4)
bx.axvline(1,color=RULE,lw=1.2,ls=(0,(4,3)))
bx.text(0.25,0.42,'evanescent:\nreflected',color=C2,fontsize=10); bx.text(1.35,0.55,'propagating:\nescapes',color=C1,fontsize=10)
bx.text(1.03,1.35,r'$\omega=\omega_{\rm ac}$: a Jordan block',color=INK,fontsize=10)
bx.set_xlabel(r'frequency $\omega/\omega_{\rm ac}$'); bx.set_ylabel(r'$|k_z|\,c/\omega_{\rm ac}$ (decay or wavenumber)')
bx.set_title("b  the Sun's acoustic cutoff is a spectral wall",color=INK,fontsize=10.5,loc='left')
save(fig,'f17_cube_cutoff')