# SCRIPT: BH-FIGS-V7
# Figures f15 (ringdown exceptional point: single excitation factors against the pair) from BH-QNM-13, and
# f11b (knots at roots of unity, now with the (3,4) torus knot 8_19) from BH-KNOT-08 and BH-KNOT-14.
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
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)
# ---------- f15 ----------
D=np.load(_os.path.join(_HERE,'bh_qnm_13.npz'))
g=1.0; w0=1-0.5j; u=np.array([1.0,0.3+0.2j]); v=np.array([0.7,-0.4+0.1j])
def modes(eta):
    d=1j*g*(1+eta); s=np.sqrt(complex(d*d+g*g)); out=[]
    for w in (w0+s,w0-s):
        r=np.array([g,w-(w0+d)]); out.append((w,(u@r)*(r@v)/(r@r)))
    return out
etas=np.logspace(-8,-0.5,120); gaps=[]; Bp=[]; Bm=[]; Bs=[]
for e in etas:
    (wp,bp),(wm,bm)=modes(e); gaps.append(abs(wp-wm)); Bp.append(abs(bp)); Bm.append(abs(bm)); Bs.append(abs(bp+bm))
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.loglog(gaps,Bp,color=C2,lw=2.4); ax.loglog(gaps,Bm,color=C4,lw=2.0,ls=(0,(5,3)))
ax.loglog(gaps,Bs,color=C1,lw=2.6)
ax.text(1.2e-2,6e3,r'each mode alone: $|B_\pm|\simeq|u\cdot Nv|/\mathrm{gap}$',color=INK,fontsize=10)
ax.text(2e-4,0.9,r'the pair: $B_++B_-=u\cdot v$, at every gap',color=INK,fontsize=10)
ax.set_ylim(0.3,3e4)
ax.set_xlabel('gap between the two overtones'); ax.set_ylabel('size of the excitation factor')
ax.set_title('a  near the exceptional point the single factors diverge',color=INK,fontsize=10.5,loc='left')
ts=D['ts']; single=D['single']; h3=D['h3']; hj=D['hj']
bx.plot(ts,single.real,color=C2,lw=1.6)
bx.plot(ts,(h3-single).real,color=C4,lw=1.6,ls=(0,(5,3)))
bx.plot(ts,h3.real,color=C1,lw=2.6)
bx.plot(ts,hj.real,color=INK,lw=1.0,ls=(0,(2,2)))
bx.set_ylim(-14,14)
bx.text(13,9.5,'overtone +, alone',color=C2,fontsize=10); bx.text(13,-11,'overtone −, alone',color=C4,fontsize=10)
bx.text(13,2.2,r'the pair: finite, $\to e^{-i\omega_0 t}(u\cdot v-i t\,u\cdot Nv)$',color=INK,fontsize=10)
bx.set_xlabel('time $t$'); bx.set_ylabel('real part of the ringdown signal')
bx.set_title(r'b  the two huge terms cancel to a finite ringdown ($\eta=10^{-3}$)',color=INK,fontsize=10.5,loc='left')
save(fig,'f15_qnm_pair')
# ---------- f11b ----------
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)
K=json.load(open(_os.path.join(_HERE,'bh_knot_14.json')))['c640']
fig=plt.figure(figsize=(12.6,4.4))
ax=fig.add_axes([0.05,0.14,0.27,0.72]); bx=fig.add_axes([0.385,0.14,0.27,0.72]); cx=fig.add_axes([0.72,0.14,0.27,0.72])
for a in (ax,bx,cx): clean(a)
ax.semilogx(Ns,raw,'o-',color=C4,lw=2,ms=5); ax.semilogx(Ns,cor,'o-',color=C1,lw=2,ms=5); ax.axhline(VOL,color=MUT,lw=1,ls=(0,(4,3)))
ax.text(30,2.62,r'$\frac{2\pi}{N}\log\langle 4_1\rangle_N$',color=INK,fontsize=11)
ax.text(30,2.12,'with $N^{3/2}$ and $3^{-1/4}$ removed',color=INK,fontsize=9.5)
ax.set_ylim(1.9,3.3); ax.text(150,1.93,'Vol$(4_1)=2.029883$',color=MUT,fontsize=9.5)
ax.set_xlabel('$N$'); ax.set_title('a  figure-eight: real, grows by its volume',color=INK,fontsize=10.5,loc='left')
bx.semilogx(Ns,Fr,'o-',color=C2,lw=2,ms=5); 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'$|\langle 3_1\rangle_N|/N^{3/2}$')
bx.set_title('b  trefoil T(2,3): one constant',color=INK,fontsize=10.5,loc='left')
for blk,col,mk in ((640,C3,'o'),(656,C5,'s'),(672,C4,'^')):
    rs=np.arange(16); cs=[K[str(blk+r)] for r in rs]
    cx.plot(rs,cs,mk,color=col,ms=7 if mk=='o' else 5,mfc='none' if mk!='o' else col,mew=1.6,label=f'N = {blk}+r')
cx.set_xlabel('r = N mod 16'); cx.set_ylabel(r'$|\langle 8_{19}\rangle_N|/N^{3/2}$'); cx.set_xticks([0,4,8,12,15])
cx.legend(frameon=False,fontsize=9,loc='upper left',labelcolor=INK)
cx.set_title(r'c  8$_{19}$ = T(3,4): sixteen constants',color=INK,fontsize=10.5,loc='left')
save(fig,'f11b_knots')