# SCRIPT: BH-ISCO-10
# Schwarzschild test-particle orbits, G = c = M = 1, proper time tau.
# Radial motion: (dr/dtau)^2 + V(r) = E^2, V = (1 - 2/r)(1 + L^2/r^2), so r'' = -V'(r)/2.
# Small wobble about a circular orbit r0: d/dtau (dr, dr') = M (dr, dr'), M = [[0, 1], [-V''(r0)/2, 0]].
# Kill conditions, written before the run:
#  I1  along the family of circular orbits, V''(r0) changes sign at r0 = 6 (the ISCO) with a simple zero:
#      eigenvalues of M are +-i w (oscillation) for r0 > 6 and +-lam (growth/decay) for r0 < 6.
#  I2  at r0 = 6 (L^2 = 12) M is a nonzero nilpotent matrix: a Jordan block.
#  I3  the unfolding coordinate k = sqrt(-V''/2) is real on one side and imaginary on the other, and the gap
#      |lam+ - lam-| / sqrt(|r0 - 6|) tends to a constant: the spectral-wall signature of v5 Section 7.
#  I4  separatrix: for L^2 = 13, orbits with energy just below the unstable circular orbit's energy advance in
#      azimuth per radial period by Dphi = (Omega/lam) ln(1/dE) + const, where Omega and lam are the unstable
#      orbit's angular rate and growth rate. Kill if the fitted slope differs from Omega/lam by more than 2 percent.
import sympy as sp, numpy as np
from scipy.integrate import solve_ivp
ok=True
def report(name,cond,val=''):
    global ok
    print(('PASS ' if cond else 'FAIL ')+name, val); ok=ok and cond
r,L2=sp.symbols('r L2',positive=True)
V=(1-2/r)*(1+L2/r**2)
Vp=sp.diff(V,r); Vpp=sp.diff(V,r,2)
# circular orbits: Vp=0 -> L2 = r^2/(r-3)
L2c=sp.solve(sp.Eq(Vp,0),L2)[0]; print('circular orbits: L^2 =',sp.simplify(L2c))
W=sp.simplify(Vpp.subs(L2,L2c)/2)       # M = [[0,1],[-W,0]]
print('V\'\'/2 on the circular family =',sp.factor(W))
zero=sp.solve(sp.Eq(sp.numer(sp.together(W)),0),r); print('zeros of V\'\'/2:',zero)
report('I1 simple sign change of V\'\' at r0 = 6', 6 in zero and sp.diff(sp.numer(sp.together(W)),r).subs(r,6)!=0)
for r0 in [12.0,7.0,6.1,5.9,5.0,4.0]:
    w=float(W.subs(r,r0)); ev=np.linalg.eigvals(np.array([[0,1],[-w,0]]))
    print(f'   r0={r0:5.2f}: V\'\'/2 = {w:+.6e}  eigenvalues = {ev[0]:.5f}, {ev[1]:.5f}')
M6=np.array([[0,1],[-float(W.subs(r,6)),0]])
report('I2 Jordan block at the ISCO', np.allclose(M6@M6,0) and np.abs(M6).max()>0, str(M6.tolist()))
rat=[]
for d in [1e-2,1e-3,1e-4,1e-5]:
    for s in [+1,-1]:
        w=float(W.subs(r,6+s*d)); ev=np.linalg.eigvals(np.array([[0,1],[-w,0]]))
        rat.append((s,d,abs(ev[0]-ev[1])/np.sqrt(d), ev[0].real!=0))
for s,d,g,real in rat: print(f'   side {"outside" if s>0 else "inside "} d={d:g}: gap/sqrt(d) = {g:.6f}  eigenvalues {"real" if real else "imaginary"}')
outs=[g for s,d,g,_ in rat if s>0]; ins=[g for s,d,g,_ in rat if s<0]
report('I3 square-root law, imaginary outside and real inside', abs(outs[-1]-outs[-2])<1e-4 and abs(ins[-1]-ins[-2])<1e-4 and all(not rl for s,d,g,rl in rat if s>0) and all(rl for s,d,g,rl in rat if s<0), f'outside -> {outs[-1]:.6f}, inside -> {ins[-1]:.6f}')
# I4 separatrix whirl law
L2v=13.0
ru=float(min(sp.solve(sp.Eq(Vp.subs(L2,L2v),0),r)))
Vf=sp.lambdify(r,V.subs(L2,L2v)); Vpf=sp.lambdify(r,Vp.subs(L2,L2v))
Es2=Vf(ru); lam=np.sqrt(-float(Vpp.subs({L2:L2v,r:ru}))/2); Om=np.sqrt(L2v)/ru**2   # dphi/dtau = L/r^2
print(f'unstable circular orbit r_u = {ru:.6f}, E_s^2 = {Es2:.9f}, lam = {lam:.6f}, Omega = {Om:.6f}, Omega/lam = {Om/lam:.6f}')
def dphi_per_period(dE):
    E2=Es2-dE
    # turning points: periapsis just outside r_u region? bound orbit between r_p > r_u and r_a; find roots of E2-V
    rs=np.linspace(ru*1.0001,200,400000); f=E2-Vf(rs); idx=np.where(np.diff(np.sign(f))!=0)[0]
    rp,ra=rs[idx[0]],rs[idx[1]]
    from scipy.optimize import brentq
    rp=brentq(lambda x:E2-Vf(x),rs[idx[0]],rs[idx[0]+1],xtol=1e-14); ra=brentq(lambda x:E2-Vf(x),rs[idx[1]],rs[idx[1]+1],xtol=1e-14)
    # dphi/dr = (L/r^2)/sqrt(E2 - V); integrate with substitution to handle endpoint sqrt singularities
    from scipy.integrate import quad
    g=lambda x:(np.sqrt(L2v)/x**2)/np.sqrt(max(E2-Vf(x),1e-300))
    I,_=quad(g,rp,ra,limit=2000,epsabs=1e-12,epsrel=1e-12)
    return 2*I, rp, ra
xs=[]; ys=[]
for dE in [1e-3,1e-4,1e-5,1e-6,1e-7,1e-8]:
    dp,rp,ra=dphi_per_period(dE); xs.append(np.log(1/dE)); ys.append(dp)
    print(f'   dE={dE:g}: periapsis {rp:.6f}, apoapsis {ra:.4f}, azimuth per radial period {dp:.6f}  (whirls ~ {dp/(2*np.pi):.3f})')
slope=(ys[-1]-ys[-3])/(xs[-1]-xs[-3])
report('I4 whirl law: slope of azimuth vs ln(1/dE) equals Omega/lam', abs(slope/(Om/lam)-1)<0.02, f'slope {slope:.6f} vs Omega/lam {Om/lam:.6f}')
print('ALL PASS' if ok else 'SOME FAILED')