gr_unphysical_flip.py

#!/usr/bin/env python



import rebound 
import reboundx
import numpy as np
import astropy.constants as const

##Constants
G=const.G.cgs.value
M_sun=const.M_sun.cgs.value
c=const.c.cgs.value
pc=const.pc.cgs.value
##Converting speed of light to code units. Units are pc, solar masses, and G=1
cc=c/(G*M_sun/pc)**0.5

##Initializing simulation. 
end_r=np.array([ 2.43429405e-05, -6.34018455e-05, -2.87889574e-04])
end_v=np.array([-43247.28772886,  42257.14066544, 152805.55464892])
sim=rebound.Simulation()
sim.add(m=4e6)
sim.add(m=100.05, x=end_r[0], y=end_r[1], z=end_r[2], vx=end_v[0], vy=end_v[1], vz=end_v[2])
sim.move_to_com()
rebx = reboundx.Extras(sim)
gr=rebx.add("gr")
gr.params["c"]=cc

##Various quantities of interest. Rg is the gravitational radius of the central black hole. Rp is the pericenter
rg=sim.particles[0].m/cc**2.
print("Rg: {0}".format(rg))
print("Initial radius/Rg: {0}".format(np.linalg.norm(end_r/rg)))
orbs=sim.calculate_orbits(primary=sim.particles[0])
aa=orbs[0].a
rp=(1.-orbs[0].e)*aa
print("Initial sma:",aa)
print("Rp/Rg:",rp/rg)


##Integrate simulation
p_in=2.0*np.pi*(aa**3.0/sim.particles[0].m)**0.5
sim.integrate(p_in)
##Final semimajor axis: negative!
orbs=sim.calculate_orbits(primary=sim.particles[0])
print("Final sma:",orbs[0].a)