Adding Post-Newtonian general relativity corrections

It's easy to add post-newtonian corrections to your REBOUND simulations with REBOUNDx. Let's start with a simulation without GR:



In [1]:

    
import rebound
sim = rebound.Simulation()
sim.add(m=1., hash="star") # Sun
sim.add(m=1.66013e-07,a=0.387098,e=0.205630, hash="planet") # Mercury-like
sim.move_to_com() # Moves to the center of momentum frame
ps = sim.particles

sim.integrate(10.)
print("pomega = %.16f"%sim.particles[1].pomega)









    



pomega = 0.0000000000000000

As expected, the pericenter did not move at all. Now let's add GR



In [2]:

    
import reboundx
rebx = reboundx.Extras(sim)
gr = rebx.load_force("gr")
rebx.add_force(gr)

The GR effects need you to set the speed of light in the right units. The constants module has a set of constants in REBOUND's default units of AU, solar masses and yr/$2\pi$ (such that G=1). If you want to use other units, you'd need to calculate c.



In [3]:

    
from reboundx import constants
gr.params["c"] = constants.C

By default, the gr and gr_potential effects assume that the massive particle is at index 0 (with gr_full all particles are "sources" so this is not an issue). If the massive particle has a different index, or you think it might move from index 0 in the particles array (e.g. due to a custom merger routine), you can attach a gr_source flag to it to identify it as the massive particle with:



In [4]:

    
ps["star"].params["gr_source"] = 1

Now we integrate as normal. We monitor the total Hamiltonian. Unlike other forces where we can calculate a separate potential, here and with gr_full the forces are velocity dependent, which means the momentum is not just mv in a Hamiltonian framework. So rather than using sim.calculate_energy and adding a potential, gr_hamiltonian calculates the full thing (classical Hamiltonian + gr).



In [5]:

    
deltat = 100.
E0 = rebx.gr_hamiltonian(gr)
sim.integrate(sim.t + deltat)
Ef = rebx.gr_hamiltonian(gr)

print("pomega = %.16f"%sim.particles[1].pomega)
juliancentury = 628.33195 # in yr/2pi
arcsec = 4.8481368e-06 # in rad
print("Rate of change of pomega = %.4f [arcsec / Julian century]"% (sim.particles[1].pomega/deltat*juliancentury/arcsec))
print("Relative error on the relativistic Hamiltonian = {0}".format(abs(Ef-E0)/abs(E0)))









    



pomega = 0.0000332591781488
Rate of change of pomega = 43.1048 [arcsec / Julian century]
Relative error on the relativistic Hamiltonian = 1.3578502228769841e-15

As expected there was pericenter precession. The literature value is 42.98 arcsec / century.

Variants

Above we added the gr effect, but there are two other implementations gr_potential and gr_full. Before running any serious simulations, you should read the more detailed descriptions at http://reboundx.readthedocs.io/en/latest/effects.html to see which implementation is appropriate for your application.