```
In [1]:
```import rebound
sim = rebound.Simulation()
sim.add(m=1., x=1., vz = 2.)

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. This is a useful set for theoretical calculations, and gives a logical behavior as the mass ratio increase, e.g., in the case of a circumbinary planet. Let's set up a binary,

```
In [2]:
```sim.add(m=1., a=1.)
sim.status()

```
```

```
In [3]:
```sim.add(m=1.e-3, a=100.)

```
In [4]:
```sim.add(primary=sim.particles[1], a=0.01)

`sim.calculate_orbits()`

. **Note that REBOUND will always output angles in the range $[-\pi,\pi]$, except the inclination which is always in $[0,\pi]$.**

```
In [5]:
```orbits = sim.calculate_orbits()
for orbit in orbits:
print(orbit)

```
```

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:

```
In [6]:
```print(sim.particles[3].calculate_orbit(sim, primary=sim.particles[1]))

```
```

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

**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,

```
In [7]:
```sim = rebound.Simulation()
sim.add(m=1.)
sim.add(a=1., e=0., inc=0.1, Omega=0.3, omega=0.1)
orbits = sim.calculate_orbits()
print(orbits[0])

```
```

```
In [8]:
```print(orbits[0].theta)

```
```

```
In [20]:
```sim = rebound.Simulation()
sim.add(m=1.)
sim.add(a=1., e=0., inc=0.1, Omega=0.3, theta = 0.4)
orbits = sim.calculate_orbits()
print(orbits[0].theta)

```
```

```
In [12]:
```import rebound
sim = rebound.Simulation()
sim.add(m=1.)
sim.add(a=1., e=0.2, Omega=0.1)
orbits = sim.calculate_orbits()
print(orbits[0])

```
```

```
In [13]:
```print(orbits[0].pomega)

```
```

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

```
In [14]:
```import rebound
sim = rebound.Simulation()
sim.add(m=1.)
sim.add(a=1., e=0.2, pomega=0.1)
orbits = sim.calculate_orbits()
print(orbits[0])

```
```

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

Finally, we can initialize particles 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,

```
In [18]:
```sim = rebound.Simulation()
sim.add(m=1.)
sim.add(a=1., e=0.1, Omega=0.3, M = 0.1)
sim.add(a=1., e=0.1, Omega=0.3, l = 0.4)
orbits = sim.calculate_orbits()
print(orbits[0].l)
print(orbits[1].l)

```
```

```
In [19]:
```import rebound
sim = rebound.Simulation()
sim.add(m=1.)
sim.add(a=1., e=0.1, omega=1.)
orbits = sim.calculate_orbits()
print(orbits[0])

```
```

**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:

```
In [5]:
```import random
import numpy as np
def simulation(par):
e,f = par
e = 10**e
f = 10**f
sim = rebound.Simulation()
sim.add(m=1.)
a = 1.
inc = random.random()*np.pi
Omega = random.random()*2*np.pi
sim.add(m=0.,a=a,e=e,inc=inc,Omega=Omega, f=f)
o=sim.calculate_orbits()[0]
if o.f < 0: # avoid wrapping issues
o.f += 2*np.pi
err = max(np.fabs(o.e-e)/e, np.fabs(o.f-f)/f)
return err
random.seed(1)
N = 100
es = np.linspace(-16.,-1.,N)
fs = np.linspace(-16.,-1.,N)
params = [(e,f) for e in es for f in fs]
pool=rebound.InterruptiblePool()
res = pool.map(simulation, params)
res = np.array(res).reshape(N,N)
res = np.nan_to_num(res)
%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib import ticker
from matplotlib.colors import LogNorm
import matplotlib
f,ax = plt.subplots(1,1,figsize=(7,5))
extent=[fs.min(), fs.max(), es.min(), es.max()]
ax.set_xlim(extent[0], extent[1])
ax.set_ylim(extent[2], extent[3])
ax.set_xlabel(r"true anomaly (f)")
ax.set_ylabel(r"eccentricity")
im = ax.imshow(res, norm=LogNorm(), vmax=1., vmin=1.e-16, aspect='auto', origin="lower", interpolation='nearest', cmap="RdYlGn_r", extent=extent)
cb = plt.colorbar(im, ax=ax)
cb.solids.set_rasterized(True)
cb.set_label("Relative Error")

```
```

```
In [6]:
```def simulation(par):
e,theta = par
e = 10**e
theta = 10**theta
sim = rebound.Simulation()
sim.add(m=1.)
a = 1.
inc = random.random()*np.pi
Omega = random.random()*2*np.pi
omega = random.random()*2*np.pi
sim.add(m=0.,a=a,e=e,inc=inc,Omega=Omega, theta=theta)
o=sim.calculate_orbits()[0]
if o.theta < 0:
o.theta += 2*np.pi
err = max(np.fabs(o.e-e)/e, np.fabs(o.theta-theta)/theta)
return err
random.seed(1)
N = 100
es = np.linspace(-16.,-1.,N)
thetas = np.linspace(-16.,-1.,N)
params = [(e,theta) for e in es for theta in thetas]
pool=rebound.InterruptiblePool()
res = pool.map(simulation, params)
res = np.array(res).reshape(N,N)
res = np.nan_to_num(res)
f,ax = plt.subplots(1,1,figsize=(7,5))
extent=[thetas.min(), thetas.max(), es.min(), es.max()]
ax.set_xlim(extent[0], extent[1])
ax.set_ylim(extent[2], extent[3])
ax.set_xlabel(r"true longitude (\theta)")
ax.set_ylabel(r"eccentricity")
im = ax.imshow(res, norm=LogNorm(), vmax=1., vmin=1.e-16, aspect='auto', origin="lower", interpolation='nearest', cmap="RdYlGn_r", extent=extent)
cb = plt.colorbar(im, ax=ax)
cb.solids.set_rasterized(True)
cb.set_label("Relative Error")

```
```

**Hyperbolic & Parabolic Orbits**

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

```
In [20]:
```sim.add(a=-0.2, e=1.4)
sim.status()

```
```

```
In [29]:
```sim = rebound.Simulation()
sim.add(m=1.)
q = 0.1
a=-1.e14
e=1.+q/np.fabs(a)
sim.add(a=a, e=e)
print(sim.calculate_orbits()[0])

```
```

**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,

```
In [19]:
```sim = rebound.Simulation()
sim.add(m=1.)
sim.add(a=1.,inc=np.pi,e=0.1, Omega=0., pomega=1.)
print(sim.calculate_orbits()[0])

```
```

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 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.

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.

```
In [ ]:
```