REBOUND is a gravitational N-body integrator. But you can also use it to integrate systems with additional, non-gravitational forces.
This tutorial gives you a very quick overview of how that works. Implementing additional forces in python as below will typically be a factor of a few slower than a C implementation. For a library that has C implementations for several commonly used additional effects (with everything callable from Python), see REBOUNDx.
We'll start by adding two particles, the Sun and an Earth-like planet, to REBOUND.
import rebound sim = rebound.Simulation() sim.integrator = "whfast" sim.add(m=1.) sim.add(m=1e-6,a=1.) sim.move_to_com() # Moves to the center of momentum frame
We could integrate this system, and the planet would go around the star at a fixed orbit with $a=1$ forever. Let's add an additional constant force that acts on the planet and points in one direction $F_x = m\cdot c$, where $m$ is the planet's mass and $c$ a constant. This is called the Stark problem. In python, we can describe this with the following function
ps = sim.particles c = 0.01 def starkForce(reb_sim): ps.ax += c
Next, we need to tell REBOUND about this function.
sim.additional_forces = starkForce
Now we can just integrate as usual. Let's keep track of the eccentricity as we integrate as it will change due to the additional force.
import numpy as np Nout = 1000 es = np.zeros(Nout) times = np.linspace(0.,100.*2.*np.pi,Nout) for i, time in enumerate(times): sim.integrate(time, exact_finish_time=0) # integrate to the nearest timestep so WHFast's timestep stays constant es[i] = sim.particles.e
And let's plot the result.
%matplotlib inline import matplotlib.pyplot as plt fig = plt.figure(figsize=(15,5)) ax = plt.subplot(111) plt.plot(times, es);
You can see that the eccentricity is oscillating between 0 and almost 1.
Note that the function
starkForce(reb_sim) above receives the argument
reb_sim when it is called. This is a pointer to the simulation structure. Instead of using the global
ps variable to access particle data, one could also use
reb_sim.contents.particles. This could be useful when one is running multiple simulations in parallel or when the particles get added and removed (in those cases
particles might change). The
contents has the same meaning as a
-> in C, i.e. follow the memory address. To find out more about pointers, check out the ctypes documentation.
The previous example assumed a conservative force, i.e. we could describe it as a potential as it is velocity independent. Now, let's assume we have a velocity dependent force. This could be a migration force in a protoplanetary disk or PR drag. We'll start from scratch and add the same two particles as before.
sim = rebound.Simulation() sim.integrator = "ias15" sim.add(m=1.) sim.add(m=1e-6,a=1.) sim.move_to_com() # Moves to the center of momentum frame
But we change the additional force to be
ps = sim.particles tau = 1000. def migrationForce(reb_sim): ps.ax -= ps.vx/tau ps.ay -= ps.vy/tau ps.az -= ps.vz/tau
We need to let REBOUND know that our force is velocity dependent. Otherwise, REBOUND will not update the velocities of the particles.
sim.additional_forces = migrationForce sim.force_is_velocity_dependent = 1
Now, we integrate as before. But this time we keep track of the semi-major axis instead of the eccentricity.
Nout = 1000 a_s = np.zeros(Nout) times = np.linspace(0.,100.*2.*np.pi,Nout) for i, time in enumerate(times): sim.integrate(time) a_s[i] = sim.particles.a fig = plt.figure(figsize=(15,5)) ax = plt.subplot(111) ax.set_xlabel("time") ax.set_ylabel("semi-major axis") plt.plot(times, a_s);
The semi-major axis decays exponentially on a timescale
In the above example, REBOUND is calling a python function at every timestep. This can be slow. Note that you can also set
rebound.additional_forces to a c function pointer. This let's you speed up the simulation significantly. However, you need to write your own c function/library that knows how to calculate the forces. Or, you use Dan Tamayo's new migration library (in preparation).