# Integrating eccentric Comets with MERCURIUS (iPython)¶

MERCURIUS ia a hybrid integration scheme which combines the WHFAST and IAS15 algorithms. It smoothly transitions between the two integrators, similar to what the hybrid integrator in the MERCURY package is doing.

```
import rebound
import numpy as np
```

*semi-active*bodies, which means they can influence/be influenced by other active bodies, but are invisible to each other. This is done by setting

`testparticle_type = 1`

. Setting
`testparticle_type = 0`

would indicate that we are adding comets as
*test*bodies. * merging bodies when a collision is triggered, conserving momentum and mass. * removing particles that leave our pre-defined box. * tracking the energy lost due to ejections or collisions.

```
sim = rebound.Simulation()
np.random.seed(42)
#integrator options
sim.integrator = "mercurius"
sim.dt = 1
sim.testparticle_type = 1
#collision and boundary options
sim.collision = "direct"
sim.collision_resolve = "merge"
sim.collision_resolve_keep_sorted = 1
sim.boundary = "open"
boxsize = 200.
sim.configure_box(boxsize)
sim.track_energy_offset = 1
#simulation time
tmax = 1e4
```

Now that the preliminary setup is complete, it’s time to add some
particles to the system! When using the MERCURIUS integrator it is
important to add active bodies first and semi-active bodies later. The
`sim.N_active`

variable distinguishes massive bodies from
semi-active/test bodies.

```
#massive bodies
sim.add(m=1., r=0.005) # Sun
a_neptune = 30.05
sim.add(m=5e-5,r=2e-4,a=a_neptune,e=0.01) # Neptune
sim.N_active = sim.N
```

Now, let’s create some comets! For this simple example we are assuming that all comets have the same mass and radius.

```
# semi-active bodies
n_comets = 100
a = np.random.random(n_comets)*10 + a_neptune
e = np.random.random(n_comets)*0.009 + 0.99
inc = np.random.random(n_comets)*np.pi/2.
m = 1e-10
r = 1e-7
for i in xrange(0,n_comets):
rand = np.random.random()*2*np.pi
sim.add(m=m, r=r, a=a[i], e=e[i], inc=inc[i], Omega=0, omega=rand, f=rand)
```

We need to move to the COM frame to avoid drifting out of our simulation box. Also it is always good practice to monitor the change in energy over the course of a simulation, which requires us to calculate it before and after the simulation.

```
sim.move_to_com()
E0 = sim.calculate_energy()
```

We can visualize our setup using `rebound.OrbitPlot`

```
%matplotlib inline
fig = rebound.OrbitPlot(sim,Narc=300)
```

Alternatively, we can also use the WebGL Widget to get an interactive visualization of the simulation.

```
sim.getWidget(size=(500,300),scale=1.8*a_neptune)
```

Finally, let’s simulate our system for and check that our final relative energy error is small. The energy error is a key measure of whether the integration was performed accurately or not.

```
sim.integrate(tmax)
dE = abs((sim.calculate_energy() - E0)/E0)
print(dE)
```

```
4.604783622830939e-07
```