# SCRIPT: BH-QNM-13
# Two merging ringdown overtones near an exceptional point: single excitation factors against the pair.
# Model: the two-mode block of a quasinormal spectrum, complex symmetric as quasinormal modes are under
# their bilinear (unconjugated) product,
#     H(eta) = [[w0 + d, g], [g, w0 - d]],  d = i g (1 + eta),  g = 1,  w0 = 1 - 0.5 i (both modes damped for |eta| <= 0.1),
# with eigenfrequencies w0 +- sqrt(d^2 + g^2); the exceptional point is eta = 0.
# Excitation factor of mode +-:  B_+- = (u.r)(r.v)/(r.r), with r the right eigenvector and "." the
# bilinear product (no conjugation). Waveform h(t) = u.exp(-iHt)v = sum_+- B_+- exp(-i w_+- t).
# Why this tests the claim: B_+ + B_- = u.(P_+ + P_-)v, and P_+ + P_- is the pair projector, which is
# the identity on the block for every eta; so if the cancellation is real the sum is u.v exactly, while
# each B alone must grow like 1/gap. The waveform at the exceptional point must be the Jordan form
# exp(-i w0 t)(u.v - i t u.N v), N = H(0) - w0, a term linear in time and no divergence.
# Kill conditions, written before the run:
#  K1  |B_+-| * |gap| tends to |u.N v|, N = H(0) - w0 the Jordan coupling, with error of order gap
#      (single factors diverge like |coupling|/gap, the law of the wall). Kill if not.
#  K2  B_+ + B_- equals u.v to 1e-9 relative at every eta, and arg(B_+/B_-) -> pi. Kill if not.
#  K3  max over t in [0,40] of |h(t)| stays bounded (below 10) as eta -> 0; the exact h at eta = 0 equals the
#      Jordan form to 1e-9; the two-mode sum approaches it with error falling at least like the gap (observed: like gap^2). KILL CONDITION FOR THE CLAIM: the pair sum diverges.
#  K4  the contour (Riesz) integral (1/2 pi i) oint u.(z-H)^-1 v exp(-izt) dz around both modes equals h(t)
#      at eta = 0 and beside it to 1e-9. Kill if not.
#  K5  a loop eta = 0.05 exp(i phi) exchanges the two frequencies and the two excitation factors after one
#      turn and restores them after two, while the ledger log det(c - H) returns single-valued. Kill if not.
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
from scipy.linalg import expm
ok=True
def report(name,cond,val=''):
    global ok
    print(('PASS ' if cond else 'FAIL ')+name, val); ok=ok and cond
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 H(eta):
    d=1j*g*(1+eta); return np.array([[w0+d,g],[g,w0-d]],dtype=complex)
def modes(eta):
    Hm=H(eta); s=np.sqrt(complex((Hm[0,0]-w0)**2+g*g)); out=[]
    for w in (w0+s,w0-s):
        r=np.array([g,w-Hm[0,0]]) # (H-w) r = 0 from the first row
        out.append((w,(u@r)*(r@v)/(r@r)))
    return out
rows=[]
for eta in [1e-1,1e-2,1e-3,1e-4,1e-5,1e-6]:
    (wp,Bp),(wm,Bm)=modes(eta); gap=abs(wp-wm)
    rows.append((eta,gap,abs(Bp),abs(Bm),abs(Bp)*gap,Bp+Bm,np.angle(Bp/Bm)/np.pi))
    print(f'   eta={eta:.0e}: gap={gap:.3e}  |B+|={abs(Bp):.3e}  |B-|={abs(Bm):.3e}  |B+|*gap={abs(Bp)*gap:.6f}  B++B-={Bp+Bm:.6f}  arg(B+/B-)={np.angle(Bp/Bm)/np.pi:+.6f} pi')
uv=u@v
cJ=abs(u@(H(0)-w0*np.eye(2))@v)
report('K1 single excitation factors grow like |coupling|/gap', all(abs(r[4]-cJ)<r[1] for r in rows) and rows[-1][2]>1e2, f'|B+|*gap -> {rows[-1][4]:.6f} against |u.N v| = {cJ:.6f}; |B+| = {rows[-1][2]:.3e} at eta=1e-6')
report('K2 B+ + B- = u.v (pair projector) and opposite phases', all(abs(r[5]-uv)<1e-9*abs(uv) for r in rows) and abs(abs(rows[-1][6])-1)<1e-3, f'u.v = {uv:.6f}; arg(B+/B-) = {rows[-1][6]:+.6f} pi')
ts=np.linspace(0,40,401)
def h_modes(eta,t):
    return sum(B*np.exp(-1j*w*t) for w,B in modes(eta))
def h_exact(eta,t):
    return u@expm(-1j*H(eta)*t)@v
N0=H(0)-w0*np.eye(2)
jordan=lambda t: np.exp(-1j*w0*t)*(uv-1j*t*(u@N0@v))
print(f'   nilpotent check: |N0^2| = {np.abs(N0@N0).max():.1e}, u.N0.v = {u@N0@v:.6f}')
mx=[]
for eta in [1e-2,1e-4,1e-6]:
    hm=np.array([h_modes(eta,t) for t in ts]); mx.append(np.abs(hm).max())
    single=max(np.abs([modes(eta)[0][1]*np.exp(-1j*modes(eta)[0][0]*t) for t in ts]))
    print(f'   eta={eta:.0e}: max|h| = {mx[-1]:.6f}   max|single mode term| = {single:.3e}')
hj=np.array([jordan(t) for t in ts]); he=np.array([h_exact(0.0,t) for t in ts])
err=np.abs(he-hj).max()/np.abs(hj).max()
h6=np.array([h_modes(1e-6,t) for t in ts]); err6=np.abs(h6-hj).max()/np.abs(hj).max()
h4=np.array([h_modes(1e-4,t) for t in ts]); err4=np.abs(h4-hj).max()/np.abs(hj).max()
print(f'   two-mode sum vs Jordan: {err4:.2e} at eta=1e-4, {err6:.2e} at eta=1e-6 (gap ratio 10)')
report('K3 pair sum bounded through the EP; Jordan form t e^{-i w0 t} at the EP', max(mx)<10 and err<1e-9 and err4/err6>5, f'max|h| = {max(mx):.4f}; exact-vs-Jordan {err:.1e}; modes(eta=1e-6)-vs-Jordan {err6:.1e}')
def riesz(eta,t,R=0.6,M=4000):
    th=np.linspace(0,2*np.pi,M,endpoint=False); z=w0+R*np.exp(1j*th); dz=1j*R*np.exp(1j*th)*(2*np.pi/M)
    Hm=H(eta); s=0
    for zi,dzi in zip(z,dz):
        s+=(u@np.linalg.solve(zi*np.eye(2)-Hm,v))*np.exp(-1j*zi*t)*dzi
    return s/(2j*np.pi)
e4=max(abs(riesz(eta,t)-h_exact(eta,t))/abs(h_exact(eta,t)) for eta in [0.0,1e-3] for t in [0.0,5.0,20.0])
report('K4 contour integral equals the waveform, at and beside the EP', e4<1e-9, f'max rel err {e4:.1e}')
c=3.0; r0=0.05; nst=2000
def loop(turns):
    (wp,Bp),(wm,Bm)=modes(r0); cur_w=wp; cur_B=Bp; led=np.log(np.linalg.det(c*np.eye(2)-H(r0))); led0=led
    for k in range(1,turns*nst+1):
        eta=r0*np.exp(2j*np.pi*k/nst); ms=modes(eta)
        j=int(np.argmin([abs(w-cur_w) for w,_ in ms])); cur_w,cur_B=ms[j]
        d=np.log(np.linalg.det(c*np.eye(2)-H(eta)))
        led=led+np.angle(np.exp(d-led))*1j+ (d.real-led.real)
    return (wp,Bp,wm,Bm),cur_w,cur_B,led-led0
(wp,Bp,wm,Bm),w1,B1,l1=loop(1); _,w2,B2,l2=loop(2)
report('K5 one loop exchanges the pair, two restore it; ledger single-valued', abs(w1-wm)<1e-6 and abs(B1-Bm)<1e-6*abs(Bm) and abs(w2-wp)<1e-6 and abs(B2-Bp)<1e-6*abs(Bp) and abs(l1)<1e-9 and abs(l2)<1e-9,
       f'after 1 turn |w-w_-| = {abs(w1-wm):.1e}, after 2 |w-w_+| = {abs(w2-wp):.1e}; ledger change {abs(l1):.1e}, {abs(l2):.1e}')
np.savez(_os.path.join(_HERE,'bh_qnm_13.npz'),rows=np.array([[r[0],r[1],r[2],r[3]] for r in rows]),ts=ts,
         hj=hj,h6=h6,single=np.array([modes(1e-3)[0][1]*np.exp(-1j*modes(1e-3)[0][0]*t) for t in ts]),
         h3=np.array([h_modes(1e-3,t) for t in ts]))
print('ALL PASS' if ok else 'SOME FAILED')