import sympy as sp
r = sp.symbols('r', positive=True)
# Units G=M=c=1, so u = 1/r, psi = 2/r.
# Full optical 4-metric: ds^2 = -e^{-psi} dt^2 + e^{psi}(dr^2 + r^2 dphi^2)   (c=1)
psi = 2/r
gtt = -sp.exp(-psi)
grr = sp.exp(psi)
gphph = sp.exp(psi)*r**2

E,L = sp.symbols('E L', positive=True)
# Conserved: E = -gtt * tdot = e^{-psi} tdot ; L = gphph phidot = e^{psi} r^2 phidot
# Timelike normalization: gtt tdot^2 + grr rdot^2 + gphph phidot^2 = -1
# => -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
# => -e^{-psi}(E e^{psi})^2 + e^{psi} rdot^2 + e^{psi} r^2 (L e^{-psi}/r^2)^2 = -1
# => -E^2 e^{psi} + e^{psi} rdot^2 + L^2 e^{-psi}/r^2 = -1
# => e^{psi} rdot^2 = E^2 e^{psi} - L^2 e^{-psi}/r^2 - 1
# => rdot^2 = E^2 - [L^2 e^{-2psi}/r^2 + e^{-psi}] = E^2 - Veff
Veff = L**2*sp.exp(-2*psi)/r**2 + sp.exp(-psi)
Veff = sp.simplify(Veff)
print("Veff =", Veff)

# Circular orbits: dVeff/dr = 0 gives L^2(r); stable boundary (ISCO): also d2Veff/dr2=0
dV = sp.diff(Veff, r)
d2V = sp.diff(Veff, r, 2)
# Solve dV=0 for L^2
L2sol = sp.solve(sp.Eq(dV,0), L**2)
print("L^2(r) circular:", [sp.simplify(s) for s in L2sol])
L2 = sp.simplify(L2sol[0])
# ISCO: substitute L^2 into d2V=0
d2V_sub = sp.simplify(d2V.subs(L**2, L2))
print("d2V at circular (set=0):", d2V_sub)
isco = sp.solve(sp.Eq(sp.numer(sp.together(d2V_sub)),0), r)
print("ISCO candidates:", isco)
for s in isco:
    sv = sp.nsimplify(s, [sp.sqrt(5)])
    print("   r_isco =", s, "=", sp.simplify(sv), " approx", sp.N(s))
print("golden (3+sqrt5):", sp.N(3+sp.sqrt(5)))
print("golden phi^2 = ", sp.N((1+sp.sqrt(5))/2)**2, " ; 3+sqrt5 = 2*phi^2? ", sp.N(2*((1+sp.sqrt(5))/2)**2))
