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()
So this worked as expected. We went down to 2 particles, and 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]);
Uh oh. The problem here is that we kept updating 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!