In this example, we study highly eccentric comets which interact with a Neptune mass planet.

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.

```
In [1]:
```import rebound
import numpy as np

First let's choose the basic properties required for the MERCURIUS integrator to run correctly. In particular, we are:

- adding comets as
*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.

```
In [2]:
```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

`sim.N_active`

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

```
In [3]:
```#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

```
In [4]:
```# 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)

```
In [5]:
```sim.move_to_com()
E0 = sim.calculate_energy()

We can visualize our setup using `rebound.OrbitPlot`

```
In [9]:
```%matplotlib inline
fig = rebound.OrbitPlot(sim,Narc=300)

```
```

```
In [10]:
```sim.getWidget(size=(500,300),scale=1.8*a_neptune)

```
```

```
In [11]:
```sim.integrate(tmax)
dE = abs((sim.calculate_energy() - E0)/E0)
print(dE)

```
```

```
In [ ]:
```