import sympy as sp
import mpmath as mp
mp.mp.dps=40
r=sp.symbols('r',positive=True)
risco=3+sp.sqrt(5)
E2=(r-1)*sp.exp(-2/r)/(r-2)
L2=r**2*sp.exp(2/r)/(r-2)
E_isco=sp.sqrt(E2.subs(r,risco))
eta=1-E_isco
print("E_isco =", float(E_isco), " eta(DFD) =", float(eta*100),"%  eta(GR)=",float((1-mp.sqrt(mp.mpf(8)/9))*100),"%")
print("L_isco(DFD)=",float(sp.sqrt(L2.subs(r,risco)))," L_isco(GR)=sqrt(12)=",float(mp.sqrt(12)))
# E2(r)=(r-1)e^{-2/r}/(r-2): verify exact circular energy form by re-derivation residual
psi=2/r
Veff=sp.exp(-2*psi)/r**2 * (r**2*sp.exp(2/r)/(r-2)) + sp.exp(-psi)
print("E2 residual:", sp.simplify(Veff-(r-1)*sp.exp(-2/r)/(r-2)))
# QPO ratio exact recompute residual
u=sp.symbols('u',positive=True)
print("(kappa/Omega)^2 exact = 1-6u+4u^2 ; at u_isco:",float((1-6*u+4*u**2).subs(u,(3-sp.sqrt(5))/4)))
