Removing particles from the simulation

Let's set up a simple simulation with 10 bodies, and have the simulation assign them unique hashes so we can keep track of them (see UniquelyIdentifyingParticles.ipynb).


In [1]:
import rebound
import numpy as np

sim = rebound.Simulation()
sim.add(m=1.)
sim.particles[0].hash = sim.generate_unique_hash()
for i in range(1,10):
    sim.add(a=i)
    sim.particles[i].hash = sim.generate_unique_hash()
sim.move_to_com()

print("Particle hashes:{0}".format([sim.particles[i].hash for i in range(sim.N)]))


Particle hashes:[10904, 10905, 10906, 10907, 10908, 10909, 10910, 10911, 10912, 10913]

We could also add custom names to particles:


In [2]:
sim.add(a=10, hash="Saturn")
print("Particle hashes:{0}".format([sim.particles[i].hash for i in range(sim.N)]))


Particle hashes:[10904, 10905, 10906, 10907, 10908, 10909, 10910, 10911, 10912, 10913, 4066125545]

Now let's do a simple example where we do a short initial integration to isolate the particles that interest us for a longer simulation:


In [3]:
Noutputs = 1000
xs = np.zeros((sim.N, Noutputs))
ys = np.zeros((sim.N, Noutputs))
times = np.linspace(0.,50*2.*np.pi, Noutputs, endpoint=False)
for i, time in enumerate(times):
    sim.integrate(time)
    xs[:,i] = [sim.particles[j].x for j in range(sim.N)]
    ys[:,i] = [sim.particles[j].y for j in range(sim.N)]
%matplotlib inline
import matplotlib.pyplot as plt
fig,ax = plt.subplots(figsize=(15,5))
for i in range(sim.N):
    plt.plot(xs[i,:], ys[i,:])
ax.set_aspect('equal')


At this stage, we might be interested in particles that remained within some semimajor axis range, particles that were in resonance with a particular planet, etc. Let's imagine a simple (albeit arbitrary) case where we only want to keep particles that had $x < 0$ at the end of the preliminary integration. Let's first print out the particle hashes and x positions.


In [4]:
print("Hash\tx")
for i in range(sim.N):
    print("{0}\t{1}".format(sim.particles[i].hash, xs[i,-1]))


Hash	x
10904	0.0
10905	0.9510565162930091
10906	-1.0717399536588612
10907	-2.2765351809117464
10908	0.15703926303973234
10909	-4.897155109586999
10910	-4.824394540939856
10911	-2.2862837234997975
10912	2.111033731282993
10913	5.290067270630363
4066125545	-8.776421396714463

4066125545 is the hash corresponding to the string "Saturn" we added above. Next, let's use the remove() function to filter out particlse. As an argument, we pass the corresponding index in the particles array.


In [5]:
for i in reversed(range(1,sim.N)):
    if xs[i,-1] > 0:
        sim.remove(i)
print("Number of particles after cut = {0}".format(sim.N))
print("Hashes of remaining particles = {0}".format([p.hash for p in sim.particles]))


Number of particles after cut = 7
Hashes of remaining particles = [10904, 10906, 10907, 10909, 10910, 10911, 4066125545]

By default, the remove() function removes the i-th particle from the particles array, and shifts all particles with higher indices down by 1. This ensures that the original order in the particles array is preserved (e.g., to help with output).

By running through the planets in reverse order above, we are guaranteed that when a particle with index i gets removed, the particle replacing it doesn't need to also be removed (we already checked it).

If you have many particles and many removals (or you don't care about the ordering), you can save the reshuffling of all particles with higher indices with the flag keepSorted=0:


In [6]:
sim.remove(2, keepSorted=0)
print("Number of particles after cut = {0}".format(sim.N))
print("Hashes of remaining particles = {0}".format([p.hash for p in sim.particles]))


Number of particles after cut = 6
Hashes of remaining particles = [10904, 10906, 4066125545, 10909, 10910, 10911]

We see that the particles array is no longer sorted by hash. Note that the default keepSorted=1 only keeps things sorted (i.e., if they were sorted by hash to start with). If you custom-assign hashes with strings (or manually) the default will simply preserve the original order.

Because in general particles can change positions in the particles array, a more robust way of referring to them (rather than with their index) is through their hash, which won't change. You can pass sim.remove either the hash directly, or if you pass a string, it will be automatically converted to its corresponding hash:


In [7]:
sim.remove(hash="Saturn")
print("Number of particles after cut = {0}".format(sim.N))
print("Hashes of remaining particles = {0}".format([p.hash for p in sim.particles]))


Number of particles after cut = 5
Hashes of remaining particles = [10904, 10906, 10909, 10910, 10911]

In [ ]: