import sympy as sp
r = sp.symbols('r', positive=True)
psi = 2/r
# Coordinate orbital frequency Omega = dphi/dt for circular geodesic.
# phidot = L e^{-psi}/r^2, tdot = E e^{psi}. Omega = phidot/tdot = (L/E) e^{-2psi}/r^2
# Use L^2=r^2 e^{2/r}/(r-2), E^2=(r-1)e^{-2/r}/(r-2).
L2 = r**2*sp.exp(2/r)/(r-2)
E2 = (r-1)*sp.exp(-2/r)/(r-2)
Omega2 = sp.simplify((L2/E2)*sp.exp(-4/r)/r**4)  # (L/E)^2 e^{-4psi}... wait psi=2/r, e^{-2psi}=e^{-4/r}
# Omega = (L/E) e^{-2psi}/r^2 ; e^{-2psi}=e^{-4/r}
Omega2 = sp.simplify( (L2/E2) * sp.exp(-8/r)/r**4 )
print("Omega^2(r) =", Omega2)

# Radial epicyclic frequency. For metric -e^{-psi}dt^2 + e^{psi}(dr^2+r^2 dphi^2),
# rdot^2 = E^2 - Veff with Veff = L^2 e^{-2psi}/r^2 + e^{-psi} (this is (dr/dtau)^2 effectively? careful)
# Actually the eqn was: grr rdot^2 = E^2 e^{psi}... let me redo as proper radial eq.
# We had: e^{psi} rdot^2 = E^2 e^{psi} - L^2 e^{-psi}/r^2 - 1  where rdot=dr/dtau
# => rdot^2 = E^2 - e^{-psi}( L^2 e^{-psi}/r^2 + 1 ) ... let me recompute cleanly.
tau = sp.symbols('tau')
E,L = sp.symbols('E L', positive=True)
# normalization: -e^{-psi} tdot^2 + e^{psi} rdot^2 + e^{psi} r^2 phidot^2 = -1
# tdot=E e^{psi}, phidot = L e^{-psi}/r^2
expr = -sp.exp(-psi)*(E*sp.exp(psi))**2 + sp.exp(psi)*r**2*(L*sp.exp(-psi)/r**2)**2 + 1  # = -e^{psi} rdot^2
# so e^{psi} rdot^2 = -expr ; rdot^2 = -expr * e^{-psi}
rdot2 = sp.simplify(-expr*sp.exp(-psi))
print("rdot^2 (dr/dtau)^2 =", rdot2)
# = E^2 - [L^2 e^{-2psi}/r^2 + e^{-psi}]  -> matches Veff_tau = L^2 e^{-2psi}/r^2 + e^{-psi}
Vt = L**2*sp.exp(-2*psi)/r**2 + sp.exp(-psi)
# radial epicyclic (proper) freq^2: kappa_tau^2 = (1/2) d^2 Vt/dr^2 at circular
# but observed kappa relative to Omega: ratio kappa/Omega is coordinate-invariant for QPO test.
d2Vt = sp.diff(Vt, r, 2)
L2 = r**2*sp.exp(2/r)/(r-2)
d2Vt_c = sp.simplify(d2Vt.subs(L**2, L2))
print("d2Vt/dr2 at circular =", d2Vt_c)
# kappa^2 (proper) = (1/2) d2Vt; tdot=E e^{psi}, so kappa_coord^2 = kappa_tau^2 / tdot^2
# ratio (kappa/Omega)^2 = [ (1/2) d2Vt / tdot^2 ] / Omega^2
# tdot^2 = E^2 e^{2psi}
kappa_tau2 = sp.simplify(sp.Rational(1,2)*d2Vt_c)
tdot2 = E2*sp.exp(2*psi).subs(r,r)
import sympy
tdot2 = ((r-1)*sp.exp(-2/r)/(r-2))*sp.exp(4/r)
kappa_coord2 = sp.simplify(kappa_tau2 / tdot2)
ratio2 = sp.simplify(kappa_coord2/Omega2)
ratio2s = sp.simplify(ratio2)
print("(kappa/Omega)^2 =", ratio2s)
# evaluate at ISCO (should be 0 there) and general
risco=3+sp.sqrt(5)
print("(kappa/Omega)^2 at ISCO:", sp.simplify(ratio2s.subs(r,risco)))
# the simplest invariant: kappa^2/Omega^2 = 1 - r_isco-form? print symbolic
print("simplified ratio2:", sp.factor(sp.simplify(ratio2s)))
