import sympy as sp
r,u = sp.symbols('r u', positive=True)
# Omega^2 = e^{-4/r}/(r^2 (r-1))  [G=M=c=1].  In u=1/r:
Omega2 = sp.exp(-4/r)/(r**2*(r-1))
# Orbital ang freq (coordinate). Kepler GR: Omega^2 = 1/r^3. Compare DFD small-u:
ser = sp.series(Omega2.rewrite(sp.exp), r, sp.oo, 6)
print("Omega^2 large-r (weak field) expansion:")
print(sp.simplify(ser))
# periastron / nodal precession ratio (kappa/Omega)^2 = 1-6u+4u^2
ratio2 = 1 - 6*u + 4*u**2
roots = sp.solve(sp.Eq(ratio2,0), u)
print("roots of (kappa/Omega)^2=0:", roots, [sp.N(x) for x in roots])
print("u_isco=(3-sqrt5)/4=", sp.N((3-sp.sqrt(5))/4))
# periastron precession per orbit: Delta phi_per = 2pi(Omega/kappa - 1)
# GR: (kappa/Omega)^2 = 1-6u. DFD: 1-6u+4u^2. 
# weak-field precession:
kO_GR = sp.sqrt(1-6*u)
kO_DFD = sp.sqrt(1-6*u+4*u**2)
prec_GR = sp.series(2*sp.pi*(1/kO_GR-1), u, 0, 3)
prec_DFD = sp.series(2*sp.pi*(1/kO_DFD-1), u, 0, 3)
print("GR periastron advance per orbit (in u):", prec_GR)
print("DFD periastron advance per orbit (in u):", prec_DFD)
# leading term should match 6pi u (=> 6piGM/c^2 a(1-e^2)); difference at O(u^2)
print("difference DFD-GR:", sp.simplify(sp.series(2*sp.pi*(1/kO_DFD-1) - 2*sp.pi*(1/kO_GR-1), u,0,3)))
