hannorein / rebound

💫 An open-source multi-purpose N-body code.

Welcome to REBOUND

REBOUND is an N-body integrator, i.e. a software package that can integrate the motion of particles under the influence of gravity. The particles can represent stars, planets, moons, ring or dust particles. REBOUND is very flexible and can be customized to accurately and efficiently solve many problems in astrophysics.

REBOUND 24 Meeting

Join us for the first virtual 2-day REBOUND meeting, aimed to bring REBOUND developers and users together. 1 day of science talks, 1 day of technical/hands-on talks, discussions on REBOUND's future + social events!

Abstract submission and registration are now open at https://hannorein.github.io/rebound24/.

Features

• No dependencies on external libraries.
• Runs natively on Linux, MacOS, and Windows.
• Symplectic integrators WHFast, SEI, LEAPFROG, EOS.
• Hybrid symplectic integrators for planetary dynamics with close encounters MERCURIUS
• High order symplectic integrators for integrating planetary systems SABA, WH Kernel methods.
• High accuracy non-symplectic integrator with adaptive time-stepping IAS15.
• Can integrate arbitrary user-defined ODEs that are coupled to N-body dynamics for tides, spin, etc
• Support for collisional/granular dynamics, various collision detection routines
• The computationally intensive parts of the code are written entirely in C, conforming to the ISO standard C99, and can be used as a thread-safe shared library
• Easy-to-use Python module, installation in 3 words: pip install rebound
• Real-time, 3D visualization, for both C and Python.
• Extensive set of example problems for both C and Python. You can run examples directly from your browser without the need to download or install anything.
• Parallelized WHFast512 integrator for super fast integrations of planetary systems with SIMD AVX512 instructions
• Parallelized with OpenMP (for shared memory systems)
• Parallelized with MPI is supported for some special use cases only (using an essential tree for gravity and collisions)
• The code is 100% open-source. All features are included in the public repository on github.

Try out REBOUND

You can try out REBOUND without installing it. Simply head over to readthedocs.org. All the C examples have been compiled with emscripten and can run directly in your browser.

One minute installation

You can install REBOUND with pip if you want to only use the python version of REBOUND:

pip install rebound

Then, you can run a simple REBOUND simulation such as

import rebound
sim = rebound.Simulation()
sim.add(m=1.0)
sim.add(m=1.0e-3, a=1.0)
sim.integrate(1000.)
sim.status()

If you want to use the C version of REBOUND simply copy and paste this line into your terminal (it won't do anything bad, we promise):

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

Documentation

The full documentation with many examples, changelogs and tutorials can be found at

https://rebound.readthedocs.org

If you have trouble installing or using REBOUND, please open an issue on github and we'll try to help as much as we can.

There are also short YouTube videos describing various aspects of REBOUND available at https://www.youtube.com/channel/UCNmrCzxcmWVTBwtDPPLxkkw .

Related projects

Additional physics

To easily incorporate additional physics modules such as migration forces, GR effects and spin into your REBOUND simulations, see REBOUNDx at https://github.com/dtamayo/reboundx

Analytical and semianalytical tools

If you're interested in comparing numerical simulations to analytical and semianalytical tools for celestial mechanics, see Celmech at https://github.com/shadden/celmech

Ephemeris-quality integrations of test particles

To generate ephemeris-quality integrations of test particles in the Solar System with a precision on par with JPL's small body integrator, see ASSIST at https://github.com/matthewholman/assist

Papers

There are several papers describing the functionality of REBOUND.

1. Rein & Liu 2012 (Astronomy and Astrophysics, Volume 537, A128) describes the code structure and the main feature including the gravity and collision routines for many particle systems. http://adsabs.harvard.edu/abs/2012A%26A...537A.128R

2. Rein & Tremaine 2011 (Monthly Notices of the Royal Astronomical Society, Volume 415, Issue 4, pp. 3168-3176) describes the Symplectic Epicycle integrator for shearing sheet simulations. https://ui.adsabs.harvard.edu/abs/2011MNRAS.415.3168R

3. Rein & Spiegel 2015 (Monthly Notices of the Royal Astronomical Society, Volume 446, Issue 2, p.1424-1437) describes the versatile high order integrator IAS15 which is now part of REBOUND. http://adsabs.harvard.edu/abs/2015MNRAS.446.1424R

4. Rein & Tamayo 2015 (Monthly Notices of the Royal Astronomical Society, Volume 452, Issue 1, p.376-388) describes WHFast, the fast and unbiased implementation of a symplectic Wisdom-Holman integrator for long term gravitational simulations. http://adsabs.harvard.edu/abs/2015MNRAS.452..376R

5. Rein & Tamayo 2016 (Monthly Notices of the Royal Astronomical Society, Volume 459, Issue 3, p.2275-2285) develop the framework for second order variational equations. https://ui.adsabs.harvard.edu/abs/2016MNRAS.459.2275R

6. Rein & Tamayo 2017 (Monthly Notices of the Royal Astronomical Society, Volume 467, Issue 2, p.2377-2383) describes the Simulationarchive for exact reproducibility of N-body simulations. https://ui.adsabs.harvard.edu/abs/2017MNRAS.467.2377R

7. Rein & Tamayo 2018 (Monthly Notices of the Royal Astronomical Society, Volume 473, Issue 3, p.3351–3357) describes the integer based JANUS integrator. https://ui.adsabs.harvard.edu/abs/2018MNRAS.473.3351R

8. Rein, Hernandez, Tamayo, Brown, Eckels, Holmes, Lau, Leblanc & Silburt 2019 (Monthly Notices of the Royal Astronomical Society, Volume 485, Issue 4, p.5490-5497) describes the hybrid symplectic integrator MERCURIUS. https://ui.adsabs.harvard.edu/abs/2019MNRAS.485.5490R

9. Rein, Tamayo & Brown 2019 (Monthly Notices of the Royal Astronomical Society, Volume 489, Issue 4, November 2019, Pages 4632-4640) describes the implementation of the high order symplectic integrators SABA, SABAC, SABACL, WHCKL, WHCKM, and WHCKC. https://ui.adsabs.harvard.edu/abs/

Acknowledgments

If you use this code or parts of this code for results presented in a scientific publication, we would greatly appreciate a citation. please cite REBOUND. The simplest way to find the citations relevant to the specific setup of your REBOUND simulation is:

sim = rebound.Simulation()
-your setup-
sim.cite()

Contributors

• Hanno Rein, University of Toronto, [email protected]
• Dan Tamayo, Harvey Mudd College, [email protected]
• David S. Spiegel, Institute for Advanced Study Princeton, [email protected]
• Garett Brown, University of Toronto, [email protected]
• Shangfei Liu, Kavli Institute for Astronomy and Astrophysics at Peking University, [email protected]
• Ari Silburt, Penn State University, [email protected]
• and many others! Check the git history to find out who contributed to the code.

REBOUND is open source and you are invited to contribute to this project!

License

REBOUND is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

REBOUND is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with REBOUND. If not, see http://www.gnu.org/licenses/.

rebound's Issues

Bulirsch-Stoer and particle ID/name

(1) I tried to calculate comet orbit (including cometary outgassing and GR-correction as additional force) using RADAU15, but didnt converged (close encounter with the Sun).
In mercury6 by John E. Chambers, it is said that radau15 usually reliable, except for very close encounters or very eccentric (e.g. Sun grazing) orbits.

so, should we make BS integrator? or rebound design is not good for BS-integrator because in BS we have to update acceleration many times..?

(2) Sometimes we need to save detail of encounter/collision or ejected particles in simulation, but in rebound particles are just removed with no ID or name (last particle in "particles" is move to the deleted particle place), what do you think about this?

thank you in advance,

Finish/polish python doc strings

Left to do:

• Units section in simulation.py
• Init routines in particle.py
• Units module units.py

Implement MPI domain decomposition and parallelization

Fix boundary issue for sweep.

There seems to be an issue with collision detection and the sweep method. See overstability example.

WHFast and test particles

We need to think about how/if we can make WHFast work (better) with test particles. Currently it fails a couple of tests that IAS15 passes.

J2000 reference frame

There are two options for the coordinate system in Horizons:

• Earth mean equator and equinox of reference epoch ("frame"):
• Ecliptic and mean equinox of reference epoch ("ecliptic")

Currently we are using "ecliptic". Would "frame" make more sense as the dafault?

See: ftp://ssd.jpl.nasa.gov/pub/ssd/Horizons_doc.pdf

Old swift based WH integrator not conserving energy

Can't get the energy conservation to work for the old swift based WH integrator. Maybe not worth fixing. WHFast should be superior in all possible ways. Unit test is in test_integrator.py.

integrate()

I hate to bring this up again, but I think there might still be a bug for some cases where the integrate function does not integrate up to where it should be. This might only affect borderline cases such as

t = 0
tmax = t + dt

There also seems to be an issue with

t = 1
dt = 0.01
tmax = 0

Add a core utility reb_tools_am()?

to calculate total angular momentum of the system?

list freeglut3-dev as the package to install for Debian (and Ubuntu)

Hi Hanno,

Inspired by your talk today, I cloned rebound. All I needed to compile it was to install Debian package freeglut3-dev. Likely, the same would hold for Ubuntu. I thought you might want to note that in the readme.

All the best,

Marten

Inefficient implementation of "get_min_ratio" in integrator_hybrid.c

The current implementation of the double for loop in get_min_ratio is not optimal, with computations that are done identically in the inner loop and could be performed in the outer loop, saving quite some time. Cheking it on the given example in solar_system with REB_INTEGRATOR_HYBRID, I get an improvement of 9% of the overall time.

Current implementation is:

for (int i=1; i<_N_active; i++){
    struct reb_particle pi = particles[i];
for (int j=1; j<_N_real; j++){
    if (i==j) continue;
    const double dxj = p0.x - particles[j].x;
    const double dyj = p0.y - particles[j].y;
    const double dzj = p0.z - particles[j].z;
    const double r0j2 = dxj*dxj + dyj*dyj + dzj*dzj;

    const double dx = pi.x - particles[j].x;
    const double dy = pi.y - particles[j].y;
    const double dz = pi.z - particles[j].z;
    const double rij2 = dx*dx + dy*dy + dz*dz;

    const double F0j = p0.m/r0j2;
    const double Fij = pi.m/rij2;

    const double ratio = F0j/Fij;

    if (ratio<min_ratio){
        min_ratio = ratio;
    }
}
}

My proposed implementation:

for (int i=1; i<_N_active; i++){
    struct reb_particle pi = particles[i];
    const double dxi = p0.x - pi.x;
    const double dyi = p0.y - pi.y;
    const double dzi = p0.z - pi.z;
    const double r0i2 = dxi*dxi + dyi*dyi + dzi*dzi;

    const double F0i = p0.m/r0i2;
for (int j=1; j<_N_real; j++){
    if (i==j) continue;

    const double dx = pi.x - particles[j].x;
    const double dy = pi.y - particles[j].y;
    const double dz = pi.z - particles[j].z;
    const double rij2 = dx*dx + dy*dy + dz*dz;

    const double Fij = pi.m/rij2;

    const double ratio = F0i/Fij;

    if (ratio<min_ratio){
        min_ratio = ratio;
    }
}
}

It's mathematically identical, but much faster. Since this is not my code, and I'm not familiar with GITHUB, I'd rather let you do the update of the code.

Regards, Jean-Marc Petit

Add softening to WHM

Documentation

I'm in the process of adding documentation to the code. I plan to use the following documentation of tools for automatically generated html pages from doc strings in the source code:

• doxygen
  • For documenting the C part, the shared library librebound.
  • This creates xml output
• sphinx-dox
  • For documenting the python module
  • This creates HTML output
• breathe
  • Allows sphinx to include the doxygen xml output, thus all documentation can be in one place and style.

reb_collision_search() double counting collisions?

It appears to me as I go through the source code, and when I run my code, that reb_collision_search() in the case of REB_COLLISION_DIRECT counts each individual collision twice, and thus tries to resolve it twice.

In collision.c I find:
line 67: loop over all the particles (indexed by i, and i runs over all particles)
line 79: loop over all the particles again (indexed by j, and j runs over all particles)

It seems to me that if a particle given by i=2 is found to collide with another particle given by j=4, that when i=4 a collision will be found when j=2 as the loop is currently set up. In this way, two collisions are found, but both collisions are between the same two particles.

I wanted to make sure my understanding of the code as it is now set up is correct.

A suggested edit is to change line 79 from:
for (int j=0;j<N;j++){
to
for (int j=i+1;j<N;j++){

then no collisions will get double counted (and the statement on line 81, to prevent a particle from colliding with itself, would no longer be necessary).

Units module updates all instances

If one updates the units on one simulation objects, it gets updates on all simulation objects.

>>> import rebound
>>> sim1 = rebound.Simulation()
>>> sim2 = rebound.Simulation()
>>> print sim1.units, sim2.units
{'length': None, 'mass': None, 'time': None} {'length': None, 'mass': None, 'time': None}
>>> sim1.add("Earth")
Searching NASA Horizons for 'Earth'... Found: Earth-Moon Barycenter (3).
>>> print sim1.units, sim2.units
{'length': 'au', 'mass': 'msun', 'time': 'yr2pi'} {'length': 'au', 'mass': 'msun', 'time': 'yr2pi'}

t+=dt

Make sure that t+=dt is at the right position. This might be problematic as check_boundaries determines the box position using t.

Installing Rebound

Good Morning,

I'm having trouble installing rebound. Python is what I use and have installed, how can I install rebound on Python? The website guideline is to copy a command in my terminal, but that gives me error and doesn't install rebound. Could you help me with that?

Thank you so much,
NiDi

CL_INVALID_WORK_GROUP_SIZE (error -34) using gravity_opencl.c

Hello

When I try to run the provided examples/opencl problem on our graphic hardware, the program always crash during the first thirty seconds of simulated time with the following error message:

N_tot= 4096 t= 21.540000 cpu= 0.077756 [s]
Line No: 221 clEnqueueNDRangeKernel failed.
Error= -54

This errors usually indicates that the choosen work group (local) size is either invalid, or not a divider of the global work size. However, after adding some printf debugging to gravity_opencl.c, these size are always fixed to 32 and 4096 for all iterations.

Is this a known problem? If I try to indicate an automatic local size, by replacing &local_size in the source with NULL, the computation does not crash anymore, but runs much more slowly (about 1sec/timestep). Changing the work group size to 64 or 128 does not prevent the crash from ocurring, either.

PS : I modified the makefile a bit to enable compilation on Linux systems (export LIB=-lOpenCL). At the moment, this is only a hack, but I would be interested in submitting patches for Linux and Xeon Phi support, since I have the opportunity to test on these platforms.

Documentation of Simulation

Add more documentation to the main Simulation class in python.

Calculate quadrupole contribution to the gravity.

Spaces vs tabs

I would like to move to a consistent style where we use either tabs or spaces. I am tempted towards spaces. Four spaces for each tab, to be specific. Any strong feelings about this?

dt_last_success

dt_last_success may not be initialized which causes a problem if the first timestep fails

Please consider creating a Debian/Ubuntu package

Dear Hanno,

I am a member of the Debian Astronomy Team. Our aim is to package astronomical software for Debian, from where it will also migrate to derivative distributions like Ubuntu or Mint.

Would you consider packaging rebound for Debian? This would ease the installation of the software and help people keeping it updated. So, it can broaden the user base of it.

If you would package it, we would be glad to help you when you run into problems. Once the package is ready, we could upload it to the Debian archive.

We have a mailing list, which I would recommend to subscribe, since this is our main information exchange. Also, we will have a Debian Sprint in Heidelberg on August 14th 2015 where you are welcome.
If you have any questions, please contact either me directly, or our mailing list.

Memory leak in re_create_simulation_from_binary (src/input.c:116)

I've checked the source code using cppcheck and there is potential memory leak in re_create_simulation_from_binary function.

In line 116 (if the binary file cannot be open) function returns NULL but the memory pointed to by r is not freed.

Implementation of gravity_tree.c

The tree should be used to calculate the force on every particle.

Look at using Binders for REBOUND

http://mybinder.org

Thanks to @minrk for the pointer!

FFT

Write an FFT(/tree hybrid) solver for gravity.

Name

We need to come up with a simple, unique and catchy name.

OpenMP and compensated summation.

Fix bug involving OpenMP and compensated summation. Currently multiple threads might update same particle.

can't write data if usleep enabled

When I set the usleep value of my simulation to 1 or more, the reb_output_orbit function (called in rebound.c) doesn't create the file on which I'm trying to write data. But when I set usleep to a negative value (no OpenGL) it works fine.

Make use of the restrict keyword

We should make sure that the compiler can optimize the particle structure efficiently as it is used all over the place. One way to do this is with the restrict keyword in C99.

See e.g. http://en.wikipedia.org/wiki/Restrict

Exact Finish Time

Currently, the exact finish time variable is inconsistent between python and c. Once again. Certainly my fault.

The question is how to choose the best default behaviour. Right now c has it set to 0 and python to 1.

Idea?

Makefile error, C version

I ran the code

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

However, when it reaches the makefile, I get the following error:

make -C ../../src/
make[1]: Entering directory '/home/pi/rebound/src'
Compiling source file rebound.c ...
cc -c -std=c99 -Wpointer-arith -D_GNU_SOURCE -O3 -march=native -Wall -g -Wno-unknown-pragmas -fPIC -DLIBREBOUND -DOPENGL -o rebound.o rebound.c
*** Error in `cc': double free or corruption (top): 0x0086cfd8 ***
Makefile:12: recipe for target 'rebound.o' failed
make[1]: *** [rebound.o] Aborted
make[1]: Leaving directory '/home/pi/rebound/src'
Makefile:12: recipe for target 'librebound' failed
make: *** [librebound] Error 2

Running just sudo ./make gives this error as well. I am attempting this install on a Raspberry Pi 3.

Fixing memory leak and complexity

There was minor memory leak in the reb_free_simulation() function that actually did not free the simulation but only pointer therein.

I also looked at the complexity of the python Simulation class. I think this was overly complicated. Now Simulation is just a subclass of ctypes.Structure and one can access field variable directly without a setter/getter.

This raises two other issues, however:

1. If we want a setter on a field, e.g. for simualtion.integrator, then we need to make the fields variable private, e.g. by calling it _integrator.
2. The documentation for the getter/setters is gone; naturally, as the getter/setters are gone.

Development branch is on https://github.com/hannorein/rebound/tree/simsim. I'm not totally happy with it yet and there might be an example that does not work yet (although almost all changes should be internal and not be visible from the outside).

Globals free version

I'm working on a thread-safe and global variable free version of the code. The current development is on the branch noglobals2.

compiling -- integrator routines

I am finding some undefined symbols during the link, having to do with how the integrator routines are chosen/defined. In current version, I can't find anywhere where integrator_part1() or other choices are set after the integrator is chosen. I can work around this by adding defs, but probably there is some little piece of code missing from the current zip?

Here is a piece of the compile output
Undefined symbols for architecture x86_64:
"_integrator_part1", referenced from:
_iterate in main.o
"_integrator_part2", referenced from:
_iterate in main.o
"_integrator_update_acceleration", referenced from:
_integrator_ias15_step in integrator_ias15.o
_integrator_whfast_corrector_Z in integrator_whfast.o

Fix formatting issue of docstrings

Some doc strings are not rendered nicely. See e.g.:

Orbit conversion

One thing that's a little confusing right now is that the Orbit class has a different convention for the attributes compared to the argument that the Particle class takes in the constructor.

Orbit class has:

• mean longitude
• true anomaly

Whereas Particle class takes:

• mean anomaly
• true anomaly

I also think the true anomalies in those two classes might be defined differently (at least for circular orbits).

Fix high e variational equation

There might be a floating point issue for large timesteps and high eccentricity test particle orbits in the variational equations.

Warnings of Intel compiler

Several warning messages like:
remark #1572: floating-point equality and inequality comparisons are unreliable.
remark #1419: external declaration in primary source file.
...

Finish MPI parallelization of collisions_sweep.c

Wrong test for switching between SYMPLECTIC and HIGHORDER in integrator_hybrid.c

Although the initialisation of reb_simulation.ri_hybrid.switch_ratio claims it's expressed in Hill's radius (rebound.c), the actual use of it shows that the test is done on comparing the central body force to the force between 2 approaching bodies. This is a total misconception. Since the integrator uses a propagator that follows the Keplerian motion (use of Stumpf functions), one needs to change to HIGHORDER when the perturbation becomes large with respect to the normal Keplerian motion, i.e. when the force between the 2 bodies is comparable to the "gradient" of the central force times the distance between the bodies. That's exactly what the Hill radius is, the distance where the 2 effects are of same amplitude. For the Earth-Moon system for example, one can definitely say the gravitation beteween Earth and Moon is the dominant effect in the local frame, but it turns out that the force acting on the Moon from the Sun it twice as large as the force acting on the Moon from the Earth.

So instead of comparing the ratio of the forces, one should rather check the distance compared to the Hill radius. The default value of ri_hybrid.switch_ratio should be ~5, and get_min_ratio should return the minimum ratio of the distances to the sum of the Hill's radii of the 2 particles.

So I propose the following implementation.

get_min_ratio in integrator_hybrid.c:

    for (int i=1; i<_N_active; i++){
        const double dxi = p0.x - particles[i].x;
        const double dyi = p0.y - particles[i].y;
        const double dzi = p0.z - particles[i].z;
        const double r0i2 = dxi*dxi + dyi*dyi + dzi*dzi;

        particles[i].rh = r0i2*pow((particles[i].m/(p0.m*3.)), 2./3.);
    }
    for (int i=1; i<_N_active; i++){
        struct reb_particle pi = particles[i];
    for (int j=1; j<_N_real; j++){
        if (i==j) continue;

        const double dx = pi.x - particles[j].x;
        const double dy = pi.y - particles[j].y;
        const double dz = pi.z - particles[j].z;
        const double rij2 = dx*dx + dy*dy + dz*dz;

        const double ratio = rij2/(pi.rh+particles[j].rh);

        if (ratio<min_ratio){
            min_ratio = ratio;
        }
    }
    }

This requires adding a field to the reb_particle structure, to hold the squared Hill radius in rebound.h:

/**
 * @brief Structure representing one REBOUND particle.
 * @details This structure is used to represent one particle. 
 * If this structure is changed, the corresponding python structure
 * needs to be changes as well. Also update the equivalent declaration 
 * for MPI in communications_mpi.c.
 */
struct reb_particle {
    double x;           ///< x-position of the particle. 
    double y;           ///< y-position of the particle. 
    double z;           ///< z-position of the particle. 
    double vx;          ///< x-velocity of the particle. 
    double vy;          ///< y-velocity of the particle. 
    double vz;          ///< z-velocity of the particle. 
    double ax;          ///< x-acceleration of the particle. 
    double ay;          ///< y-acceleration of the particle. 
    double az;          ///< z-acceleration of the particle. 
    double m;           ///< Mass of the particle. 
    double r;           ///< Radius of the particle. 
    double rh;          ///< Hill radius of the particle, squared. 
    double lastcollision;       ///< Last time the particle had a physical collision.
    struct reb_treecell* c;     ///< Pointer to the cell the particle is currently in.
    int id;             ///< Unique id to identify particle.
};

And finally, change the default value of ri_hybrid.switch_ratio to 5 in reb_create_simulation in rebound.c:

r->ri_hybrid.switch_ratio = 5; // 5 Hill radii  

As I said in my other issue, since this is not my code, and I'm not familiar with GITHUB, I'd rather let you do the update of the code.

Regards, Jean-Marc Petit

Python - segfault on adding particles with gravity = "tree"

The following minimal example

import rebound

sim = rebound.Simulation()

sim.gravity = "tree"

print "adding particle #1"
sim.add(m=1)
print "adding particle #2"
sim.add(m=1)
print "added 2 particles"

prints

adding particle #1
adding particle #2
Segmentation fault (core dumped)

The issue only presents itself when using the "tree" module and the segfault only happens when adding the second particle. Adding/changing parameters to add() such as particle position or velocity has no effect on the issue.

This is on Ubuntu 15.10 and just installed through pip install rebound.

License

@jobovy brought up the issue that some people might favour a non-GPL license for REBOUND.

A non-GPL license would, for example, allow people to include the integrators in REBOUND in other open source projects that are not under a GPL license (astropy, etc). This seems like a desirable thing to do, I think.

The only problem I have is that I would like to avoid companies using our software for commercial purposes without contributing to the open source community. As an example, I'm thinking of http://www.spaceroutes.com, a company optimizing spacecraft trajectories. They have 14 patents on this topic which I find ethically problematic, to say the least.

I'll leave this issue open for a while. Please post your comments.

OPENMP going crazy in Version 2

Combining the openmp example and the solar_system example, I get strange or bad results.

1. Using REB_INTEGRATOR_HYBRID, the rest being default values, running with 4 threads, I get a message that the integrator goes into HIGHORDER, meaning there are "close encounters", while in the non-OPENMP case, there is no such thing. If running with only 1 thread, things go well, although a little slower than without OPENMP.

2. Could this be due to the integrator ? I decided to use REB_INTEGRATOR_LEAPFORG which is the one used in the OPENMP example. Here, the program runs, and produces reasonnable results. But runnig with 4 threads is 3 to 4 times SLOWER than with 1 thread, depending on REB_GRAVITY_BASIC (3 times) or (REB_GRAVITY_COMPENSATED (4 times). What is going on ?

I copy my example code at end of this message.

Regards, Jean-Marc Petit

PS: in the current example for openmp, the speed up is computed by compraing the time with 1 thread and N threads. But running with 1 thread and OPENMP on is longer than running without OPENMP altogether. So the speeding factor is over estimated.

problem.c code:

/**
 * Solar System
 *
 * This example integrates all planets of the Solar
 * System. The data comes from the NASA HORIZONS system. 
 */
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>
#include <time.h>
#include <string.h>
#include <sys/time.h>
#include "rebound.h"

double ss_pos[10][3] = 
{
    {3.256101656448802E-03  , -1.951205394420489E-04 , -1.478264728548705E-04},  
    {-1.927589645545195E-01 , 2.588788361485397E-01  , 3.900432597062033E-02 }, 
    {-5.976537074581466E-01 , 3.918678996109574E-01  , 3.990356741282203E-02 }, 
    {-7.986189029000561E-01 , -6.086873314992410E-01 , -1.250824315650566E-04}, 
    {7.897942807177173E-01  , 1.266671734964037E+00  , 7.092292179885432E-03 }, 
    {-4.314503046344270E+00 , 3.168094294126697E+00  , 8.331048545353310E-02 }, 
    {-4.882304833383455E+00 , -8.689263067189865E+00 , 3.453930436208210E-01 }, 
    {1.917757033372740E+01  , 5.671738750949031E+00  , -2.273858614425555E-01},  
    {2.767031517959636E+01  , -1.150331645280942E+01 , -4.008018419157927E-01},  
    {7.765250227278298E+00  , -3.190996242617413E+01 , 1.168394015703735E+00 }, 

};
double ss_vel[10][3] = 
{
    {3.039963463108432E-06 ,  6.030576499910942E-06 ,  -7.992931269075703E-08}, 
    {-2.811550184725887E-02,  -1.586532995282261E-02,  1.282829413699522E-03 }, 
    {-1.113090630745269E-02,  -1.703310700277280E-02,  4.089082927733997E-04 },
    {1.012305635253317E-02 ,  -1.376389620972473E-02,  3.482505080431706E-07 }, 
    {-1.135279609707971E-02,  8.579013475676980E-03 ,  4.582774369441005E-04 }, 
    {-4.555986691913995E-03,  -5.727124269621595E-03,  1.257262404884127E-04 }, 
    {4.559352462922572E-03 ,  -2.748632232963112E-03,  -1.337915989241807E-04}, 
    {-1.144087185031310E-03,  3.588282323722787E-03 ,  2.829006644043203E-05 }, 
    {1.183702780101068E-03 ,  2.917115980784960E-03 ,  -8.714411604869349E-05}, 
    {3.112825364672655E-03 ,  1.004673400082409E-04 ,  -9.111652976208292E-04},
};

double ss_mass[10] =
{
    1.988544e30,
    3.302e23,
    48.685e23,
    6.0477246e24,
    6.4185e23,
    1898.13e24,
    5.68319e26,
    86.8103e24,
    102.41e24,
    1.4639248e+22,
};

void heartbeat(struct reb_simulation* r);
double e_init;
double l_init;
double tmax;

void run_sim(){
    struct reb_simulation* const r = reb_create_simulation();
    // Setup constants
    r->dt           = 4;                // in days
    tmax            = 7.3e7;            // 20 Myr
    r->G            = 1.4880826e-34;        // in AU^3 / kg / day^2.
    r->ri_whfast.safe_mode  = 0;        // Turn off safe mode. Need to call reb_integrator_synchronize() before outputs. 
    r->ri_whfast.corrector  = 11;       // 11th order symplectic corrector
    r->integrator       = REB_INTEGRATOR_LEAPFROG;
    r->gravity      = REB_GRAVITY_BASIC;
    //r->integrator     = REB_INTEGRATOR_WHFAST;
    //r->integrator     = REB_INTEGRATOR_WH;
    //r->integrator     = REB_INTEGRATOR_HYBRID;
    //r->gravity        = REB_GRAVITY_COMPENSATED;
    //r->integrator     = REB_INTEGRATOR_IAS15;     // Alternative non-symplectic integrator
    r->heartbeat        = heartbeat;

    // Initial conditions
    for (int i=0;i<10;i++){
        struct reb_particle p = {0};
        p.x  = ss_pos[i][0];        p.y  = ss_pos[i][1];        p.z  = ss_pos[i][2];
        p.vx = ss_vel[i][0];        p.vy = ss_vel[i][1];        p.vz = ss_vel[i][2];
        p.m  = ss_mass[i];
        reb_add(r, p); 
    }
    reb_move_to_com(r);
    reb_tools_energy(r, &e_init, &l_init);
    system("rm -f energy.txt");
    reb_integrate(r, tmax);
    reb_free_simulation(r);
}

int main(int argc, char* argv[]){
    // Get the number of processors
    int np = omp_get_num_procs();
    // Set the number of OpenMP threads to be the number of processors  
    omp_set_num_threads(np);

    struct timeval tim;
    gettimeofday(&tim, NULL);
    double timing1 = tim.tv_sec+(tim.tv_usec/1000000.0);
    run_sim();
    system("mv -f energy.txt energy_4thr.txt");

    gettimeofday(&tim, NULL);
    double timing2 = tim.tv_sec+(tim.tv_usec/1000000.0);
    omp_set_num_threads(1);
    run_sim();
    system("mv -f energy.txt energy_1thr.txt");

    gettimeofday(&tim, NULL);
    double timing3 = tim.tv_sec+(tim.tv_usec/1000000.0);
    printf("Elapse time 1= %- 9f\n", timing2-timing1);
    printf("Elapse time 2= %- 9f\n", timing3-timing2);
    // Output speedup
    //printf("\n\nOpenMP speed-up: %.3fx (perfect scaling would give %dx)\n",(timing3-timing2)/(timing2-timing1),np);
}

void heartbeat(struct reb_simulation* r){
    if (reb_output_check(r, 7300000.)){
        reb_output_timing(r, tmax);
        reb_integrator_synchronize(r);
        FILE* f = fopen("energy.txt","a");
        double e;
        double l;
        reb_tools_energy(r, &e, &l);
        fprintf(f,"%e %e %e\n",r->t, fabs((e-e_init)/e_init), fabs((l-l_init)/l_init));
        fclose(f);
    }
}

Add Travis-CI support

Make a simple travis script to check if all examples compile.

Divergence of Integration on Separate Computers

When I run the code below (a simple simulation of the Kepler-11 system) on two different computers, the output varies on the order of ~10^-9.

import numpy as np
import rebound
from numpy import ndarray 
import matplotlib.pyplot as plt


model='MigaIVa'
integrator = 'ias15'

tmax = 1000. # Final time in years
Noutputs = 1001 # Number of outputs

rebound.integrator=integrator
rebound.G=4.*np.pi**2


K11 = rebound.Particle( m=0.950 )
rebound.add( K11 )
rebound.add( primary=K11, m=7.456e-6, a=0.091113, MEAN=3.122075,
             e=0.04423, omega=0.3604, inc=1.54 )  #K11b
rebound.add( primary=K11, m=1.070e-5, a=0.106505, MEAN=3.658224,
             e=0.01719, omega=0.9726, inc=1.55 )  #K11c
rebound.add( primary=K11, m=1.779e-5, a=0.154243, MEAN=0.711965,
             e=0.00633, omega=2.4566, inc=1.59 )  #K11d
rebound.add( primary=K11, m=3.426e-5, a=0.193940, MEAN=5.559193,
             e=0.00258, omega=4.1323, inc=1.55 )  #K11e
rebound.add( primary=K11, m=2.378e-5, a=0.249511, MEAN=1.598297,
             e=0.01073, omega=6.2107, inc=1.56 )  #K11f
rebound.add( primary=K11, m=7.952e-5, a=0.463991, MEAN=5.868932,
             e=0., inc=1.57, Omega=0. )  #K11g


rebound.move_to_com()

ps = rebound.particles
rebound.dt=0.01*np.sqrt(ps[1].calculate_orbit().a**3)

NumPlanet = rebound.N-1


print rebound.calculate_com()

x = np.zeros((NumPlanet,Noutputs))
times = np.linspace(0.,tmax,Noutputs)


for i,time in enumerate(times):

    rebound.integrate(time)

    for j in range(NumPlanet):
        x[j,i] = ps[j+1].x

Data = open('SimpleIntData.txt', 'a')
[Planet,Time] = np.ma.shape(x)
for t in range(Time):
    output = str(t)+' '
    for p in range(Planet):
        output = output+str(x[p,t])+' '
    Data.write(output+'\n')
Data.write('\n')
Data.close()

Collision resolve can get messed up if multiple mergers per timestep

When a collision is detected, collision.c currently stores the indices of the particles involved in the collision, specifically the indices correspond to their positions in the particles array. In the case of

-multiple collisions in a single timestep, and
-particle mergers

I think I have found an issue. I give a specific example to highlight this. Let's say we have four particles, referred to by their original indices in the particles array: 0,1, 2, and 3. Let's say after a time step there is a collision between 0 and 1 as well as a collision between 2 and 3. The collision between 0 and 1 is resolved first and results in a merger. The information for particle 0 is updated accordingly and particle 1 is removed. If keepSorted=1 in reb_remove(), particles 2 and 3 will move down to fill the gap left by particle 1 (final order: 0, 2, 3. Now it's time to resolve the collision between particles 2 and 3. The collision resolve function will now be looking in the memory at where the particles used to be located in the particles array, not where they are now. Where particle 3 used to be in the array is no longer technically included in the simulation particle, but the particles array still has memory allocated for particle 3's old spot and particle 3's information is still stored there. The collision resolve function will now be trying to resolve a collision between particle 3 (at particle 2's old spot) and particle 3 (at particle 3's old spot, now off the portion of the array that the simulation cares about).

So, it seems that making keepSorted=0 in reb_remove() for the initial removal of particle 1 would fix this. However, it will not. If the collision between particle 2 and 3 leads to a merger, then particle 2's information will be updated based on the merger information, and then reb_remove() will be called on particle 3's old position, which is outside of the array currently recognized by the simulation. reb_remove() will print out the following error:
"\nIndex %d passed to particles_remove was out of range (N=%d). Did not remove particle.\n"
and particle 3 (which is at particle 1's old spot) won't actually be removed from the simulation.

My proposal to fix this:

-Assign each particle a unique ID and rewrite the reb_collision_search() function to give particle ID's to the collision struct. These ID's will follow the particles wherever they go in the array.

The simple fix I propose above would make it so particle ID's would have to be assigned. REBOUND currently does not require particle ID's to be assigned (as far as I know) so the fix proposed above would probably not be ideal for a general fix. Any other ideas?

python documentation page missing?

I've been using rebound and noticed that on the documentation page, the API python documentation suddenly disappeared. The page is now blank. Is that easy to restore? It would be really useful for me to have for reference. Thank you.

Inhomogeneous particles distribution

Particles are gathering at the radial inner border of the box when using the tree.