Integrators
Numerical integrators are the backbone of any Nbody package. A numerical integrator evolves particles forward in time, one timestep at a time. To do that, the integrator needs to know the current position and velocity coordinates of the particles, and the equations of motion which come in the form of a set of ordinary differential equations.
Because an exact solution to these differential equations is in general unknown, each integrator attempts to approximate the true solution numerically. Different integrators do this differently and each of them has some advantages and some disadvantages. Each of the builtin integrators of REBOUND is described in this section.
IAS15
IAS15 stands for Integrator with Adaptive Stepsize control, 15th order. It is a very high order, nonsymplectic integrator which can handle arbitrary forces (including those who are velocity dependent). It is in most cases accurate down to machine precision (16 significant decimal digits). The IAS15 implementation in REBOUND can integrate variational equations. The algorithm is described in detail in Rein & Spiegel 2015 and also in the original paper by Everhart 1985.
IAS15 is the default integrator of REBOUND, so if you want to use it, you don't need to do anything. However, you can also set it explicitly:
struct reb_simulation* r = reb_simulation_create();
r>integrator = REB_INTEGRATOR_IAS15;
sim = rebound.Simulation()
sim.integrator = "ias15"
The setting for IAS15 are stored in the reb_integrator_ias15
structure.
epsilon
(double
)
IAS15 is an adaptive integrator. It chooses its timesteps automatically. This parameter controls the accuracy of the integrator. The default value is \(10^{9}\). Setting this parameter to 0 turns off adaptive timestepping and a constant timestep will is used. Turning off adaptive timestepping is rarely useful.
Important
It is tempting to change
epsilon
to achieve a speedup at the loss of some accuracy. However, that makes rarely sense. The reason is that IAS15 is a very high (15th!) order integrator. Suppose we increase the timestep by a factor of 10. This will increase the error by a factor of \(10^{15}\). In other words, a simulation that previously was converged to machine precision will now have an error of order unity. min_dt
(double
)
This sets the minimum allowed timestep. The default value is 0. Set this to a finite value if the adaptive timestep becomes excessively small, for example during close encounters or because of finite floating point precision. Use with caution and make sure the simulation results still make physically sense as you might be in danger of ignoring small timescales in the problem. The following code sets the smallest timestep to \(10^{3}\) time units:
struct reb_simulation* r = reb_simulation_create(); r>ri_ias15.min_dt = 1e3;
sim = rebound.Simulation() sim.ri_ias15.min_dt = 1e3
adaptive_mode
(unsigned int
) This flag determines how the adaptive timestep is chosen. The previous name of this flag was
epsilon_global
. The default is 2 which corresponds to the timestep criterion described in Pham, Rein, and Spiegel (2024). This should be optimal in almost all cases. If set to 0, the fractional error is estimated viamax(acceleration_error/acceleration)
and the timestep criterion of Rein and Spiegel (2015) is used. If set to 1, IAS15 estimates the fractional error viamax(acceleration_error)/max(acceleration)
where the maximum is taken over all particles. As before, the timestep criterion of Rein and Spiegel (2015) is used. This was the default until January 2024. If set to 3, then the criterion of Aarseth 1985 is used.
All other members of this structure are only for internal IAS15 use.
WHFast
WHFast is an implementation of the symplectic WisdomHolman integrator. It is the best choice for systems in which there is a dominant central object and perturbations to the Keplerian orbits are small. It supports first and second symplectic correctors as well as the kernel method of Wisdom et al. 1996 with various different kernels. The basic implementation of WHFast is described in detail in Rein & Tamayo 2015. The higher order aspects of it are described in Rein, Tamayo & Brown 2019. WHFast also supports first order variational equations which can be used in chaos estimators (Rein & Tamayo 2016). The user can choose between Jacobi and Democratic Heliocentric coordinates.
The following code enables the WHFast integrator. Because WHFast is not an adaptive integrator, you also need to set a timestep. Typically, this should be a small fraction (a few percent) of the smallest dynamical timescale in the problem.
struct reb_simulation* r = reb_simulation_create();
r>integrator = REB_INTEGRATOR_WHFAST;
r>dt = 0.1;
sim = rebound.Simulation()
sim.integrator = "whfast"
sim.dt = 0.1
The setting for WHFast are stored in the reb_integrator_whfast
structure, which itself is part of the simulation structure.
unsigned int corrector

This variable turns on/off different first symplectic correctors for WHFast. By default, it is set to zero and symplectic correctors are turned off.
First symplectic correctors remove error terms up to \(O(\epsilon \cdot dt^p)\), where \(p\) is the order of the symplectic corrector, and \(\epsilon\) is the mass ratio in the system. The following first correctors are implemented in REBOUND:
Order Number of stages 0 Correctors turned off (default) 3 2 5 4 7 6 11 10 17 16 For most cases you want to choose the 17th order corrector. You only want to consider lower order correctors if frequent outputs are required and speed is an issue. Symplectic correctors are turned on as follows.
r>ri_whfast.corrector = 17; r>ri_whfast.safe_mode = 0;
sim.ri_whfast.corrector = 17 sim.ri_whfast.safe_mode = 0
Note that the above code also turns off the safe mode. Most likely, you want to do that too (see below for a description of the safe mode).
unsigned int corrector2

This variable turns on/off second symplectic correctors for WHFast. By default, second symplectic correctors are off (0). Set to 1 to use second symplectic correctors.
Info
The nomenclature can be a bit confusing. First symplectic correctors are different from second symplectic correctors. And in REBOUND first symplectic correctors have different orders (see above). Second symplectic correctors on the other hand can only be turned on or off. See Rein, Tamayo & Brown 2019 for more on high order symplectic integrators.
unsigned int kernel

This variable determines the kernel of the WHFast integrator. The following options are currently supported:
 The standard WisdomHolman kick step. This is the default.
 Exact modified kick. This works for Newtonian gravity only. Not additional forces.
 The composition kernel.
 Lazy implementer's modified kick. This is often the best option.
Check Rein, Tamayo & Brown 2019 for details on what these kernel methods are. The syntax to use them is
r>ri_whfast.kernel = REB_WHFAST_KERNEL_DEFAULT; // or r>ri_whfast.kernel = REB_WHFAST_KERNEL_MODIFIEDKICK; // or r>ri_whfast.kernel = REB_WHFAST_KERNEL_COMPOSITION; // or r>ri_whfast.kernel = REB_WHFAST_KERNEL_LAZY;
sim.ri_whfast.kernel = "default" # or sim.ri_whfast.kernel = "modifiedkick" # or sim.ri_whfast.kernel = "composition" # or sim.ri_whfast.kernel = "lazy"
unsigned int coordinates

WHFast supports different coordinate systems. Default are Jacobi Coordinates. Other options are democratic heliocentric coordinates, and the WHDS coordinates (Hernandez & Dehnen, 2017) The syntax to use them is
r>ri_whfast.coordinates = REB_WHFAST_COORDINATES_JACOBI; // or r>ri_whfast.coordinates = REB_WHFAST_COORDINATES_DEMOCRATICHELIOCENTRIC; // or r>ri_whfast.coordinates = REB_WHFAST_COORDINATES_WHDS;
sim.ri_whfast.coordinates = "jacobi" # or sim.ri_whfast.coordinates = "democraticheliocentric" # or sim.ri_whfast.coordinates = "whds"
unsigned int recalculate_coordinates_this_timestep
 Setting this flag to one will recalculate the internal coordinates from the particle structure in the next timestep. After the timestep, the flag gets set back to 0. If you want to change particles after every timestep, you also need to set this flag to 1 before every timestep. Default is 0.
unsigned int safe_mode
 If this flag is set (the default), WHFast will recalculate the internal coordinates (Jacobi/heliocentric/WHDS) and synchronize every timestep, to avoid problems with outputs or particle modifications between timesteps. Setting it to 0 will result in a speedup, but care must be taken to synchronize and recalculate the internal coordinates when needed. See also the AdvWHFast.ipynb tutorial.
unsigned int keep_unsynchronized
 This flag determines if the inertial coordinates generated are discarded in subsequent timesteps (cached Jacobi/heliocentric/WHDS coordinates are used instead). The default is 0. Set this flag to 1 if you require outputs and bitwise reproducibility
All other members of the reb_integrator_whfast
structure are for internal use only.
GraggBulirschStoer (BS)
The GraggBulirschStoer integrator (short BS for BulirschStoer) is an adaptive integrator which uses Richardson extrapolation and the modified midpoint method to obtain solutions to ordinary differential equations.
The version in REBOUND is based on the method described in Hairer, Norsett, and Wanner 1993 (see section II.9, page 224ff), specifically the JAVA implementation available in the Hipparchus package. The Hipparchus as well as the REBOUND version are adaptive in both the timestep and the order of the method for optimal performance. The BS implementation in REBOUND can integrate first and second order variational equations.
The BS integrator is particularly useful for short integrations where only medium accuracy is required. For long integrations a symplectic integrator such as WHFast performs better. For high accuracy integrations the IAS15 integrator performs better. Because BS is adaptive, it can handle close encounters. Currently a collision search is only performed after every timestep, i.e. not after a subtimestep.
The following code enables the BS integrator and sets both the relative and absolute tolerances to 0.0001 (the default is \(10^{8}\)):
struct reb_simulation* r = reb_simulation_create();
r>integrator = REB_INTEGRATOR_BS;
r>ri_bs.eps_rel = 1e4;
r>ri_bs.eps_abs = 1e4;
sim = rebound.Simulation()
sim.integrator = "BS"
sim.ri_bs.eps_rel = 1e4
sim.ri_bs.eps_abs = 1e4
The BS integrator tries to keep the error of each coordinate \(y\) below \(\epsilon_{abs} + \epsilon_{rel} \cdot \lefty\right\). Note that this applies to both position and velocity coordinates of all particles which implies that the code units you're choosing for the integration matter. If you need fine control over the scales used internally, you can set the getscale
function pointer in r>ri_bs.nbody_ode
(this is currently undocumented, search the source code for getscale
to find out more).
Info
The code does not guarantee that the errors remain below the tolerances. In particular, note that BS is not a symplectic integrator which results in errors growing linearly in time (phase errors grow quadratically in time). It requires some experimentation to find the tolerances that offer the best compromise between accuracy and speed for your specific problem.
You can limit the timestep with both a maximum and minimum timestep:
r>ri_bs.min_dt = 1e5;
r>ri_bs.max_dt = 1e2;
sim.ri_bs.min_dt = 1e5
sim.ri_bs.max_dt = 1e2
Compared to the other integrators in REBOUND, BS can be used to integrate arbitrary ordinary differential equations (ODEs), not just the Nbody problem. We expose an ODEAPI in REBOUND which allows you to make use of this. Userdefined ODEs are always integrated with BS. You can choose to integrate the Nbody equations with BS as well, or any of the other integrators.
If you choose BS for the Nbody equations, then BS will treat all ODEs (Nbody + all userdefined ones) as one big system of coupled ODEs. This means your timestep will be set by either the Nbody problem or the userdefined ODEs, whichever involves the shorter timescale.
If you choose IAS15 or WHFast for the Nbody equation but also have userdefined ODEs, then they cannot be treated as one big coupled system of ODEs anymore. In that case the Nbody integration is done first. Then the userdefined ODEs are advanced to the exact same time as the Nbody system using BS using whatever timestep is required to achieve the tolerance set in the ri_bs
struct. During the integration of the userdefined ODEs, the coordinates of the particles in the Nbody simulation are assumed to be fixed at their final position and velocity. This introduces an error. However, if the system evolves adiabatically (the timescales in the userdefined ODEs are much longer than in the Nbody problem), then the error will be small.
The following code sets up a REBOUND simulation in which a harmonic oscillator is driven by the phase of a planet orbiting a star:
void derivatives(struct reb_ode* const ode, double* const yDot, const double* const y, const double t){
struct reb_orbit o = reb_orbit_from_particle(ode>r>G, ode>r>particles[1], ode>r>particles[0]);
const double omega = 1;
double forcing = sin(o.f);
yDot[0] = y[1];
yDot[1] = omega*omega*y[0] + forcing;
}
void run(){
struct reb_simulation* r = reb_simulation_create();
reb_simulation_add_fmt(r, "m", 1.);
reb_simulation_add_fmt(r, "m a e", 1e3, 1., 0.1);
r>integrator = REB_INTEGRATOR_BS;
struct reb_ode* ho = reb_ode_create(r,2); // Add an ODE with 2 dimensions
ho>derivatives = derivatives; // Right hand side of the ODE
ho>y[0] = 1; // Initial conditions
ho>y[1] = 0;
}
import numpy as np
def derivatives(ode, yDot, y, t):
omega = 1.0
sim_pointer = ode.contents.r
orbit = sim_pointer.contents.particles[1]
forcing = np.sin(orbit.f)
yDot[0] = y[1]
yDot[1] = omega*omega*y[0] + forcing
sim = rebound.Simulation()
sim.add(m=1)
sim.add(m=1e3, a=1, e=0.1)
sim.integrator = "BS"
ho = sim.create_ode(length=2) # Add an ODE with 2 dimensions
ho.derivatives = derivatives # Right hand side of the ODE
ho.y[0] = 1.0 # Initial conditions
ho.y[1] = 0.0
Mercurius
MERCURIUS is a hybrid symplectic integrator very similar to MERCURY (Chambers 1999).
It uses WHFast for long term integrations but switches over smoothly to IAS15 for close encounters.
The MERCURIUS implementation is described in Rein et al 2019.
The following code enables MERCURIUS and sets the critical radius to 4 Hill radii
struct reb_simulation* r = reb_simulation_create();
r>integrator = REB_INTEGRATOR_MERCURIUS;
r>ri_mercurius.r_crit_hill = 4.;
sim = rebound.Simulation()
sim.integrator = "mercurius"
sim.ri_mercurius.r_crit_hill = 4.
The reb_integrator_mercurius
structure contains the configuration and data structures used by the hybrid symplectic MERCURIUS integrator.
double (*L) (const struct reb_simulation* const r, double d, double dcrit)

This is a function pointer to the force switching function. If NULL (the default), the MERCURY switching function will be used. The argument
d
is the distance between two particles. The argumentdcrit
is the maximum critical distances of the two particles. The return value is a scalar between 0 and 1. If this function always returns 1, then the integrator effectively becomes the standard WisdomHolman integrator.The following switching functions are available:

Mercury switching function
This is the same polynomial switching function as used in MERCURY.
 Smooth switching functionsdouble reb_integrator_mercurius_L_mercury(const struct reb_simulation* const r, double d, double dcrit);
These two polynomials switching functions are 4 and 5 times differentiable. Using smooth switching functions can improve the accuracy. For a detailed discussion see Hernandez 2019.
double reb_integrator_mercurius_L_C4(const struct reb_simulation* const r, double d, double dcrit); double reb_integrator_mercurius_L_C5(const struct reb_simulation* const r, double d, double dcrit);

Infinitely differentiable switching function
This is an infinitely differentiable switching function.
double reb_integrator_mercurius_L_infinity(const struct reb_simulation* const r, double d, double dcrit);
The switching function can be set using this syntax:
struct reb_simulation* r = reb_simulation_create(); r>ri_mercurius.L = reb_integrator_mercurius_L_infinity;
sim = rebound.Simulation() sim.ri_mercurius.L = "infinity"

double r_crit_hill
 The critical switchover radii of particles are calculated automatically based on multiple criteria. One criterion calculates the Hill radius of particles and then multiplies it with the
r_crit_hill
parameter. The parameter is in units of the Hill radius. The default value is 3. unsigned int recalculate_coordinates_this_timestep
 Setting this flag to one will recalculate heliocentric coordinates from the particle structure at the beginning of the next timestep. After a single timestep, the flag gets set back to 0. If one changes a particle manually after a timestep, then one needs to set this flag to 1 before the next timestep.
unsigned int recalculate_r_crit_this_timestep
 Setting this flag to one will recalculate the critical switchover distances dcrit at the beginning of the next timestep. After one timestep, the flag gets set back to 0. If you want to recalculate
dcrit
at every timestep, you also need to set this flag to 1 before every timestep. unsigned int safe_mode
 If this flag is set to 1 (the default), the integrator will recalculate heliocentric coordinates and synchronize after every timestep to avoid problems with outputs or particle modifications between timesteps. Setting this flag to 0 will result in a speedup, but care must be taken to synchronize and recalculate coordinates manually if needed.
TRACE
TRACE is a hybrid timereversible integrator, based on the algorithm described in Hernandez & Dehnen 2023. It uses WHFast for long term integrations but switches timereversibly to BS or IAS15 for all close encounters. TRACE is appropriate for systems with a dominant central mass that will occasionally have close encounters. A paper describing the TRACE implementation is in preparation.
The following code enables TRACE and sets the critical radius to 4 Hill radii
struct reb_simulation* r = reb_create_simulation();
r>integrator = REB_INTEGRATOR_TRACE;
r>ri_trace.r_crit_hill = 4;
sim = rebound.Simulation()
sim.integrator = "trace"
sim.ri_trace.r_crit_hill = 4
The reb_integrator_trace
structure contains the configuration and data structures used by the hybrid symplectic TRACE integrator.
int (*S) (const struct reb_simulation* const r, const unsigned int i, const unsigned int j)

This is a function pointer to the switching function for close encounters between noncentral bodies. If NULL (the default), the default switching function will be used. The arguments
i
andj
are the indices of the two particles considered. The return value is either 0 or 1. A return value of 1 means a close encounter has been flagged. If the return values of both this function and the central switching function below are always 0, then the integrator effectively becomes the standard WisdomHolman integrator.
Default switching function
This is a similar (but slightly modified) switching function used in MERCURY. It uses a modified Hill radius criteria, with heliocentric distance replacing semimajor axis
int reb_integrator_trace_switch_default(const struct reb_simulation* const r, const unsigned int i, const unsigned int j);
The switching function can be manually set using this syntax:
struct reb_simulation* r = reb_create_simulation(); r>ri_trace.S = reb_integrator_trace_switch_default;
sim = rebound.Simulation() sim.ri_trace.S = "default"

int (*S_peri) (const struct reb_simulation* const r, const unsigned int j)

This is a function pointer to the switching function for close encounters involving the central body. If NULL (the default), the default switching function will be used. The argument
j
is the index of the noncentral particle considered. The return value is either 0 or 1. A return value of 1 means a close encounter has been flagged.
Default switching function
This switching function checks if a body is close to its pericenter by considering a timescale derived from highorder derivatives of the particle's herliocentric position, inspired by Pham, Rein, and Spiegel 2024.
 Fdot switching functionint reb_integrator_trace_switch_peri_default(const struct reb_simulation* const r, const unsigned int j);
This switching function checks if a body is close to its pericenter by measuring the rate of change of true anomaly, inspired by Wisdom 2015
int reb_integrator_trace_switch_peri_fdot(const struct reb_simulation* const r, const unsigned int j);

Distance switching function
This switching function checks for close encounters with a simple heliocentric distance check
int reb_integrator_trace_switch_peri_distance(const struct reb_simulation* const r, const unsigned int j);
The switching function can be manually set using this syntax:
struct reb_simulation* r = reb_create_simulation(); r>ri_trace.S_peri = reb_integrator_trace_switch_peri_default; // default r>ri_trace.S_peri = reb_integrator_trace_switch_peri_fdot; r>ri_trace.S_peri = reb_integrator_trace_switch_peri_distance; r>ri_trace.S_peri = reb_integrator_trace_switch_peri_none; // Turn off pericenter switching
sim = rebound.Simulation() sim.ri_trace.S_peri = "default" # Following Pham et al 2024 sim.ri_trace.S_peri = "fdot" sim.ri_trace.S_peri = "distance" sim.ri_trace.S_peri = "none" # Turn off pericenter switching

double r_crit_hill
 The critical switchover radii of noncentral particles are calculated based on a modified Hill radii criteria. This modified Hill radius for each particle is calculated and then multiplied by the
hillfac
parameter. The parameter is in units of the modified Hill radius. This value is used by thedefault
switching function. The default value is 4. double peri_crit_eta
 The criteria for a pericenter approach with the central body. This criteria is used in the
default
pericenter switching condition. It flags a particle as in a close pericenter approach if the ratio of the timestep to the condition described in Pham, Rein, and Spiegel 2024. The default value is 1. double peri_crit_fdot
 The criteria for a pericenter approach with the central body. This criteria is used in the
fdot
pericenter switching condition. It flags a particle as in a close pericenter approach if the effective period at pericenter (in units of the timestep) is less than theperi_crit_fdot
parameter. The default value is 17. double peri_crit_distance

The criteria for a pericenter approach with the central body. This criteria is used in the
distance
pericenter switching condition. It flags a particle as in a close pericenter approach if its heliocentric distance is less than theperi_crit_distance
parameter.These switching criteria can be manually set using this syntax:
struct reb_simulation* r = reb_create_simulation(); r>ri_trace.peri_crit_eta = 0.5; // or r>ri_trace.peri_crit_fdot = 16; // or r>ri_trace.peri_crit_distance = 1;
sim = rebound.Simulation() sim.ri_trace.peri_crit_eta = 0.5 # or sim.ri_trace.peri_crit_fdot = 16 # or sim.ri_trace.peri_crit_distance = 1 # or
unsigned int peri_mode
 This variable determines how TRACE integrates close approaches with the central star. The following options are currently supported:
 Integrating the entire system with BS. This is the default.
 Integrating only the Kepler Step with BS.
 Integrating the entire system with IAS15.
Check Lu, Hernandez & Rein 2024 for details on what these kernel methods are. The syntax to use them is
r>ri_trace.peri_mode = REB_TRACE_PERI_PARTIAL_BS; // or r>ri_trace.peri_mode = REB_TRACE_PERI_FULL_BS; // or r>ri_trace.peri_mode = REB_TRACE_PERI_FULL_IAS15; // or
sim.ri_trace.peri_mode = "PARTIAL_BS" # or sim.ri_trace.peri_mode = "FULL_BS" # or sim.ri_trace.peri_mode = "FULL_IAS15" # or
SABA
SABA are symplectic integrators developed by Laskar & Robutel 2001 and Blanes et al. 2013. The implementation in REBOUND supports SABA1, SABA2, SABA3, and SABA4 as well as the corrected versions SABAC1, SABAC2, SABAC3, and SABAC4. Different correctors can be selected. In addition, the following methods with various generalized orders are supported: SABA(8,4,4), SABA(8,6,4), SABA(10,6,4). See Rein, Tamayo & Brown 2019 for details on how these methods work.
The reb_integrator_saba
structure contains the configuration and data structures used by the SABA integrator family.
unsigned int type

This parameter specifies which SABA integrator type is used. The following SABA integrators are supported:
Numerical value C constant name Description 0x0 REB_SABA_1
SABA1 (WisdomHolman) 0x1 REB_SABA_2
SABA2 0x2 REB_SABA_3
SABA3 0x3 REB_SABA_4
SABA4 0x100 REB_SABA_CM_1
SABACM1 (Modified kick corrector) 0x101 REB_SABA_CM_2
SABACM2 (Modified kick corrector) 0x102 REB_SABA_CM_3
SABACM3 (Modified kick corrector) 0x103 REB_SABA_CM_4
SABACM4 (Modified kick corrector) 0x200 REB_SABA_CL_1
SABACL1 (lazy corrector) 0x201 REB_SABA_CL_2
SABACL2 (lazy corrector) 0x202 REB_SABA_CL_3
SABACL3 (lazy corrector) 0x203 REB_SABA_CL_4
SABACL4 (lazy corrector) 0x4 REB_SABA_10_4
SABA(10,4), 7 stages 0x5 REB_SABA_8_6_4
SABA(8,6,4), 7 stages 0x6 REB_SABA_10_6_4
SABA(10,6,4), 8 stages, default 0x7 REB_SABA_H_8_4_4
SABAH(8,4,4), 6 stages 0x8 REB_SABA_H_8_6_4
SABAH(8,6,4), 8 stages 0x9 REB_SABA_H_10_6_4
SABAH(10,6,4), 9 stages SABA(10,6,4) is the default integrator. It has a generalized order of \(O(\epsilon dt^{10} + \epsilon^2 dt^6 + \epsilon^3 dt^4)\).
Below is an example on how to enable the SABA integrators in REBOUND and set a specific type.
struct reb_simulation* r = reb_simulation_create(); r>integrator = REB_INTEGRATOR_SABA; r>ri_saba.type = REB_SABA_10_6_4;
One can also use the following shorthand:sim = rebound.Simulation() sim.integrator = "saba" sim.ri_saba.type = "(10,6,4)"
sim = rebound.Simulation() sim.integrator = "SABA(10,6,4)"
unsigned int safe_mode
 This flag has the same functionality as in WHFast. Default is 1. Setting this to 0 will provide a speedup, but care must be taken with synchronizing integration steps and modifying particles.
unsigned int keep_unsynchronized
 This flag determines if the inertial coordinates generated are discarded in subsequent timesteps (cached Jacobi coordinates are used instead). The default is 0. Set this flag to 1 if you require outputs and bitwise reproducibility
JANUS
Janus is a bitwise timereversible highorder symplectic integrator using a mix of floating point and integer arithmetic. It is described in Rein & Tamayo 2018.
The following code shows how to enable JANUS and set the length and velocity scales.
struct reb_simulation* r = reb_simulation_create();
r>integrator = REB_INTEGRATOR_JANUS;
r>ri_janus.scale_pos = 1e10;
r>ri_janus.scale_vel = 1e10;
sim = rebound.Simulation()
sim.integrator = "janus"
sim.ri_janus.scale_pos = 1e10
sim.ri_janus.scale_vel = 1e10
The reb_integrator_janus
structure contains the configuration and data structures used by the bibwise reversible JANUS integrator.
double scale_pos
 Scale of the problem. Positions get divided by this number before the conversion to an integer. Default: \(10^{16}\).
double scale_vel
 Scale of the problem. Velocities get divided by this number before the conversion to an integer. Default: \(10^{16}\).
unsigned int order
 The order of the scheme. Default is 6.
unsigned int recalculate_integer_coordinates_this_timestep
 If this flag is set, then JANUS will recalculate the integer coordinates from floating point coordinates at the next timestep.
All other members of this structure are only for internal use and should not be changed manually.
Embedded Operator Splitting Method (EOS)
This is the Embedded Operator Splitting (EOS) methods described in Rein 2019.
The reb_integrator_eos
structure contains the configuration and data structures used by EOS.
unsigned int phi0
 Outer operator splitting scheme (see below for options)
unsigned int phi1
 Inner operator splitting scheme (see below for options)
unsigned int n
 Number of subtimesteps. Default: 2.
unsigned int safe_mode
 If set to 0, always combine drift steps at the beginning and end of
phi0
. If set to 1,n
needs to be bigger than 1.
The following operator splitting methods for phi0
and phi1
are supported in the EOS integrator.
Numerical value  Constant name  Description 

0x00  REB_EOS_LF 
2nd order, standard leapfrog 
0x01  REB_EOS_LF4 
4th order, three function evaluations 
0x02  REB_EOS_LF6 
6th order, nine function evaluations 
0x03  REB_EOS_LF8 
8th order, seventeen function evaluations, see Blanes & Casa (2016), p91 
0x04  REB_EOS_LF4_2 
generalized order (4,2), two force evaluations, McLachlan 1995 
0x05  REB_EOS_LF8_6_4 
generalized order (8,6,4), seven force evaluations 
0x06  REB_EOS_PLF7_6_4 
generalized order (7,6,4), three force evaluations, pre and postprocessors 
0x07  REB_EOS_PMLF4 
4th order, one modified force evaluation, pre and postprocessors, Blanes et al. (1999) 
0x08  REB_EOS_PMLF6 
6th order, three modified force evaluations, pre and postprocessors, Blanes et al. (1999) 
The following code shows how to enable EOS and set the embedded methods.
struct reb_simulation* r = reb_simulation_create();
r>integrator = REB_INTEGRATOR_EOS;
r>ri_eos.phi0 = REB_EOS_LF4;
r>ri_eos.phi1 = REB_EOS_LF4;
r>ri_eos.n = 6;
sim = rebound.Simulation()
sim.integrator = "eos"
sim.ri_eos.phi0 = "LF4"
sim.ri_eos.phi1 = "LF4"
sim.ri_eos.n = 6
Leapfrog
REB_INTEGRATOR_LEAPFROG
This is the standard leap frog integrator. It is second order and symplectic. No configuration is available (the timestep is set in the simulation structure).
Symplectic Epicycle Integrator (SEI)
REB_INTEGRATOR_SEI
Symplectic Epicycle Integrator (SEI), mixed variable symplectic integrator for the shearing sheet, second order, Rein & Tremaine 2011. The reb_integrator_sei
structure contains the configuration and data structures used by the Symplectic Epicycle Integrator (SEI).
double OMEGA

Epicyclic/orbital frequency. This can be set as follows:
struct reb_simulation* r = reb_simulation_create(); r>integrator = REB_INTEGRATOR_SEI; r>ri_sei.OMEGA = 1.0;
sim = rebound.Simulation() sim.integrator = "sei" sim.ri_sei.OMEGA = 1.0
double OMEGAZ
 Epicyclic frequency in vertical direction. Defaults to
OMEGA
if not set.
All other members of this structure are only for internal use and should not be changed manually.
No integrator
Sometimes it might make sense to simply not advance any particle positions or velocities. By selecting this integrator, one can still perform integration steps, but particles will not move.
Here is how to do that:
struct reb_simulation* r = reb_simulation_create();
r>integrator = REB_INTEGRATOR_NONE;
sim = rebound.Simulation()
sim.integrator = "none"
WHFast512
WHFast512 is a symplectic WisdomHolman integrator. It is using Single Instruction Multiple Data (SIMD) parallelism and 512bit Advanced Vector Extensions (AVX512) to speed up the integration of planetary systems by up to 4.7x compared to the standard version of WHFast.
Important
To use WHFast512 you need to compile and run REBOUND on a computer that has a CPU which supports AVX512 instructions. You will see an error message if you try to use WHFast512 but have not compiled REBOUND with the AVX512 flag. We describe below how to do this for both the C and python versions of REBOUND below.
To find out if your CPU supports AVX512 instructions, check for the AVX512 flags by running
cat /proc/cpuinfo  grep avx512
Note that you can read Simulationarchives of simulations which used WHFast512 on machines that do not support AVX512 instruction. If a synchronization is required, it will be performed with the standard WHFast integrator.
To turn on the AVX512 flag, go to the Makefile in problem directory. Add this line at the top:
export AVX512=1
export AVX512=0
march=native
flag to the compiler. This will optimize your code (and enable AVX512 instruction) for the specific CPU you're using.
export OPT=march=native
make clean
make
To use WHFast512 from python, you need to compile REBOUND with AVX512 instructions enabled. They are disabled by default and enabled with the AVX512 environment variable. To install the latest release of REBOUND on pypi use:
export AVX512=1
pip install rebound
export AVX512=1
pip install e .
Once you have compiled REBOUND with AVX512 enabled, you can use WHFast512 like any other integrator:
struct reb_simulation* r = reb_simulation_create();
r>integrator = REB_INTEGRATOR_WHFAST512;
sim = rebound.Simulation()
sim.integrator = "whfast512"
See also this example on how to use WHFast512. If you are interested in integrating 2 or 4 planet systems in parallel, see this example.
To allow for the best performance, WHFast512 has certain limitations that WHFast does not have.
 The number of particles cannot exceed 9 (1 star and 8 planets) and needs to be constant.
 Although you can use WHFast512 with any number of planets (up to 8), the performance is best if the system has either 2, 4, or 8 planets.
 The gravitational constant needs to be exactly equal to 1. Note that you can always rescale your system such that G=1.
 The integrator always combines the first and last drift step (
safe_mode=0
for WHFast).  No variational or test particles are supported (although a particle can have mass 0).
 MEGNO and other chaos indicators are not supported.
 WHFast512 always uses democratic heliocentric coordinates. Jacobi coordinates are not supported.
 The timestep needs to be constant and the
exact_finish_time
flag needs to be set to 0. To change the timestep, first synchronize the simulation, then callreb_simulation_reset_integrator()
.  The masses of all particles need to be constant. To change the masses, first synchronize the simulation, then call
reb_simulation_reset_integrator()
.  Additional forces (other than the GR potential) and REBOUNDx are not supported.
The setting for WHFast512 are stored in the reb_integrator_whfast512
structure, which itself is part of the simulation structure.
The following settings are available:
unsigned int keep_unsynchronized
 This flag determines if democratic heliocentric coordinates are reused after subsequent calls to
reb_simulation_integrate()
. The default is 0. This makes WHFast512 recalculate democratic heliocentric coordinates at the beginning of eachreb_simulation_integrate()
call. Set this flag to 1 if you want to continue an integration using unsynchronized democratic heliocentric coordinates. This is useful if you require outputs (and therefore synchronization) but don't want the integration to be affected by the output to allow for bitwise reproducibility. unsigned int gr_potential
 This flag determines if an additional \(1/r^2\) potential is included in the force calculation. The default is 0. Set to 1 to turn on the potential. This can be used to mimic general relativistic precession. Note that this feature assumes units of AU and year/2pi.
unsigned int N_systems
 This flag determines how many systems are integrated in parallel. Possible values are 1, 2, or 4. By default this is set to 1 which means WHFast512 is integrating only one system at a time. If your system has fewer than 5 planets, then you can use WHFast512 to integrate 2 systems in parallel. If your system has fewer than 3 planets, then you can use WHFast512 to integrate 4 systems in parallel. If multiple systems are integrated at the same time, particles must be added in the following order: Star 1, Planet, Planet, Star 2, Planet, Planet, ... For more information see the this example.