Orbital Elements

Note: All angles for orbital elements are in radians

We can add particles to a simulation by specifying cartesian components:

Any components not passed automatically default to 0. REBOUND can also accept orbital elements.

Reference bodies

As a reminder, there is a one-to-one mapping between (x,y,z,vx,vy,vz) and orbital elements, and one should always specify what the orbital elements are referenced against (e.g., the central star, the system's barycenter, etc.). The differences betwen orbital elements referenced to these centers differ by $\sim$ the mass ratio of the largest body to the central mass. By default, REBOUND always uses Jacobi elements, which for each particle are always referenced to the center of mass of all particles with lower index in the simulation.

For the painstaking user: When separating out the center of mass degree of freedom and reducing the N body problem to N-1 Kepler problems and interaction terms, there are a number of possible Hamiltonian splittings (see e.g., Hernandez & Dehnen 2017), and different possible choices for the primary mass in each of the separate Kepler problems. REBOUND takes this primary mass to be the total mass of all the particles in the simulation. If particles are added from the inside out, this gives logical behavior in the limit of a hierarchical system, even for large masses (one can think of it as setting up our new particle in a 2-body orbit around all the interior mass concentrated at the interior particles' center of mass).

Let's set up a stellar binary,

We always have to pass a semimajor axis (to set a length scale), but any other elements are by default set to 0. Notice that our second star has the same vz as the first one due to the default Jacobi elements. Now we could add a distant planet on a circular orbit,

This planet is set up relative to the binary center of mass (again due to the Jacobi coordinates), which is probably what we want. But imagine we now want to place a test mass in a tight orbit around the second star. If we passed things as above, the orbital elements would be referenced to the binary/outer-planet center of mass. We can override the default by explicitly passing a primary (any instance of the Particle class):

All simulations are performed in Cartesian elements, so to avoid the overhead, REBOUND does not update particles' orbital elements as the simulation progresses. However, you can always access any orbital element through, e.g., sim.particles[1].inc (see the diagram, and table of orbital elements under the Orbit structure at http://rebound.readthedocs.org/en/latest/python_api.html). This will calculate that orbital element individually--you can calculate all the particles' orbital elements at once with sim.calculate_orbits(). REBOUND will always output angles in the range $[-\pi,\pi]$, except the inclination which is always in $[0,\pi]$.

Notice that there is always one less orbit than there are particles, since orbits are only defined between pairs of particles. We see that we got the first two orbits right, but the last one is way off. The reason is that again the REBOUND default is that we always get Jacobi elements. But we initialized the last particle relative to the second star, rather than the center of mass of all the previous particles.

To get orbital elements relative to a specific body, you can manually use the calculate_orbit method of the Particle class:

though we could have simply avoided this problem by adding bodies from the inside out (second star, test mass, first star, circumbinary planet).

When you access orbital elements individually, e.g., sim.particles[1].inc, you always get Jacobi elements. If you need to specify the primary, you have to do it with sim.calculate_orbit() as above.

Edge cases and orbital element sets

Different orbital elements lose meaning in various limits, e.g., a planar orbit and a circular orbit. REBOUND therefore allows initialization with several different types of variables that are appropriate in different cases. It's important to keep in mind that the procedure to initialize particles from orbital elements is not exactly invertible, so one can expect discrepant results for elements that become ill defined. For example,

The problem here is that $\omega$ (the angle from the ascending node to pericenter) is ill-defined for a circular orbit, so it's not clear what we mean when we pass it, and we get spurious results for both $\omega$ and $f$, since the latter is also undefined as the angle from pericenter to the particle's position. However, the true longitude $\theta$, the broken angle from the $x$ axis to the ascending node = $\Omega + \omega + f$, and then to the particle's position, is always well defined:

To be clearer and ensure we get the results we expect, we could instead pass theta to specify the longitude we want, e.g.

Here we have a planar orbit, in which case the line of nodes becomes ill defined, so $\Omega$ is not a good variable, but we pass it anyway! In this case, $\omega$ is also undefined since it is referenced to the ascending node. Here we get that now these two ill-defined variables get flipped. The appropriate variable is pomega ($\varpi = \Omega + \omega$), which is the angle from the $x$ axis to pericenter:

We can specify the pericenter of the orbit with either $\omega$ or $\varpi$:

Note that if the inclination is exactly zero, REBOUND sets $\Omega$ (which is undefined) to 0, so $\omega = \varpi$.

Finally, we can specify the position of the particle along its orbit using mean (rather than true) longitudes or anomalies (for example, this might be useful for resonances). We can either use the mean anomaly $M$, which is referenced to pericenter (again ill-defined for circular orbits), or its better-defined counterpart the mean longitude l $= \lambda = \Omega + \omega + M$, which is analogous to $\theta$ above,

REBOUND calculates the mean longitude in such a way that it smoothly approaches $\theta$ in the limit of $e\rightarrow0$:

In summary, you can specify the phase of the orbit through any one of the angles M, f, theta or l=$\lambda$. Additionally, one can instead use the time of pericenter passage T. This time should be set in the appropriate time units, and you'd initialize sim.t to the appropriate time you want to start the simulation.

Accuracy

As a test of accuracy and demonstration of issues related to the last section, let's test the numerical stability by intializing particles with small eccentricities and true anomalies, computing their orbital elements back, and comparing the relative error. We choose the inclination and node longitude randomly:

We see that the behavior is poor, which is physically due to $f$ becoming poorly defined at low $e$. If instead we initialize the orbits with the true longitude $\theta$ as discussed above, we get much better results:

Hyperbolic & Parabolic Orbits

REBOUND can also handle hyperbolic orbits, which have negative $a$ and $e>1$:

Currently there is no support for exactly parabolic orbits, but we can get a close approximation by passing a nearby hyperbolic orbit where we can specify the pericenter = $|a|(e-1)$ with $a$ and $e$. For example, for a 0.1 AU pericenter,

Retrograde Orbits

Orbital elements can be counterintuitive for retrograde orbits, but REBOUND tries to sort them out consistently. This can lead to some initially surprising results. For example,

We passed $\Omega=0$ and $\varpi=1.$. For prograde orbits, $\varpi = \Omega + \omega$, so we'd expect $\omega = 1$, but instead we get $\omega=-1$. If we think about things physically, $\varpi$ is the angle from the $x$ axis to pericenter, measured in the positive direction (counterclockwise) defined by $z$. $\Omega$ is always measured in this same sense, but $\omega$ is always measured in the orbital plane in the direction of the orbit. For retrograde orbits, this means that $\omega$ is measured in the opposite sense to $\Omega$, so $\varpi = \Omega - \omega$, which is why we got $\omega = -1$.

Similarly, the retrograde version of $\theta = \Omega + \omega + f$ is $\theta = \Omega - \omega - f$, and l = $\lambda = \Omega + \omega + M$ becomes $\lambda = \Omega - \omega - M$. REBOUND chooses these conventions based on whether $i < \pi/2$, which means that if you were tracking $\varpi$ for nearly polar orbits, you would get unphysical jumps if the orbits crossed back and forth between prograde and retrograde. Of course, $\varpi$ is not a good angle at such high inclinations, and only has physical meaning when the orbital plane nearly coincides with the reference plane.

Exceptions

Adding a particle or getting orbital elements from particles in a simulation should never yield NaNs in any of the structure fields. Please let us know if you find a case that does.

In cases where it would return a NaN, REBOUND will raise a ValueError. The only cases that should do so when adding a particle are 1) passing an eccentricity of exactly 1. 2) passing a negative eccentricity. 3) Passing $e>1$ if $a>0$. 4) Passing $e<1$ if $a<0$. 5) Passing a longitude or anomaly for a hyperbolic orbit that's beyond the range allowed by the asymptotes defined by the hyperbola. You will also get errors if you try to initialize particles with orbital elements manually with rebound.Particle().

When obtaining orbital elements from a Particle structure, REBOUND will raise a ValueError if 1) the primary's mass is zero, or 2) the particle's and primary's position are the same.

Negative inclinations

While inclinations are only defined in the range $[0,\pi]$, you can also pass negative inclinations when adding particles in REBOUND. This is interpreted as referencing $\Omega$ and $\omega$ to the descending, rather than the ascending node. So for example, if one set up particles with the same $\Omega$ and a range of inclinations distributed around zero, one would obtain what one might expect, i.e. a set of orbits that are all rotated around the same line of nodes.

Jacobi masses

There is a classical Hamiltonian splitting for the N-body problem (see e.g., Wisdom & Holman 1991) that when expanded to first order in the planet/star mass ratio, gives an interaction Hamiltonian with the same form as the disturbing function for an exterior perturber. This makes it particularly attractive for analytic or semi-analytic studies. In this splitting, the masses of the primaries for each planet take on a particular form. One can add particles using these jacobi masses with the jacobi_masses flag. By default this flag is false and the primary mass is the total mass of all particles in the simulation (see the top of this ipynb).

The jacobi mass and default mass assigned by REBOUND always agree for the first particle, but differ for all the rest

We can calculate orbital elements using jacobi masses by using the same flag in sim.calculate_orbits