import sympy as sp
r = sp.symbols('r', positive=True)
psi = 2/r
# From before: rdot^2 = E^2 - Veff, Veff = L^2 e^{-2psi}/r^2 + e^{-psi}
# Circular: rdot=0 and dVeff/dr=0.
# L^2(r) = r^2 e^{2/r}/(r-2)  (from solve). E^2 = Veff at circular.
L2 = r**2*sp.exp(2/r)/(r-2)
Veff = L2*sp.exp(-4/r)/r**2 + sp.exp(-2/r)
E2 = sp.simplify(Veff)
print("E^2(r) circular =", E2)
# specific energy E (per unit rest mass, c=1). Binding energy fraction = 1 - E
phi2 = (1+sp.sqrt(5))/2
risco = 3+sp.sqrt(5)
E2_isco = sp.simplify(E2.subs(r, risco))
print("E^2 at ISCO =", E2_isco, " => E =", sp.simplify(sp.sqrt(E2_isco)), " approx", sp.N(sp.sqrt(E2_isco)))
L2_isco = sp.simplify(L2.subs(r, risco))
print("L^2 at ISCO =", L2_isco, " => L =", sp.N(sp.sqrt(L2_isco)), "  exact:", sp.simplify(sp.sqrt(L2_isco)))
# accretion efficiency eta = 1 - E_isco  (energy radiated per unit rest mass falling from infinity to ISCO)
E_isco = sp.sqrt(E2_isco)
eta = sp.simplify(1 - E_isco)
print("eta = 1 - E_isco =", eta, " approx", sp.N(eta))
print("   GR Schwarzschild eta = 1 - sqrt(8/9) =", sp.N(1-sp.sqrt(sp.Rational(8,9))))
# Try to express E_isco in golden-ratio closed form
print("E_isco numeric:", sp.N(E_isco, 30))
# test e^{-1/r} form: at r=3+sqrt5, 2/r = 2/(3+sqrt5) = 2(3-sqrt5)/4 = (3-sqrt5)/2
arg = sp.simplify(2/risco)
print("psi at ISCO = 2/r_isco =", arg, "=", sp.nsimplify(arg,[sp.sqrt(5)]), " approx", sp.N(arg))
