Sometimes we are not interested in particles that get too far from the central body. Here we will define a radius beyond which we remove particles from the simulation. Let's set up an artificial situation with 2 planets, and the inner one simply moves radially outward with $v > v_{escape}$.

```
In [1]:
```import rebound
import numpy as np
def setupSimulation():
sim = rebound.Simulation()
sim.integrator = "ias15" # IAS15 is the default integrator, so we don't need this line
sim.add(m=1., id=0)
sim.add(m=1e-3,x=1.,vx=2., id=1)
sim.add(m=1e-3,a=1.25,M=np.pi/2, id=2)
sim.move_to_com()
return sim

We have assigned each particle an ID for later reference (see IDs.ipynb for more information).

```
In [2]:
```sim = setupSimulation()
sim.status()

```
```

Now let's run a simulation for 20 years (in default units where $G=1$, and thus AU, yr/2$\pi$, and $M_\odot$, see Units.ipynb for how to change units), and set up a 50 AU sphere beyond which we remove particles from the simulation. We can do this by setting the `exit_max_distance`

flag of the simulation object. If a particle's distance (from the origin of whatever inertial reference frame chosen) exceeds `sim.exit_max_distance`

, an exception is thrown.

If we simply call `sim.integrate()`

, the program will crash due to the unhandled exception when the particle escapes, so we'll create a `try`

-`except`

block to catch the exception.

```
In [3]:
```sim = setupSimulation() # Resets everything
sim.exit_max_distance = 50.
Noutputs = 1000
times = np.linspace(0,20.*2.*np.pi,Noutputs)
xs = np.zeros((3,Noutputs))
ys = np.zeros((3,Noutputs))
for i,time in enumerate(times):
try:
sim.integrate(time)
except rebound.Escape as error:
print(error)
max_d2 = 0.
for p in sim.particles:
d2 = p.x*p.x + p.y*p.y + p.z*p.z
if d2>max_d2:
max_d2 = d2
mid = p.id
sim.remove(id=mid)
for j in range(sim.N):
xs[j,i] = sim.particles[j].x
ys[j,i] = sim.particles[j].y

```
```

Let's check that the particle 1 was correctly removed from the simulation:

```
In [4]:
```sim.status()

```
```

`particles[1]`

(which had `id = 1`

before) has evidently been replaced with `particles[2]`

). By default, `remove()`

preserves the ordering in the `particles`

array (see IDs.ipynb for more info). Now let's plot what we got:

```
In [5]:
```%matplotlib inline
import matplotlib.pyplot as plt
fig,ax = plt.subplots(figsize=(15,5))
for i in range(3):
ax.plot(xs[i,:], ys[i,:])
ax.set_aspect('equal')
ax.set_xlim([-2,10]);

```
```

`xs[1]`

with `particles[1].x`

after `particle[1]`

was removed. This means that following the removal, `xs[1]`

all of a sudden started getting populated by the values for the particle with `ID=2`

. This is why the radial green trajectory (horizontal line along the $x$ axis) all of a sudden jumps onto the roughly circular orbit corresponding to the outer particle with `ID=2`

(originally red). One way to fix these problems is:

```
In [6]:
```sim = setupSimulation() # Resets everything
sim.exit_max_distance = 50.
Noutputs = 1000
times = np.linspace(0,20.*2.*np.pi,Noutputs)
xs = np.zeros((3,Noutputs))
ys = np.zeros((3,Noutputs))
for i,time in enumerate(times):
try:
sim.integrate(time)
except rebound.Escape as error:
print(error)
max_d2 = 0.
for p in sim.particles:
d2 = p.x*p.x + p.y*p.y + p.z*p.z
if d2>max_d2:
max_d2 = d2
mid = p.id
sim.remove(id=mid)
for j in range(sim.N):
xs[sim.particles[j].id,i] = sim.particles[j].x
ys[sim.particles[j].id,i] = sim.particles[j].y
fig,ax = plt.subplots(figsize=(15,5))
for i in range(3):
ax.plot(xs[i,:], ys[i,:])
ax.set_aspect('equal')
ax.set_xlim([-2,10]);

```
```

Much better! Since at the beginning of the integration the IDs match up with the corresponding indices in the `xs`

and `ys`

arrays, we solved problem by using the IDs as indices throughout the simulation.

As an aside, the horizontal drift of the circular orbit is a real effect: in the center of mass frame, if the Jupiter-mass planet is drifting right at some speed, the Sun must be moving at a speed lower by a factor of approximately 1000 (their mass ratio) in the opposite direction, so the Sun-particle2 system slowly drifts left. If we integrated long enough, this would mean all our particles would eventually leave our box.

If we wanted to make sure things stayed in the box, we could additionally move to new center of mass frame after each removal of a particle, but this would introduce unphysical jumps in the remaining particles' time series, since their coordinates are measured between different inertial frames. Of course, whether this matters depends on the application!