Raise exception when two particles are further apart than a set amount

[alelatt]

Hi, I wanted to know if there was a way of raising an exception when two particles get further apart than a set amount, similar to what is already possible by setting "exit_min_distance" for two particles that get nearer than a set amount. I noticed that "exit_max_distance" looks at the distance of a particle from the coordinate center instead of the distance between particles.

I'm interested in operating on an unbound particle (with the final objective of removing it) once it is exactly at a set distance from the rest of the system. Using an exception would help not having to set very small integration times to catch the time at which the particle gets further than the amount I set, thus saving a bit of computing time. Also, to get the time I'd have to notice that the particle was further than the set distance at the end of an integration timestep, having to go back to the start of that timestep and going forward once more with smaller timesteps, and I fear worsening the integration's precision by doing this back and forward.

After a particle becomes unbound and escapes, the system moves away from the center of mass and, sooner or later, it will also be further from the center than whatever distance was set with "exit_max_distance" and will start raising the exception, thus I can't use "exit_max_distance". I can set the distance to infinity after removing the particle, but it becomes a problem again if another particle becomes unbound and I have to start the procedure from scratch.

I'd be very happy if it was possible to add the exception as a new feature, but if it's not possible and there's a better way to do what I need to it would be very useful to know.

Also, I'm working in Python, but have a basic understanding of C. If it's better to modify the package myself and anyone can tell me how, I can try.

Thank you very much

Alessandro

[hannorein]

Hi Alessandro,

Let me know if I misunderstand your question, but here are some pointer.

If you want to find the exact time when two particles are further away than some distance, then you have an implicit problem. You'll need to do some going back and forth, using a bisection method or Newton's method.

If you are ok with just detecting the timestep after which the particles were at least some distance away, then it's very easy to do. I'd suggest you write a custom heartbeat function. The function is called every timestep so you might want to make it fast. If you can do it in C that would be most efficient (see this example on how to add a c heartbeat function from the python side). Here's a simple example that will throw an exception if two particles exceed some distance.

void heartbeat(struct reb_simulation* r){
   double distance = reb_particle_distance(&r->particles[1], &r->particles[2]);
   if (distance> 1.234){
       r->status = REB_STATUS_ESCAPE; // This will trigger an exception in python
   }
}

Note that depending on which integrator you use and the precise settings, there might be a need to use safe_mode=1 to get a precise answer.

You're right about the centre of mass slowing drifting after removing a particle. I recommend you call sim.move_to_com() after removing a particle to make sure the COM stays close to the origin.

[alelatt]

Thank you very much for the detailed answer! This was exactly what I needed.

[hannorein]

Glad to hear!

[alelatt]

Sorry to bother again, is there an easy way to pass the maximum distance to the heartbeat function from python or do I have to encode it in the C file as you did in your example? Ideally I'd like to be able to update that distance as the code progresses, in a similar way as I can do with "exit_min_distance" and "exit_max_distance"

[hannorein]

Yes. Have a look at the example I've linked above. In short, define a global variable in c with

double distance =0;
void heartbeat(...) {
   ...
}

and then set it from python using something like:

c_double.in_dll(clibheartbeat,"distance").value = 1.234

[alelatt]

Great, thank you! I thought that I could only read the variable from Python. Thanks again