Restricted three body problem. (C)ΒΆ

This example simulates a disk of test particles around a central object, being perturbed by a planet. It uses the heliocentric version of WHFast.

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>
#include "rebound.h"

void heartbeat(struct reb_simulation* r);

int main(int argc, char* argv[]){
    struct reb_simulation* r = reb_create_simulation();
    // Setup constants
    r->integrator   = REB_INTEGRATOR_WHFAST;
    r->ri_whfast.coordinates = REB_WHFAST_COORDINATES_DEMOCRATICHELIOCENTRIC;
    r->boundary     = REB_BOUNDARY_OPEN;
    r->softening    = 1e-6;
    r->dt           = 1.0e-2*2.*M_PI;
    r->N_active     = 2;    // Only the star and the planet have non-zero mass
    r->heartbeat    = heartbeat;

    reb_configure_box(r,8.,1,1,1); // Box with size 8 AU

    // Initial conditions for star
    struct reb_particle star = {0};
    star.m  = 1;
    reb_add(r, star);

    // Initial conditions for planet
    double planet_e = 0.;
    struct reb_particle planet = {0};
    planet.x  = 1.-planet_e;
    planet.vy = sqrt(2./(1.-planet_e)-1.);
    planet.m  = 1e-2;
    reb_add(r, planet);
    reb_move_to_com(r);

    while(r->N<10000){
        double x    = ((double)rand()/(double)RAND_MAX-0.5)*8.;
        double y    = ((double)rand()/(double)RAND_MAX-0.5)*8.;
        double a    = sqrt(x*x+y*y);
        double phi  = atan2(y,x);
        if (a<.1) continue;
        if (a>4.) continue;

        double vkep = sqrt(r->G*star.m/a);
        struct reb_particle testparticle = {0};
        testparticle.x  = x;
        testparticle.y  = y;
        testparticle.z  = 1.0e-2*x*((double)rand()/(double)RAND_MAX-0.5);
        testparticle.vx = -vkep*sin(phi);
        testparticle.vy = vkep*cos(phi);
        reb_add(r, testparticle);
    }

    reb_integrate(r, INFINITY);
}

void heartbeat(struct reb_simulation* r){
    if (reb_output_check(r, 20.*M_PI)){
        reb_output_timing(r, 0);
        reb_output_orbits(r, "orbit.txt");
    }
}

This example is located in the directory examples/restricted_threebody