For convenience, REBOUND offers simple functionality for converting units. One implicitly sets the units for the simulation through the values used for the initial conditions, but one has to set the appropriate value for the gravitational constant G, and sometimes it is convenient to get the output in different units.
The default value for G is 1, so one can:
a) use units for the initial conditions where G=1 (e.g., AU, $M_\odot$, yr/$2\pi$)
b) set G manually to the value appropriate for the adopted initial conditions, e.g., to use SI units,
In [1]:
import rebound
import math
sim = rebound.Simulation()
sim.G = 6.674e-11
c) set rebound.units:
In [2]:
sim.units = ('yr', 'AU', 'Msun')
print("G = {0}.".format(sim.G))
When you set the units, REBOUND converts G to the appropriate value for the units passed (must pass exactly 3 units for mass length and time, but they can be in any order). Note that if you are interested in high precision, you have to be quite particular about the exact units.
As an aside, the reason why G differs from $4\pi^2 \approx 39.47841760435743$ is mostly that we follow the convention of defining a "year" as 365.25 days (a Julian year), whereas the Earth's sidereal orbital period is closer to 365.256 days (and at even finer level, Venus and Mercury modify the orbital period). G would only equal $4\pi^2$ in units where a "year" was exactly equal to one orbital period at $1 AU$ around a $1 M_\odot$ star.
Adding particles
If you use sim.units at all, you need to set the units before adding any particles. You can then add particles in any of the ways described in WHFast.ipynb. You can also add particles drawing from the horizons database (see Churyumov-Gerasimenko.ipynb). If you don't set the units ahead of time, HORIZONS will return initial conditions in units of AU, $M_\odot$ and yrs/$2\pi$, such that G=1.
Above we switched to units of AU, $M_\odot$ and yrs, so when we add Earth:
In [3]:
sim.add('Earth')
ps = sim.particles
import math
print("v = {0}".format(math.sqrt(ps[0].vx**2 + ps[0].vy**2 + ps[0].vz**2)))
we see that the velocity is correctly set to approximately $2\pi$ AU/yr.
If you'd like to enter the initial conditions in one set of units, and then use a different set for the simulation, you can use the sim.convert_particle_units function, which converts both the initial conditions and G. Since we added Earth above, we restart with a new Simulation instance; otherwise we'll get an error saying that we can't set the units with particles already loaded:
In [4]:
sim = rebound.Simulation()
sim.units = ('m', 's', 'kg')
sim.add(m=1.99e30)
sim.add(m=5.97e24,a=1.5e11)
sim.convert_particle_units('AU', 'yr', 'Msun')
sim.status()
We first set the units to SI, added (approximate values for) the Sun and Earth in these units, and switched to AU, yr, $M_\odot$. You can see that the particle states were converted correctly--the Sun has a mass of about 1, and the Earth has a distance of about 1.
Note that when you pass orbital elements to sim.add, you must make sure G is set correctly ahead of time (through either 3 of the methods above), since it will use the value of sim.G to generate the velocities:
In [5]:
sim = rebound.Simulation()
print("G = {0}".format(sim.G))
sim.add(m=1.99e30)
sim.add(m=5.97e24,a=1.5e11)
sim.status()
The orbital speed of Earth is $\sim 3\times 10^4$ m/s, but since we didn't correctly set G ahead of time, we get $\sim 3\times 10^9$ m/s, so the Earth would fly off the Sun in this simulation.