While you can always check particles' states after every call to sim.integrate, you might want to check for particular conditions every timestep that you would otherwise miss. You could make the array of times to which you integrate increasingly finer, but this will slow things down, as will writing a Python function that gets called every timestep.
This notebook shows how to keep track of a particle's minimum distance from the central star using REBOUNDx to make sure everything is done quickly in C. We first create a simulation:
In [1]:
import rebound
import reboundx
import numpy as np
sim = rebound.Simulation()
sim.add(m=1., hash="Sun")
sim.add(a=1., e=0.5, f=np.pi)
rebx= reboundx.Extras(sim)
ps = sim.particles
Now we add our effect, and for each particle whose distance we want to track, we have to add a min_distance parameter. This parameter will automatically get updated to lower values any timestep where the particle's distance drops.
In [2]:
track = rebx.load_operator("track_min_distance")
rebx.add_operator(track)
ps[1].params["min_distance"] = 5.
By default, the effect calculates the distance from sim.particles[0]. You can specify a different particle through its hash (added as a min_distance_from parameter to the particle we want to track). So the following line is optional in our case:
In [3]:
ps[1].params["min_distance_from"] = ps[0].hash
We might also be interested in what the orbit looks like at these special times. We can add a min_distance_orbit parameter to our particle to store the instantaneous (heliocentric) orbit at the time corresponding to our min_distance. We have to set it to a rebound.Orbit instance, so we make one (defaults to all zeros, but will get updated).
Note that when we add objects as parameters, we're responsible for that memory, so we need to keep a reference to that object (like we do here with 'o'). If we did ps[1].params["min_distance_orbit"] == rebound.Orbit(), or did all of this within a function that later returned, python could think we're done with the rebound.Orbit() object and garbage collect it, causing undefined behavior.
In [4]:
o = rebound.Orbit()
ps[1].params["min_distance_orbit"] = o
We now integrate as usual, and see that the parameter has been updated:
In [5]:
sim.integrate(10.)
ps[1].params["min_distance"]
Out[5]:
We can check the corresponding orbit. In this two-body case most of the elements haven't changed, but we see our timestep caught it just past pericenter:
In [6]:
print(ps[1].params["min_distance_orbit"])
Let's try artificially moving our particle inward:
In [7]:
ps[1].x /= 3.
ps[1].y /= 3.
ps[1].z /= 3.
Now if we keep integrating further, both min_distance and min_distance_orbit will get updated:
In [8]:
sim.integrate(20.)
print("min_distance = {0}".format(ps[1].params["min_distance"]))
print("min_distance_orbit = {0}".format(ps[1].params["min_distance_orbit"]))
In [ ]: