Integrating arbitrary ODEs (C)

This examples shows how to integrate arbitrary ODEs with REBOUND. In this case we couple a harmonic oscillator to an N-body simulation and drive it using the orbital phase of a planet.

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

const double k = 1.; // Constants for the Harmonic Oscillator
const double m = 1.;

void derivatives(struct reb_ode* const ode, double* const yDot, const double* const y, const double t){
    const double omega = sqrt(k/m);
    struct reb_orbit o = reb_tools_particle_to_orbit(ode->r->G, ode->r->particles[1], ode->r->particles[0]);
    double forcing = sin(o.f);
    yDot[0] = y[1];
    yDot[1] = -omega*omega*y[0] + forcing;
}

int main(int argc, char* argv[]) {
    struct reb_simulation* r = reb_create_simulation();

    reb_add_fmt(r, "m", 1.);                // Central object
    reb_add_fmt(r, "m a e", 1e-3, 1., 0.1); // Jupiter mass planet
    reb_move_to_com(r);

    r->integrator = REB_INTEGRATOR_BS;  // Bulirsch-Stoer integrator
    r->ri_bs.eps_rel = 1e-8;            // Relative tolerance
    r->ri_bs.eps_abs = 1e-8;            // Absolute tolerance
    r->dt = 1e-2;

    struct reb_ode* ho = reb_create_ode(r,2);   // Add an ODE with 2 dimensions
    ho->derivatives = derivatives;              // Right hand side of the ODE
    ho->y[0] = 1;                               // Initial conditions
    ho->y[1] = 0;

    while(r->t<10){
        reb_integrate(r, r->t + 0.3);
        printf("y(%.5f) \t = %.5f \n",r->t, ho->y[0]);
    }


    reb_free_ode(ho);
    reb_free_simulation(r);

}

This example is located in the directory examples/arbitrary_ode