# SCRIPT: BH-ISCO-CUBE-17
# Four short symbolic checks used in v8, each from a stated equation.
#  C1  Inverse-cube reading of the ISCO. For a Newtonian central force F ~ r^-p, a nudged circular orbit has
#      radial epicyclic and orbital frequencies with omega_r^2/Omega^2 = 3 - p, and circular orbits are stable
#      iff p < 3 (apsidal angle pi/sqrt(3-p)). For Schwarzschild, take omega_r^2 = V''/2 from eq. (isco) of the
#      paper (proper time, M = 1) and Omega = dphi/dtau = L/r^2 with L^2 = r^2/(r-3). Define p by the same ratio.
#      Claim: 3 - p = 1 - 6M/r exactly, so p = 2 + 6M/r: Newton's 2 far away, the inverse-cube value 3 at r = 6M.
#      (The ratio is the same in coordinate time, since both frequencies carry the same dtau/dt.)
#  C2  A clock gradient is not curvature: for ds^2 = -N(x)^2 dt^2 + dx^2 the scalar curvature is R = -2 N''/N,
#      so a lapse linear in proper distance is flat (Rindler). Schwarzschild: near r = 2M the lapse is
#      rho/(4M) = kappa*rho in proper distance rho.
#  C3  Extremal Reissner-Nordstrom: proper distance to the horizon diverges logarithmically and the lapse is
#      exponential in proper distance: no fold, no surface gravity.
#  C4  Lamb's isothermal atmosphere: vertical velocity obeys w_tt = c^2 w_zz - (c^2/H) w_z. With
#      w = exp(z/2H) psi exp(-i omega t), psi'' + k_z^2 psi = 0, k_z^2 = (omega^2 - omega_ac^2)/c^2,
#      omega_ac = c/(2H). The first-order generator [[0,1],[-k_z^2,0]] squares to -k_z^2 I and is a nonzero
#      nilpotent (a Jordan block) at omega = omega_ac.
# Kill: any identity fails to simplify to zero.
import sympy as sp
ok=True
def report(name,cond,val=''):
    global ok
    print(('PASS ' if cond else 'FAIL ')+name, val); ok=ok and cond
r,M,x,t,z,c,H,om=sp.symbols('r M x t z c H omega',positive=True)
# C1 (M = 1 units of eq. isco, restored by r -> r/M)
half_Vpp=(r-6)/(r**3*(r-3)); L2=r**2/(r-3); Om2=L2/r**4
ratio=sp.simplify(half_Vpp/Om2)
p=sp.simplify(3-ratio)
report('C1 omega_r^2/Omega^2 = 1 - 6M/r, so p = 2 + 6M/r', sp.simplify(ratio-(1-6/r))==0 and sp.simplify(p-(2+6/r))==0 and p.subs(r,6)==3,
       f'ratio = {sp.factor(ratio)}, p = {sp.simplify(p)}, p(6M) = {p.subs(r,6)}')
# C2 2D curvature of -N^2 dt^2 + dx^2
N=sp.Function('N')(x)
g=sp.diag(-N**2,1); coords=[t,x]; ginv=g.inv()
def Gam(a,b,cc): return sum(ginv[a,d]*(sp.diff(g[d,b],coords[cc])+sp.diff(g[d,cc],coords[b])-sp.diff(g[b,cc],coords[d]))/2 for d in range(2))
def Riem(a,b,cc,d):
    return sp.diff(Gam(a,b,d),coords[cc])-sp.diff(Gam(a,b,cc),coords[d])+sum(Gam(a,cc,e)*Gam(e,b,d)-Gam(a,d,e)*Gam(e,b,cc) for e in range(2))
Ric=sp.Matrix(2,2,lambda b,d: sum(Riem(a,b,a,d) for a in range(2)))
R=sp.simplify(sum(ginv[b,d]*Ric[b,d] for b in range(2) for d in range(2)))
rho=sp.symbols('rho',positive=True); kap=sp.symbols('kappa',positive=True)
report('C2 R = -2 N\'\'/N; linear lapse is flat', sp.simplify(R+2*sp.diff(N,x,2)/N)==0 and sp.simplify(R.subs(N,kap*x).doit())==0, f'R = {R}')
eps=sp.symbols('epsilon',positive=True)
lapse=sp.sqrt(1-2*M/(2*M+eps)); prop=sp.integrate(1/sp.sqrt(1-2*M/(2*M+sp.Symbol('u',positive=True))),(sp.Symbol('u',positive=True),0,eps))
lead=sp.limit(lapse/sp.sqrt(eps),eps,0); propl=sp.limit(prop/sp.sqrt(eps),eps,0)
report('C2 Schwarzschild lapse = rho/(4M) near the horizon', sp.simplify(lead/propl-1/(4*M))==0, f'lapse ~ {lead} sqrt(eps), rho ~ {propl} sqrt(eps)')
# C3 extremal RN: N = 1 - M/r, proper distance int dr/N
u=sp.symbols('u',positive=True)
prop3=sp.integrate(1/(1-M/(M+u)),(u,eps,1))
report('C3 extremal RN: proper distance diverges like -M log(eps); lapse ~ eps/M ~ exp(-rho/M)', sp.limit(prop3/(-M*sp.log(eps)),eps,0)==1, f'distance ~ {sp.limit(prop3/(-sp.log(eps)),eps,0)} * log(1/eps)')
# C4 Lamb
psi=sp.Function('psi')(z)
w=sp.exp(z/(2*H))*psi
eq=sp.expand(sp.simplify((-om**2*w - c**2*sp.diff(w,z,2) + (c**2/H)*sp.diff(w,z))*sp.exp(-z/(2*H))/(-c**2)))
target=sp.diff(psi,z,2)+(om**2-(c/(2*H))**2)/c**2*psi
kz2=sp.symbols('k_z2'); Mgen=sp.Matrix([[0,1],[-kz2,0]])
report('C4 Lamb: psi\'\' + (omega^2 - omega_ac^2)/c^2 psi = 0, omega_ac = c/2H; generator nilpotent at the cutoff',
       sp.simplify(eq-target)==0 and Mgen**2==-kz2*sp.eye(2) and (Mgen.subs(kz2,0))**2==sp.zeros(2) and Mgen.subs(kz2,0)!=sp.zeros(2), '')
print('ALL PASS' if ok else 'SOME FAILED')