I have here a block of code that finds when two particles have a close encounter with each other. What about 5 particles?

Question

The reason for this code is to find the time it takes for a planetary system with 2 planets(particles) to have a close encounter between each other. I, however, would like to solve that time for a planetary system with 5 particles. This will be printed at the last line print(sim.t).

http://rebound.readthedocs.io/en/latest/examples.html Here is some documentation on the usage of REBOUND, where you can find the best explanation on some of the declarations I've made.

sim = setupSimulation()
sim.exit_min_distance = 0.01
Noutputs = 10000
year = 2.*np.pi
times = np.linspace(0.,10E+9.*year, Noutputs)
distances = np.zeros(Noutputs)
ps = sim.particles
try:
    for i,time in enumerate(times):
        sim.integrate(time)
        dp = ps[1] - ps[2]   
        distances[i] = np.sqrt(dp.x*dp.x+dp.y*dp.y+dp.z*dp.z)
except rebound.Encounter as error:
    print(error)

print(sim.t)

Line 2 defines a close encounter between the two particles (simulation will stop when this value is reached).

Line 11 takes the difference in coordinates of particles, if <= .01, prints that there is a close encounter.

I was thinking along the line of an if-then statement:

if ps[1] - ps[2] <= .01:
    dp = ps[1] - ps[2]
else ps[2] - ps[3] <= .01:
    dp = ps[2] - ps[3]

and so on...

I'd like to make sure this works before I run it simply because the simulation time is 10e+9 years, and can only be ran on a local supercomputer to obtain results in a reasonable amount of time.

You sure that shouldn't be `ps[0]`? What is `sim`? What type is `p`? — Eric, May 27 '17 at 00:32
`ps[0]` is the star that particles 1-5 are orbiting. I am using an N body integrator called rebound, where `sim` is the simulation itself. Would you like to see the code that was left out? Prior to this block I use `sim.add()` six times to add one star and 5 rotating bodies. I'm not entirely sure what you mean by your last question. — dlsj, May 27 '17 at 00:49
We probably can't answer this question without a lot more information than you've given us. To start with, you *really* need to describe the simulation library you're using (a link to its documentation would be a good start). Without knowing what `sim` does, we can't possibly know what's going on in the code you've shown. It's also not really clear what you're asking for. What output do you expect when you're running your simulator with 5 particles (and how does that differ from when there are only two)? — Blckknght, May 27 '17 at 01:12
@Blckknght Hey, so I actually hesitated to add more information since I was once told that I shouldn't add too much. I've edited the original post, leading to a link of a few examples REBOUND can be used. If I were to include the rest of my code, you'd be able to see exactly what the end goal was, but again I hesitated. The whole point of this is to find when my planetary system will have a close encounter, and return to me the exact time it will happen. — dlsj, May 27 '17 at 06:00

whp · Accepted Answer · 2017-05-27T01:40:10.520

I don't see why the if-statement approach wouldn't work for detecting if any two particles get close enough to each other. Edit: writing this example as if your code sample for 2 particles is working.

Depending on if you ever need to check non-sequential particles, you'll need to generate each combination and check it (shown in example).

It looks like you'll also need a different enumeration strategy for distances. If the indexes don't have to match the time counter, you could increment a separate counter (shown in example).

You could also use itertools to generate the if statements to avoid typos and increase clarity:

import itertools

# other setup code

counter = 0
try:
    for time in times:
        sim.integrate(time)
        for p0, p1 in itertools.combinations(range(1, 6)):  # generates combinations of the particles numbered 1 thru 5
            dp = ps[p0] - ps[p1]
            distances[counter] = np.sqrt(dp.x*dp.x+dp.y*dp.y+dp.z*dp.z)
            counter += 1
except rebound.Encounter as error:
    print(error)

I'll add a little more - to make sure your code works, you probably should write some kind of test suite for it. What happens if it resolves vs. never resolves.

You could create an object to mock sim and its methods, returning bodies in positions that you know will or will not meet your exit conditions without all the overhead of needing a supercomputer. You could just hard-code them in a small list.

Hey, thanks for the example, I should be working on non-sequential particles once I spit out the data for my current project, so thanks for the future help! However, will I not need to consider time on just a simple hard code? — dlsj, May 27 '17 at 06:05
Depends what you're testing. For integration tests, you could also set up initial conditions that you know should exit quickly (e.g. co-orbital particles moving in opposite directions, coplanar orbits, etc.) so that you can run locally. — whp, May 28 '17 at 02:21
How would you go about locating which two planets, out of the 5, had the close encounter? — dlsj, Jun 20 '17 at 02:11

I have here a block of code that finds when two particles have a close encounter with each other. What about 5 particles?

1 Answers1