Quick User Quide (C)

This section describes the C version of REBOUND.

Installation

You can download, compile and run REBOUND on almost any modern operating system within seconds. Simply copy and paste this line to your terminal and press enter:

git clone http://github.com/hannorein/rebound && cd rebound/examples/shearing_sheet && make && ./rebound

or if you do not have git installed:

wget --no-check-certificate https://github.com/hannorein/rebound/tarball/master -O- | tar xvz && cd hannorein-rebound-*/examples/shearing_sheet/ && make && ./rebound

Make sure you have a compiler suite installed. Open a terminal and type make and cc to test if your installation is complete. If you are on OSX, you can download Xcode from the AppStore (for free). Once installed, open Xcode, go to Settings, then Downloads and install the Command Line Tools.

Note: REBOUND does not work on Windows, and we currently do not have plans to support it.

Code structure

REBOUND can be used as a shared library. However, installing a system-wide shared library can sometimes be an obstacle for new users, especially if you want to change the code frequently or don’t have root access. For that reason, all the examples can be compiled by simply typing make in any of the example directories.

Let’s look at how to setup a simple REBOUND simulation:

#include "rebound.h"

int main(int argc, char* argv[]) {
        struct reb_simulation* r = reb_create_simulation();
        r->dt = 0.1;
        r->integrator = REB_INTEGRATOR_WHFAST;

        struct reb_particle p1 = {0};
        p1.m = 1.;
        reb_add(r, p1);

        struct reb_particle p2 = {0};
        p2.x = 1;
        p2.vy = 1;
        p2.m = 0.;
        reb_add(r, p2);

        reb_move_to_com(r);
        reb_integrate(r,100.);
}

In the first line we include the REBOUND header file. This file contains all the declarations of the structures and functions that we will be using.

Next, we declare the only function in our file. It is the standard C main() function. Within that, we first create a reb_simulation structure. This is the main structure that contains all the variables, pointers and particles of a REBOUND simulation. You can create multiple reb_simulation structures at the same time. REBOUND is thread-safe.

We can then set flags and variables in the reb_simulation structure. Note that the r variable is a pointer to the structure, so we use the arrow syntax r-> to set variables. The next line chooses the integrator module. Here, we use the WHFast symplectic integrator.

We then create two particles, both of which are represented by a reb_particle structure. The = {0} syntax ensures that our structs are initialized with zeros. We set the initial conditions (the ones we don’t want to be zero) and then add the particle to the simulation using the reb_add() function. Note that this function takes two arguments, the first one is the simulation to which you want to add the particle, and the second is the particle that you want to add.

Finally, we call the REBOUND function reb_move_to_com(). It moves the particles to a centre of mass reference frame (this prevents particles from drifting away from the origin). We then start the integration. Here, we integrate for 100 time units. By default REBOUND used units in which G=1, thus a particle around an m=1 mass central object at a semi-major axis of 1 needs 2pi time units for one orbit.

Note that all REBOUND functions start with the three character prefix reb_.

Next, let’s add a call-back function to the above example. This function will be called after every timestep and we can use it to output simulation data. The relevant function pointer is called heartbeat in the reb_simulation structure. We first declare and implement the function and then set the pointer in the main routine:

void heartbeat(struct reb_simulation* r){
       printf("%f\n",r->t);
}
int main(int argc, char* argv[]) {
       ...
       r->heartbeat = heartbeat;
       ...
}

As you can probably guess, this will make the program print out the current time after every timestep. Since the heartbeat function receives the reb_simulation structure, you have access to all the variables and particles within the simulation. You don’t need any global variables for that. For example, if we wanted to print out the x coordinate of the 2nd particle (the index starts at 0, so the second particle has index 1), we could use this heartbeat function.

void heartbeat(struct reb_simulation* r){
       double x = r->particles[1].x;
       printf("%f\n",x);
}

REBOUND comes with various built-in output functions that make your life easier. It can for example calculate the orbital elements for you or output to a binary file to save space. The examples are the best way to get to know these functions. You can also look at the rebound.h file in the src/ directory to get an glimpse of the available functions.

Compiling and directory structure

If you look at the examples in the examples/ directory, you see one .c file and one Makefile. All the REBOUND code itself is in the src/ directory. This setup keeps the directory in which you’re working in nice and clean. To compile one of the examples, go to the directory and type make. Then the following events happen

  • The Makefile sets up various environment variables. These determine settings like the compiler optimization flags and which libraries are included (see below).

  • Next, it calls the Makefile in the src/ directory and compiles the entire REBOUND code into a shared library.

  • It then creates a symbolic link from the current directory to the location of the share library in the src directory.

  • Finally it compiles your code, the problem.c file, into an executable file.

You can execute that file with ./rebound. After you edited a file, you can simply type make again to recompile. If you change any of the environment variables, clean the build directiory first, by executing make clean.

Possible issues when compiling REBOUND

REBOUND should be extremely easy to compile as it does not require any external libraries. You might nevertheless run into problems. Two of the most common issues are:

  • Missing compilers. Make sure you have a C compiler installed. If you are using a Mac, install the XCode package which you can download for free on the AppStore.

  • Missing glfw3 library. You can compile REBOUND with support for real-time OpenGL visualizations. This requires the glfw3 library. If you are on a Mac, then the easiest way to install the glfw3 library is with homebrew: brew tap homebrew/versions && brew install glfw3. If you are on Linux, you can install it with your package manager, for example with sudo apt-get install libglfw3-dev. Alternatively, you can disable the OpenGL visualization in the Makefile by setting OPENGL=0. Then, execute make clean and try compiling the program again. On some system the glfw library is called glfw3 instead. In that case, changing -lglfw to -lglfw3 in the file src/Makefile.defs might help.

  • Issue with march native Some users have reported issues related to the compiler flag -march=native which tries to optimize the code for the native architecture. This seems to happen with certain compilers on MacOSX. We removed this flag, which might results in slightly less optimized code. Readding the -march=native or -mtune=native flags, in the file src/Makefile.defs or in setup.py for the python version might help performance in certain cases.

API Documentation

We provide a full API documentation in a separate file. The most important REBOUND API structures and functions are listed below. Note that you can also look at the code itself. The starting point is the rebound.h file in the src/ directory. This is where the public API is defined.

The reb_simulation structure

struct reb_simulation

Main struct encapsulating one entire REBOUND simulation.

This structure contains all variables, status flags and pointers of one REBOUND simulation. To create a REBOUND simulation use the reb_create_simulation() function. This will ensure that all variables and pointers are initialized correctly.

Variables related to time, current number of particles and simulation status/control

double t

Current simulation time.

double G

Gravitational constant. Default: 1.

double softening

Gravitational softening parameter. Default: 0.

double dt

Current timestep.

double dt_last_done

Last dt used by integrator.

unsigned long long steps_done

Timesteps done.

int N

Current number of particles on this node.

int N_var

Total number of variational particles. Default: 0.

int var_config_N

Number of variational configuration structs. Default: 0.

struct reb_variational_configuration *var_config

These configuration structs contain details on variational particles.

int N_active

Number of massive particles included in force calculation (default: N). Particles with index >= N_active are considered testparticles.

int testparticle_type

Type of the particles with an index>=N_active. 0 means particle does not influence any other particle (default), 1 means particles with index < N_active feel testparticles (similar to MERCURY’s small particles). Testparticles never feel each other.

struct reb_hash_pointer_pair *particle_lookup_table

Array of pairs that map particles’ hashes to their index in the particles array.

int hash_ctr

Counter for number of assigned hashes to assign unique values.

int N_lookup

Number of entries in the particle lookup table.

int allocatedN_lookup

Number of lookup table entries allocated.

int allocatedN

Current maximum space allocated in the particles array on this node.

struct reb_particle *particles

Main particle array. This contains all particles on this node.

struct reb_vec3d *gravity_cs

Vector containing the information for compensated gravity summation.

int gravity_cs_allocatedN

Current number of allocated space for cs array.

struct reb_treecell **tree_root

Pointer to the roots of the trees.

int tree_needs_update

Flag to force a tree update (after boundary check)

double opening_angle2

Square of the cell opening angle \( \theta \).

REB_STATUS status

Set to 1 to exit the simulation at the end of the next timestep.

int exact_finish_time

Set to 1 to finish the integration exactly at tmax. Set to 0 to finish at the next dt. Default is 1.

unsigned int force_is_velocity_dependent

Set to 1 if integrator needs to consider velocity dependent forces.

unsigned int gravity_ignore_terms

Ignore the gravity form the central object (1 for WHFast, 2 for WHFastHelio, 0 otherwise)

double output_timing_last

Time when reb_output_timing() was called the last time.

unsigned long display_clock

Display clock, internal variable for timing refreshs.

int save_messages

Set to 1 to ignore messages (used in python interface).

char **messages

Array of strings containing last messages (only used if save_messages==1).

double exit_max_distance

Exit simulation if distance from origin larger than this value.

double exit_min_distance

Exit simulation if distance from another particle smaller than this value.

double usleep

Wait this number of microseconds after each timestep, useful for slowing down visualization.

int track_energy_offset

< Datastructure stores visualization related data. Does not have to be modified by the user.

Track energy change during collisions and ejections (default: 0).

double energy_offset

Energy offset due to collisions and ejections (only calculated if track_energy_offset=1).

double walltime

Walltime in seconds used by REBOUND for this simulation (integration only, not visualization, heartbeat function, etc).

uint32_t python_unit_l

Information only used for when working with units in python.

uint32_t python_unit_m

Information only used for when working with units in python.

uint32_t python_unit_t

Information only used for when working with units in python.

Variables related to ghost/root boxes

struct reb_vec3d boxsize

Size of the entire box, root_x*boxsize.

double boxsize_max

Maximum size of the entire box in any direction. Set in box_init().

double root_size

Size of a root box.

int root_n

Total number of root boxes in all directions, root_nx*root_ny*root_nz. Default: 1. Set in box_init().

int root_nx

Number of root boxes in x direction. Default: 1.

int root_ny

Number of root boxes in y direction. Default: 1.

int root_nz

Number of root boxes in z direction. Default: 1.

int nghostx

Number of ghostboxes in x direction.

int nghosty

Number of ghostboxes in y direction.

int nghostz

Number of ghostboxes in z direction.

Variables related to collision search and detection

int collision_resolve_keep_sorted

Keep particles sorted if collision_resolve removes particles during a collision.

struct reb_collision *collisions

Array of all collisions.

int collisions_allocatedN

Size allocated for collisions.

double minimum_collision_velocity

Used for hard sphere collision model.

double collisions_plog

Keep track of momentum exchange (used to calculate collisional viscosity in ring systems.

double max_radius[2]

Two largest particle radii, set automatically, needed for collision search.

long collisions_Nlog

Keep track of number of collisions.

Variables related to the chaos indicator MEGNO

int calculate_megno

Internal flag that determines if megno is calculated (default=0, but megno_init() sets it to the index of variational particles used for megno)

double megno_Ys

Running megno sum (internal use)

double megno_Yss

Running megno sum (internal use)

double megno_cov_Yt

covariance of MEGNO Y and t

double megno_var_t

variance of t

double megno_mean_t

mean of t

double megno_mean_Y

mean of MEGNO Y

long megno_n

number of covariance updates

Variables related to SimulationArchive

int simulationarchive_version

Version of the SA binary format (1=original/, 2=incremental)

long simulationarchive_size_first

(Deprecated SAV1) Size of the initial binary file in a SA

long simulationarchive_size_snapshot

(Deprecated SAV1) Size of a snapshot in a SA (other than 1st), in bytes

double simulationarchive_auto_interval

Current sampling cadence, in code units.

double simulationarchive_auto_walltime

Current sampling cadence, in wall time.

unsigned long long simulationarchive_auto_step

Current sampling cadence, in time steps.

double simulationarchive_next

Next output time (simulation tim or wall time, depending on wether auto_interval or auto_walltime is set)

unsigned long long simulationarchive_next_step

Next output step (only used if auto_steps is set)

char *simulationarchive_filename

Name of output file.

Variables describing the current module selection

enum [anonymous]

Available collision routines.

Values:

REB_VISUALIZATION_NONE = 0

No visualization (default if OPENGL compiler flag is turned off)

REB_VISUALIZATION_OPENGL = 1

OpenGL visualization (default if OPENGL compiler flag is turned on)

REB_VISUALIZATION_WEBGL = 2

WebGL visualization, only usable from Jupyter notebook widget.

enum [anonymous]

Available collision routines.

Values:

REB_COLLISION_NONE = 0

Do not search for collisions (default)

REB_COLLISION_DIRECT = 1

Direct collision search O(N^2)

REB_COLLISION_TREE = 2

Tree based collision search O(N log(N))

REB_COLLISION_MERCURIUS = 3

Direct collision search optimized for MERCURIUS.

REB_COLLISION_LINE = 4

Direct collision search O(N^2), looks for collisions by assuming a linear path over the last timestep.

REB_COLLISION_LINETREE = 5

Tree-based collision search O(N log(N)), looks for collisions by assuming a linear path over the last timestep.

enum [anonymous]

Available integrators.

Values:

REB_INTEGRATOR_IAS15 = 0

IAS15 integrator, 15th order, non-symplectic (default)

REB_INTEGRATOR_WHFAST = 1

WHFast integrator, symplectic, 2nd order, up to 11th order correctors.

REB_INTEGRATOR_SEI = 2

SEI integrator for shearing sheet simulations, symplectic, needs OMEGA variable.

REB_INTEGRATOR_LEAPFROG = 4

LEAPFROG integrator, simple, 2nd order, symplectic.

REB_INTEGRATOR_NONE = 7

Do not integrate anything.

REB_INTEGRATOR_JANUS = 8

Bit-wise reversible JANUS integrator.

REB_INTEGRATOR_MERCURIUS = 9

MERCURIUS integrator.

REB_INTEGRATOR_SABA = 10

SABA integrator family (Laskar and Robutel 2001)

REB_INTEGRATOR_EOS = 11

Embedded Operator Splitting (EOS) integrator family (Rein 2019)

enum [anonymous]

Available boundary conditions.

Values:

REB_BOUNDARY_NONE = 0

Do not check for anything (default)

REB_BOUNDARY_OPEN = 1

Open boundary conditions. Removes particles if they leave the box.

REB_BOUNDARY_PERIODIC = 2

Periodic boundary conditions.

REB_BOUNDARY_SHEAR = 3

Shear periodic boundary conditions, needs OMEGA variable.

enum [anonymous]

Available gravity routines.

Values:

REB_GRAVITY_NONE = 0

Do not calculate graviational forces.

REB_GRAVITY_BASIC = 1

Basic O(N^2) direct summation algorithm, choose this for shearing sheet and periodic boundary conditions.

REB_GRAVITY_COMPENSATED = 2

Direct summation algorithm O(N^2) but with compensated summation, slightly slower than BASIC but more accurate.

REB_GRAVITY_TREE = 3

Use the tree to calculate gravity, O(N log(N)), set opening_angle2 to adjust accuracy.

REB_GRAVITY_MERCURIUS = 4

Special gravity routine only for MERCURIUS.

REB_GRAVITY_JACOBI = 5

Special gravity routine which includes the Jacobi terms for WH integrators.

reb_simulation::[anonymous] visualization

Available collision routines.

reb_simulation::[anonymous] collision

Available collision routines.

reb_simulation::[anonymous] integrator

Available integrators.

reb_simulation::[anonymous] boundary

Available boundary conditions.

reb_simulation::[anonymous] gravity

Available gravity routines.

Integrator structs (the contain integrator specific variables and temporary data structures)

struct reb_simulation_integrator_sei ri_sei

The SEI struct.

struct reb_simulation_integrator_whfast ri_whfast

The WHFast struct.

struct reb_simulation_integrator_saba ri_saba

The SABA struct.

struct reb_simulation_integrator_ias15 ri_ias15

The IAS15 struct.

struct reb_simulation_integrator_mercurius ri_mercurius

The MERCURIUS struct.

struct reb_simulation_integrator_janus ri_janus

The JANUS struct.

struct reb_simulation_integrator_eos ri_eos

The EOS struct.

Callback functions

void (*additional_forces)(struct reb_simulation *const r)

This function allows the user to add additional (non-gravitational) forces.

void (*pre_timestep_modifications)(struct reb_simulation *const r)

This function allows the user to make changes before each timestep.

void (*post_timestep_modifications)(struct reb_simulation *const r)

This function allows the user to make changes after each timestep.

void (*heartbeat)(struct reb_simulation *r)

This function is called at the beginning of the simulation and at the end of each timestep.

void (*display_heartbeat)(struct reb_simulation *r)

This function is called at the beginning of the simulation and at the end of each timestep.

double (*coefficient_of_restitution)(const struct reb_simulation *const r, double v)

Return the coefficient of restitution. By default it is NULL, assuming a coefficient of 1.

The velocity of the collision is given to allow for velocity dependent coefficients of restitution.

int (*collision_resolve)(struct reb_simulation *const r, struct reb_collision)

Resolve collision within this function. By default it is NULL, assuming hard sphere model.

A return value of 0 indicates that both particles remain in the simulation. A return value of 1 (2) indicates that particle 1 (2) should be removed from the simulation. A return value of 3 indicates that both particles should be removed from the simulation.

void (*free_particle_ap)(struct reb_particle *p)

Free particle’s ap pointer. Called in reb_remove function.

void (*extras_cleanup)(struct reb_simulation *r)

Called in reb_free_pointers function for any necessary cleanup in external libraries that depend on the simulation structure.

Hooks for external libraries

void *extras

Pointer to connect additional (optional) libraries, e.g., reboundx.

Main REBOUND functions

group MainRebFunctions

These are the functions that typically need to be called by the user.

Functions

struct reb_simulation *reb_create_simulation(void)

Creates and initialises a REBOUND simulation.

Allocate memory for one reb_simulation structure, initialise all variables and returni the pointer to the reb_simulation sructure. This function must be called before any particles are added.

struct reb_simulation *reb_copy_simulation(struct reb_simulation *r)

Creates a deep copy of a REBOUND simulation.

All simulation data, including all particle data will be copied. Function pointers need to be set manually afterwards. Working on the copy will not affect the original simulation.

int reb_diff_simulations(struct reb_simulation *r1, struct reb_simulation *r2, int output_option)

Compares to REBOUND simulations bit by bit.

Return

If r1 and r2 are exactly equal to each other then 0 is returned, otherwise 1. The walltime parameter in the REBOUND struct is ignored in this comparison.

Parameters
  • r1: The first REBOUND simulation to be compared.

  • r2: The second REBOUND simulation to be compared.

  • output_option: Is set to 1, then the output is printed on the screen. If set to 2, only the return value indicates if the simulations are different.

void reb_init_simulation(struct reb_simulation *r)

Initialize reb_simulation structure.

Same as reb_create_simulation() but does not allocate memory for structure itself.

Parameters
  • r: Structure to be initialized (needs to be allocated externally).

void reb_step(struct reb_simulation *const r)

Performon one integration step.

You rarely want to call this function yourself. Use reb_integrate instead.

Parameters
  • r: The rebound simulation to be integrated by one step.

enum REB_STATUS reb_integrate(struct reb_simulation *const r, double tmax)

Performs the actual integration.

This function performs an integration from the current time t until time tmax.

Return

This function returns an integer, indicating the success of the integration.

Parameters
  • r: The rebound simulation to be integrated.

  • tmax: The time to be integrated to. Set this to INFINITY to integrate forever.

void reb_integrator_synchronize(struct reb_simulation *r)

Synchronize particles manually at end of timestep.

This function should be called if the WHFAST integrator is used, safe_mode is set to zero and an output is needed. This advances the positions and velocities to be synchronized. If enabled, it also applies the symplectic corrector. If safe_mode is enabled, this function has no effect.

Parameters
  • r: The rebound simulation to be synchronized

void reb_integrator_reset(struct reb_simulation *r)

Cleanup all temporarily stored integrator values.

Parameters
  • r: The rebound simulation to be considered

void reb_configure_box(struct reb_simulation *const r, const double boxsize, const int root_nx, const int root_ny, const int root_nz)

Configure the boundary/root box.

This function helps to setup the variables for the simulation box. Call this function when using open, periodic or shear periodic boundary conditions.

Parameters
  • r: The rebound simulation to be considered

  • boxsize: The size of the root box

  • root_nx: The numbe rof root boxes in the x direction.

  • root_ny: The numbe rof root boxes in the y direction.

  • root_nz: The numbe rof root boxes in the z direction.

void reb_free_simulation(struct reb_simulation *const r)

Frees up all space used by a REBOUND simulation and the reb_simulation structure itself.

The REBOUND simulation is not usable anymore after being passed to this function.

Parameters
  • r: The rebound simulation to be freed

void reb_add(struct reb_simulation *const r, struct reb_particle pt)

Adds a particle to the simulation.

This function adds the particle pt to the simulation.

Parameters
  • r: The rebound simulation to which the particle will be added

  • pt: The particle to be added. Note that this is a structure, not a reference to a structure.

void reb_remove_all(struct reb_simulation *const r)

Remove all particles.

Parameters
  • r: The rebound simulation to be considered

int reb_remove(struct reb_simulation *const r, int index, int keepSorted)

Remove a particle by the position in particles array.

Return

Returns 1 if particle was successfully removed, 0 if index passed was out of range.

Parameters
  • r: The rebound simulation to be considered

  • index: The index in the particles array of the particle to be removed.

  • keepSorted: Set to 1, then particles with indices higher than index are all shifted down one position, ensuring the ordering remains.

int reb_remove_by_hash(struct reb_simulation *const r, uint32_t hash, int keepSorted)

Remove a particle by its hash.

see examples/removing_particles_from_simulation.

Return

Returns 1 if particle successfully removed, 0 if hash was not found in the particles array.

Parameters
  • r: The rebound simulation to be considered

  • id: The hash of the particle to be removed.

  • keepSorted: If set to 1 keep the particles with indices in the particles array higher than the one with the passed id are all shifted down one position, ensuring the ordering remains.

struct reb_particle *reb_get_particle_by_hash(struct reb_simulation *const r, uint32_t hash)

Get a pointer to a particle by its hash.

see examples/uniquely_identifying_particles_with_hashes.

Return

A pointer to the particle if found, NULL otherwise.

Parameters
  • r: The rebound simulation to be considered.

  • hash: The hash of the particle to search for.

void reb_run_heartbeat(struct reb_simulation *const r)

Run the heartbeat function and check for escaping/colliding particles.

You rarely want to call this function yourself. It is used internally to call the function you set to the heartbeat variable in reb_simulation.

Parameters
  • r: The rebound simulation to be considered

double reb_integrator_mercurius_L_mercury(const struct reb_simulation *const r, double d, double dcrit)

A force switching function for the MERCURIUS integrator. This function implements the same polynomial switching function used in MERCURY.

double reb_integrator_mercurius_L_infinite(const struct reb_simulation *const r, double d, double dcrit)

A force switching function for the MERCURIUS integrator. This function implements an infinitely differentiable switching function.

int reb_collision_resolve_halt(struct reb_simulation *const r, struct reb_collision c)

Resolve collision by simply halting the integration and setting r->status=REB_EXIT_COLLISION (Default)

int reb_collision_resolve_hardsphere(struct reb_simulation *const r, struct reb_collision c)

Hardsphere collision resolving routine (default).

int reb_collision_resolve_merge(struct reb_simulation *const r, struct reb_collision c)

Merging collision resolving routine.

Merges particle with higher index into particle of lower index. Conserves mass, momentum and volume.

Tool functions

group ToolsRebFunctions

Functions

double reb_random_uniform(double min, double max)

Return uniformly distributed random variable in a given range.

Return

A random variable

Parameters
  • min: Minimum value.

  • max: Maximum value.

double reb_random_powerlaw(double min, double max, double slope)

Returns a random variable drawn form a powerlaw distribution.

Return

A random variable

Parameters
  • min: Minimum value.

  • max: Maximum value.

  • slope: Slope of powerlaw distribution.

double reb_random_normal(double variance)

Return a random number with normal distribution.

Algorithm by D.E. Knut, 1997, The Art of Computer Programmin, Addison-Wesley.

Return

A random variable

Parameters
  • variance: Variance of normal distribution.

double reb_random_rayleigh(double sigma)

Return a random variable drawn form a Rayleigh distribution.

Calculated as described on Rayleigh distribution wikipedia page

Return

A random variable

Parameters
  • sigma: Scale parameter.

void reb_move_to_hel(struct reb_simulation *const r)

Move to the heliocentric frame.

This function moves all particles to the heliocentric frame. Note that the simulation will not stay in the heliocentric frame as it is not an inertial frame. Variational particles are not affected by the function.

Parameters
  • r: The rebound simulation to be considered

void reb_move_to_com(struct reb_simulation *const r)

Move to center of momentum and center of mass frame.

This function moves all particles to the center of mass frame (sometimes also called center of momentum frame). In this frame the center of mass is at rest. It is recommended to call this function before you are doing a long term orbit integration. If the particles are slowly drifting away from the coordinate origin, numerical errors might build up.

Parameters
  • r: The rebound simulation to be considered

struct reb_particle reb_get_com(struct reb_simulation *r)

Returns the center of mass.

Return

The center of mass as a particle (mass, position and velocity correspond to the center of mass)

Parameters
  • r: The rebound simulation to be considered

struct reb_particle reb_get_com_of_pair(struct reb_particle p1, struct reb_particle p2)

Returns the center of mass of two particles.

Return

The center of mass as a particle (mass, position and velocity correspond to the center of mass)

Parameters
  • p1: One of the two particles

  • p2: One of the two particles

void reb_serialize_particle_data(struct reb_simulation *r, uint32_t *hash, double *m, double *radius, double (*xyz)[3], double (*vxvyvz)[3], double (*xyzvxvyvz)[6])

Sets arrays to particle data.

This function can be used to quickly access particle data in a serialized form. NULL pointers will not be set.

Parameters
  • r: The rebound simulation to be considered

  • hash: 1D array to to hold particle hashes

  • mass: 1D array to to hold particle masses

  • radius: 1D array to to hold particle radii

  • xyz: 3D array to to hold particle positions

  • vxvyvz: 3D array to to hold particle velocities

  • xyzvxvyvz: 3D array to to hold particle positions and velocities

void reb_set_serialized_particle_data(struct reb_simulation *r, uint32_t *hash, double *m, double *radius, double (*xyz)[3], double (*vxvyvz)[3], double (*xyzvxvyvz)[6])

Sets particle data to data provided in arrays.

This function can be used to quickly set particle data in a serialized form. NULL pointers will not be accessed.

Parameters
  • r: The rebound simulation to be considered

  • hash: 1D array to to hold particle hashes

  • mass: 1D array to to hold particle masses

  • radius: 1D array to to hold particle radii

  • xyz: 3D array to to hold particle positions

  • vxvyvz: 3D array to to hold particle velocities

  • xyzvxvyvz: 3D array to to hold particle positions and velocities

struct reb_particle reb_get_com_without_particle(struct reb_particle com, struct reb_particle p)

Takes the center of mass of a system of particles and returns the center of mass with one of the particles removed.

Return

The center of mass with particle p removed.

Parameters
  • com: A particle structure that holds the center of mass state for a system of particles (mass, position, velocity).

  • p: The particle to be removed from com.

int reb_get_particle_index(struct reb_particle *p)

Returns a particle pointer’s index in the simulation it’s in.

Return

The integer index of the particle in its simulation (will return -1 if not found in the simulation).

Parameters
  • p: A pointer to the particle

struct reb_particle reb_get_com_range(struct reb_simulation *r, int first, int last)

Returns the center of mass for particles with indices between first (inclusive) and last (exclusive).

For example, reb_get_com_range(r, 6, 9) returns COM for particles 6, 7 and 8.

Return

A reb_particle structure for the center of mass of all particles in range [first, last). Returns particle filled with zeros if passed last <= first.

Parameters
  • r: A pointer to the simulation structure.

  • first: First index in range to consider.

  • last: Will consider particles with indices < last (i.e., particle with index last not considered).

struct reb_particle reb_get_jacobi_com(struct reb_particle *p)

Returns the jacobi center of mass for a given particle.

Return

A reb_particle structure for the center of mass of all particles with lower index. Returns particles[0] if passed the 0th particle.

Parameters
  • p: A pointer to the particle

Output functions

group OutputRebFunctions

List of the built-in output functions for REBOUND

Functions

int reb_output_check(struct reb_simulation *r, double interval)

This function checks if a new output is required at this time.

This is typically used within the heartbeat function to generate equally spaced outputs.

Return

The return value is 1 if an output is required and 0 otherwise.

Parameters
  • interval: Output interval.

  • r: The rebound simulation to be considered

void reb_output_timing(struct reb_simulation *r, const double tmax)

Output status information on the screen.

Outputs the current number of particles, the time and the time difference since the last output to the screen.

Parameters
  • r: The rebound simulation to be considered

  • tmax: The maximum integration time (used to calculate the progress in percent)

void reb_output_orbits(struct reb_simulation *r, char *filename)

Append an ASCII file with orbital paramters of all particles.

The orbital parameters are calculated in Jacobi coordinates. Particles are assumed to be sorted from the inside out, the central object having index 0. Each time the function is called N-1 rows are appended to the file with name filename. Each row in the file corresponds to one particle and contains the following columns (tab separated): time, semi-major axis, eccentricity, inclination, Omega (longitude ascending node), omega (argument of pericenter), lambda (mean longitude), period, f (true anomaly).

Parameters
  • r: The rebound simulation to be considered

  • filename: Output filename.

void reb_output_binary(struct reb_simulation *r, const char *filename)

Save the reb_simualtion structure as a binary.

This function can be used to save the current status of a REBOUND simualtion and later restart the simualtion.

Parameters
  • r: The rebound simulation to be considered

  • filename: Output filename.

void reb_binary_diff(char *buf1, size_t size1, char *buf2, size_t size2, char **bufp, size_t *sizep)

This function compares two REBOUND simulations and records the difference in a buffer.

This is used for taking a SimulationArchive Snapshot.

Parameters
  • buf1: The buffer corresponding to the first rebound simulation to be compared

  • buf2: The buffer corresponding to the second rebound simulation to be compared

  • bufp: The buffer which will contain the differences.

int reb_binary_diff_with_options(char *buf1, size_t size1, char *buf2, size_t size2, char **bufp, size_t *sizep, int output_option)

Same as reb_binary_diff but with more options.

Return

0 is returned if the simulations do not differ (are equal). 1 is return if they differ.

Parameters
  • output_option: If set to 0, the differences are written to bufp. If set to 1, printed on the screen. If set to 2, then only the return value indicates any differences.

void reb_output_ascii(struct reb_simulation *r, char *filename)

Append the positions and velocities of all particles to an ASCII file.

Parameters
  • r: The rebound simulation to be considered

  • filename: Output filename.

void reb_output_binary_positions(struct reb_simulation *r, const char *filename)

Write the positions of all particles to a binary file.

Parameters
  • r: The rebound simulation to be considered

  • filename: Output filename.

void reb_output_velocity_dispersion(struct reb_simulation *r, char *filename)

Append the velocity dispersion of the particles to an ASCII file.

Parameters
  • r: The rebound simulation to be considered

  • filename: Output filename.

Particle setup functions

group SetupRebFunctions

List of the built-in setup helper functions

Derivative functions

This function calculates the first/second derivative of a Keplerian orbit.

Derivatives of Keplerian orbits are required for variational equations, in particular for optimization problems. The derivative is calculated with respect to the variables that appear in the function name. One variable implies that a first derivative is returned, two variables implies that a second derivate is returned. Classical orbital parameters and those introduced by Pal (2009) are supported. Pal coordinates have the advantage of being analytical (i.e. infinite differentiable). Classical orbital parameters may have singularities, for example when e is close to 0. Note that derivatives with respect to Cartesian coordinates are trivial and therefore not implemented as seperate functions. The following variables are supported: a, e, inc, f, omega, Omega, h, k, ix, iy and m (mass).

Return

The derivative as a particle structre. Each structure element is a derivative.

Parameters
  • G: The gravitational constant

  • primary: The primary of the Keplerian orbit

  • po: The original partical for which the derivative is to be calculated.

struct reb_particle reb_derivatives_lambda(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_h(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_k(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_k_k(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_h_h(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_lambda_lambda(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_k_lambda(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_h_lambda(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_k_h(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_a(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_a_a(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_ix(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_ix_ix(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_iy(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_iy_iy(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_k_ix(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_h_ix(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_lambda_ix(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_lambda_iy(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_h_iy(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_k_iy(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_ix_iy(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_a_ix(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_a_iy(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_a_lambda(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_a_h(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_a_k(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_m(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_m_a(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_m_lambda(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_m_h(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_m_k(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_m_ix(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_m_iy(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_m_m(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_e(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_e_e(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_inc(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_inc_inc(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_Omega(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_Omega_Omega(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_omega(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_omega_omega(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_f(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_f_f(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_a_e(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_a_inc(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_a_Omega(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_a_omega(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_a_f(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_e_inc(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_e_Omega(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_e_omega(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_e_f(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_m_e(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_inc_Omega(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_inc_omega(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_inc_f(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_m_inc(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_omega_Omega(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_Omega_f(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_m_Omega(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_omega_f(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_m_omega(double G, struct reb_particle primary, struct reb_particle po)
struct reb_particle reb_derivatives_m_f(double G, struct reb_particle primary, struct reb_particle po)

Enums

enum reb_input_binary_messages

Enum describing possible errors that might occur during binary file reading.

Values:

REB_INPUT_BINARY_WARNING_NONE = 0
REB_INPUT_BINARY_ERROR_NOFILE = 1
REB_INPUT_BINARY_WARNING_VERSION = 2
REB_INPUT_BINARY_WARNING_POINTERS = 4
REB_INPUT_BINARY_WARNING_PARTICLES = 8
REB_INPUT_BINARY_ERROR_FILENOTOPEN = 16
REB_INPUT_BINARY_ERROR_OUTOFRANGE = 32
REB_INPUT_BINARY_ERROR_SEEK = 64
REB_INPUT_BINARY_WARNING_FIELD_UNKOWN = 128
REB_INPUT_BINARY_ERROR_INTEGRATOR = 256

Functions

double reb_tools_M_to_f(double e, double M)

returns the true anomaly for a given eccentricity and mean anomaly

Return

True anomaly

Parameters
  • e: Eccentricity

  • M: Mean anomaly

struct reb_particle reb_tools_orbit2d_to_particle(double G, struct reb_particle primary, double m, double a, double e, double omega, double f)

Initialize a particle on an orbit in the xy plane.

Return

Returns a particle structure with the given orbital parameters.

Parameters
  • G: Gravitational constant.

  • primary: Particle structure for the orbit’s reference body.

  • m: Mass of the particle.

  • a: Semi-major axis of the particle.

  • e: Eccentricity of the particle.

  • omega: Pericenter of the particle.

  • f: true anomaly of the particle.

struct reb_particle reb_tools_orbit_to_particle_err(double G, struct reb_particle primary, double m, double a, double e, double i, double Omega, double omega, double f, int *err)

Initialize a particle on a 3D orbit, passing an error variable to flag why particle is set to nan. See Fig. 2.13 of Murray & Dermott Solar System Dynamics for diagram.

Error codes:

  1. Can’t set e exactly to 1.

  2. Eccentricity can never be less than zero.

  3. Bound orbit (a>0) can’t have e>1.

  4. Unbound orbit (a<0) can’t have e<1.

  5. Unbound orbit can’t have f set beyond the asymptotes defining the particle.

    Return

    Returns a particle structure with the given orbital parameters.

    Parameters
    • G: Gravitational constant.

    • primary: Particle structure for the orbit’s reference body.

    • m: Mass of the particle.

    • a: Semi-major axis of the particle.

    • e: Eccentricity of the particle.

    • i: inclination of the particle to the reference plane..

    • Omega: Longitude of the ascending node of the particle.

    • omega: argument of pericenter of the particle.

    • f: true anomaly of the particle.

    • err: Pointer to error code that wil be set by this function. Used for checking why particle was set to nans.

struct reb_particle reb_tools_orbit_to_particle(double G, struct reb_particle primary, double m, double a, double e, double i, double Omega, double omega, double f)

Initialize a particle on a 3D orbit. See Fig. 2.13 of Murray & Dermott Solar System Dynamics for diagram.

Return

Returns a particle structure with the given orbital parameters.

Parameters
  • G: Gravitational constant.

  • primary: Particle structure for the orbit’s reference body.

  • m: Mass of the particle.

  • a: Semi-major axis of the particle.

  • e: Eccentricity of the particle.

  • i: inclination of the particle to the reference plane.

  • Omega: Longitude of the ascending node of the particle.

  • omega: argument of pericenter of the particle.

  • f: true anomaly of the particle.

struct reb_orbit reb_tools_particle_to_orbit_err(double G, struct reb_particle p, struct reb_particle primary, int *err)

This function calculates orbital elements for a given particle, passing an error variable to flag why orbit is set to nan.

Error codes:

  1. Primary has no mass.

  2. Particle and primary positions are the same.

    Return

    reb_orbit struct with orbital parameters.

    Parameters
    • G: The gravitational constant.

    • p: reb_particle for which the orbit is calculated.

    • primary: Particle structure for the orbit’s reference body.

    • err: error code for checking why orbit was set to nans.

struct reb_orbit reb_tools_particle_to_orbit(double G, struct reb_particle p, struct reb_particle primary)

This function calculates orbital elements for a given particle.

Return

reb_orbit struct with orbital parameters.

Parameters
  • G: The gravitational constant.

  • p: reb_particle for which the orbit is calculated.

  • primary: Particle structure for the orbit’s reference body.

struct reb_particle reb_tools_pal_to_particle(double G, struct reb_particle primary, double m, double a, double lambda, double k, double h, double ix, double iy)

Initialize a particle on a 3D orbit. See Pal 2009 for a definition of these coordinates.

Pal describes a coordinate system for Keplerian Orbits that is analytical (i.e. infinitely differentiable) between spatial coordinates and orbital elements. See http://adsabs.harvard.edu/abs/2009MNRAS.396.1737P

Return

Returns a particle structure with the given orbital parameters.

Parameters
  • G: Gravitational constant.

  • primary: Particle structure for the orbit’s reference body.

  • m: Mass of the particle.

  • a: Semi-major axis of the particle.

  • lambda: longitude.

  • k: Eccentricity/pericenter k = e*cos(w).

  • h: Eccentricity/pericenter h = e*sin(w).

  • ix: Inclination, x component.

  • iy: Inclination, y component.

struct reb_simulation *reb_create_simulation_from_binary(char *filename)

Reads a binary file.

Also initialises the particles array with data form the binary file. This can be used to restart a simualtion.

Return

Returns a pointer to a REBOUND simulation.

Parameters
  • filename: Filename to be read.

void reb_tools_init_plummer(struct reb_simulation *r, int _N, double M, double R)

This function sets up a Plummer sphere.

Parameters
  • r: The rebound simulation to be considered

  • _N: Number of particles in the plummer sphere.

  • M: Total mass of the cluster.

  • R: Characteristic radius of the cluster.

char *reb_read_char(int argc, char **argv, const char *argument)

Reads arguments from the command line.

Return

Returns NULL if argument was not given. Return the argument otherwise.

Parameters
  • argc: Number of command line arguments.

  • argv: Array of command line arguments.

  • argument: Argument to look for.

double reb_read_double(int argc, char **argv, const char *argument, double _default)

Reads arguments as a double value from the command line.

Return

Returns _default if argument was not given. Return the argument converted to double otherwise.

Parameters
  • argc: Number of command line arguments.

  • argv: Array of command line arguments.

  • argument: Argument to look for.

  • _default: Default value.

int reb_read_int(int argc, char **argv, const char *argument, int _default)

Reads arguments as a int value from the command line.

Return

Returns _default if argument was not given. Return the argument converted to int otherwise.

Parameters
  • argc: Number of command line arguments.

  • argv: Array of command line arguments.

  • argument: Argument to look for.

  • _default: Default value.

Miscellaneous functions

group MiscRebFunctions

List of the miscellaneous helper functions for REBOUND

Functions

double reb_tools_energy(const struct reb_simulation *const r)

Calculate the total energy (potential and kinetic).

Does not work for SEI (shearing sheet simulations).

Calculate the total energy (potential and kinetic).

Return

Total energy.

Parameters
  • r: The rebound simulation to be considered.

struct reb_vec3d reb_tools_angular_momentum(const struct reb_simulation *const r)

Calculate the system’s angular momentum.

Return

The angular momentum vector as a reb_vec3d struct.

Parameters
  • r: The rebound simulation to be considered.

int reb_add_var_1st_order(struct reb_simulation *const r, int testparticle)

Add and initialize a set of first order variational particles.

Return

Returns the index of the first variational particle added

Parameters
  • r: The rebound simulation to be considered

  • testparticle: This flag determines if the set of variational particles is for a testparticle or not. If testparticle is >= 0, then only one variational particle (the test particle) will be added. If testparticle is -1, one variational particle for each real particle will be added.

int reb_add_var_2nd_order(struct reb_simulation *const r, int testparticle, int index_1st_order_a, int index_1st_order_b)

Add and initialize a set of second order variational particles.

Note that a set of second order variational particles requires two sets of first order variational equations.

Return

Returns the index of the first variational particle added

Parameters
  • r: The rebound simulation to be considered

  • testparticle: This flag determines if the set of variational particles is for a testparticle or not. If testparticle is >= 0, then only one variational particle (the test particle) will be added. If testparticle is -1, one variational particle for each real particle will be added.

  • index_1st_order_a: The index of the corresponding first variational particles.

  • index_1st_order_b: The index of the corresponding first variational particles.

void reb_tools_megno_init(struct reb_simulation *const r)

Init the MEGNO particles, enable MEGNO calculation.

Parameters
  • r: The rebound simulation to be considered

void reb_tools_megno_init_seed(struct reb_simulation *const r, unsigned int seed)

Init the MEGNO particles, enable MEGNO calculation, and specify a seed for the random number generation.

Parameters
  • r: The rebound simulation to be considered

  • seed: The seed to use for the random number generator

double reb_tools_calculate_megno(struct reb_simulation *r)

Get the current MEGNO value.

Return

Returns the current value of the MEGNO

Parameters
  • r: The rebound simulation to be considered

double reb_tools_calculate_lyapunov(struct reb_simulation *r)

Returns the largest Lyapunov characteristic number (LCN), or maximal Lyapunov exponent.

MEGNO needs to be enabled to calculate this value.

Return

Returns the current CN

Parameters
  • r: The rebound simulation to be considered

uint32_t reb_hash(const char *str)

Returns hash for passed string.

Return

hash for the passed string.

Parameters
  • str: String key.

struct reb_particle reb_particle_nan(void)

Returns a reb_particle structure with fields/hash/ptrs initialized to nan/0/NULL.

Return

reb_particle with fields initialized to nan.

void reb_exit(const char *const msg)

Print out an error message, then exit in a semi-nice way.

void reb_warning(struct reb_simulation *const r, const char *const msg)

Print or store a warning message, then continue.

void reb_error(struct reb_simulation *const r, const char *const msg)

Print or store an error message, then continue.

int reb_get_next_message(struct reb_simulation *const r, char *const buf)

Get the next warning message stored. Used only if save_messages==1.

Return

Return value is 0 if no messages are present, 1 otherwise.

Parameters
  • r: The rebound simulation to be considered

  • buf: The buffer in which the error message it copied (needs to be at least reb_max_messages_length long).