5. Examples (C)

5.1. Radiation forces on a debris disk

This example shows how to integrate a debris disk around a Sun-like star, with dust particles under the action of radiation forces using WHFast.

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

void heartbeat(struct reb_simulation* r);

double tmax = 3e12;                     // in sec, ~ 10^5 yrs

int main(int argc, char* argv[]){
    struct reb_simulation* sim = reb_simulation_create();
    // Setup constants
    double AU = 1.5e11;                 // in meters
    sim->integrator     = REB_INTEGRATOR_WHFAST;
    sim->G              = 6.674e-11;    // Use SI units
    sim->dt             = 1e8;          // At ~100 AU, orbital periods are ~1000 yrs, so here we use ~1% of that, in sec
    sim->N_active       = 1;            // The dust particles do not interact with one another gravitationally
    sim->heartbeat      = heartbeat;
    sim->usleep         = 1000;             // Slow down integration (for visualization only). Remove in your integrations!

    // sun
    struct reb_particle sun = {0};
    sun.m  = 1.99e30;                   // mass of Sun in kg
    reb_simulation_add(sim, sun);

    struct rebx_extras* rebx = rebx_attach(sim);
    struct rebx_force* rad = rebx_load_force(rebx, "radiation_forces");
    rebx_add_force(rebx, rad);

    rebx_set_param_double(rebx, &rad->ap, "c", 3.e8);

    /* The effect assumes the radiation source is particles[0].  You can set it to another one by adding a radiation_source flag to it:
     * rebx_set_param_int(rebx, &sim->particles[3].ap, "radiation_source", 1);

     * Initialize dust particles
     * We idealize a perfectly coplanar debris disk with particles that have semimajor axes between 100 and 120 AU.
     * We initialize particles with 0.01 eccentricity, random pericenters and azimuths, and uniformly distributed
     * semimajor axes in the above range.  We take all particles to have a beta parameter (ratio of radiation
     * force to gravitational force from the star) of 0.1.
     *
     ***** Only particles that have their beta parameter set will feel radiation forces. *****/

    double amin = 100.*AU;
    double awidth = 20.*AU;
    double e = 0.1;
    double inc = 0.;
    double Omega = 0.;
    double Ndust = 1000;                // Number of dust particles

    int seed = 3;                       // random number generator seed
    srand(seed);
    double a, pomega, f;
    struct reb_particle p;
    double beta = 0.1;
    for(int i=1; i<=Ndust; i++){
                                        // first we set up the orbit from sets of uniformly drawn orbital elements.
        a = amin + awidth*(double)rand() / (double)RAND_MAX;
        pomega = 2*M_PI*(double)rand() / (double)RAND_MAX;
        f = 2*M_PI*(double)rand() / (double)RAND_MAX;

        double m=0;                     // We treat the dust grains as massless.
        p = reb_particle_from_orbit(sim->G, sim->particles[0], m, a, e, inc, Omega, pomega, f);
        reb_simulation_add(sim, p);
                                        // Only particles with beta set will feel radiation forces
        rebx_set_param_double(rebx, &sim->particles[i].ap, "beta", beta);
    }

    reb_simulation_move_to_com(sim);
    reb_simulation_integrate(sim, tmax);

    /* Note that the debris disk will seem to undergo oscillations at first.
     * This is actually the correct behavior.  We set up the
     * particles with orbital elements (that assume only gravity is acting).  When we turn on radiation,
     * all of a sudden particles are moving too fast for the reduced gravity they feel (due to radiation
     * pressure), so they're all at pericenter.  All the particles therefore move outward in phase.  It
     * takes a few cycles before the differential motion randomizes the phases. */

    rebx_free(rebx);
    reb_simulation_free(sim);
}

void heartbeat(struct reb_simulation* sim){
    if(reb_simulation_output_check(sim, 1.e8)){
        //reb_simulation_output_timing(sim, tmax);
    }
}

This example is located in the directory examples/rad_forces_debris_disk

5.2. Type I migration.

This example shows how to add Type I migration. See detailed annotations and explanations for the parameters in the TypeIMigration ipython example and Implemented Effects documentation for REBOUNDx.

#include <stdio.h>
#include <string.h>
#include "rebound.h"
#include "reboundx.h"
#include "core.h"

int main(int argc, char* argv[]){
    struct reb_simulation* sim = reb_simulation_create();
    /* sim units are ('yr', 'AU', 'Msun') */
    sim->G = 4*M_PI*M_PI;

    struct reb_particle star = {0};
    star.m     = 1.;
    reb_simulation_add(sim, star);

    double m = 0.00001;
    double a1 = 0.5;
    double a2 = 0.85;
    double e = 0;
    double inc = 0.;
    double Omega = 0.;
    double omega = 0.;
    double f = 0.;

    struct reb_particle p1 = reb_particle_from_orbit(sim->G, star, m, a1, e, inc, Omega, omega, f);
    struct reb_particle p2 = reb_particle_from_orbit(sim->G, star, m, a2, e, inc, Omega, omega, f);
    reb_simulation_add(sim, p1);
    reb_simulation_add(sim, p2);
    reb_simulation_move_to_com(sim);

    sim->dt = 0.002;  //The period at inner disk edge divided by 20, for a disk edge location at 0.1 AU
    sim->integrator = REB_INTEGRATOR_WHFAST;

    struct rebx_extras* rebx = rebx_attach(sim);

    struct rebx_force* type_I_mig = rebx_load_force(rebx, "type_I_migration");
    rebx_set_param_double(rebx, &type_I_mig->ap, "ide_position", 0.3);
    rebx_set_param_double(rebx, &type_I_mig->ap, "ide_width", 0.02); //Calculated using the scale height value given below
    rebx_set_param_double(rebx, &type_I_mig->ap, "tIm_flaring_index", 0.25);
    rebx_set_param_double(rebx, &type_I_mig->ap, "tIm_surface_density_exponent", 1);
    rebx_set_param_double(rebx, &type_I_mig->ap, "tIm_surface_density_1", 0.00011255); //In units of solarmass/(AU^2) which is roughly 1000g/cm^2
    rebx_set_param_double(rebx, &type_I_mig->ap, "tIm_scale_height_1", 0.03);
    rebx_add_force(rebx, type_I_mig);

    double tmax = 1.e4;  // Note that one can calculate the timescale of the semi-major axis of the outer planet then set the integration time to twice this value

    reb_simulation_integrate(sim, tmax);
    rebx_free(rebx);
    reb_simulation_free(sim);
}

This example is located in the directory examples/type_I_migration

5.3. Adding custom post-timestep modifications and forces.

This allows the user to use the built-in functions of REBOUNDx but also include their own specialised functions.

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

/* There are two qualitatively different options for adding in effects. One can add a force, which evaluates accelerations to add
 * on top of gravitational accelerations every timestep and will be integrated numerically.
 * Alternatively one can apply operators before and/or after each REBOUND timestep that update particle properties
 * (masses, positions/velocities, or other parameters).
 *
 * Here's an example of a force, for whichi we have to update the particle accelerations in the passed particles array (not sim->particles!)
 * Function must have the same prototype as below.
 * We can have the user assign parameters to the force structure in main(), and our force function read and/or update those parameters each timestep when called.

void stark_force(struct reb_simulation* const sim, struct rebx_force* const starkforce, struct reb_particle* const particles, const int N){
    double* starkconst = rebx_get_param(sim->extras, starkforce->ap, "starkconst");    // get parameters we want user to set

    if(starkconst != NULL){
        particles[1].ax += (*starkconst);           // make sure you += not =, which would overwrite other accelerations
    }
}

/* Here's a very simple example of an operator to change the planet's orbit
 *
 * For a post_timestep_modification we update the particle states (positions, velocities, masses etc.)
 * Function must have the same prototype as below, passing simulation and operator pointers (which can be used to hold parameters),
 * and the timestep over which the operator should act.

void simple_drag(struct reb_simulation* const sim, struct rebx_operator* const dragoperator, const double dt){
    double* dragconst = rebx_get_param(sim->extras, dragoperator->ap, "dragconst");    // get parameters we want user to set

    if(dragconst != NULL){
        sim->particles[1].vx *= 1. - (*dragconst)*dt;
        sim->particles[1].vy *= 1. - (*dragconst)*dt;
        sim->particles[1].vz *= 1. - (*dragconst)*dt;
    }
}


int main(int argc, char* argv[]){
    struct reb_simulation* sim = reb_simulation_create();
    struct reb_particle p = {0};
    p.m     = 1.;
    reb_simulation_add(sim, p);

    double m = 0.;
    double a1 = 1.;
    double e = 0.4;
    double inc = 0.;
    double Omega = 0.;
    double omega = 0.;
    double f = 0.;

    struct reb_particle p1 = reb_particle_from_orbit(sim->G, p, m, a1, e, inc, Omega, omega, f);
    reb_simulation_add(sim,p1);
    reb_simulation_move_to_com(sim);

    struct rebx_extras* rebx = rebx_attach(sim);  // first initialize rebx

    /* We now add our custom functions. We do this in the same way as we add REBOUNDx forces and operators.
     * We'll still get back a force and operator struct, respectively, with a warning that the name wasn't found
     * in REBOUNDx and that we're reponsible for setting their type flags and function pointers.
     */

    struct rebx_force* stark = rebx_create_force(rebx, "stark_force");

    /* We first set a flag for whether our force only depends on particle positions (REBX_FORCE_POS)
     * or whether it depends on velocities (or velocities and positions, REBX_FORCE_VEL)
     */
    stark->force_type = REBX_FORCE_VEL;
    stark->update_accelerations = stark_force;  // set the function pointer to what we wrote above
    rebx_add_force(rebx, stark);                // Now it's initialized, add to REBOUNDx

    struct rebx_operator* drag = rebx_create_operator(rebx, "simple_drag"); // Now we create our custom operator

    /* We first set the operator_type enum for whether our operator modifies dynamical variables (positions, velocities
     * or masses, REBX_OPERATOR_UPDATER), or whether it's just passively recording the state of the simulation, or
     * updating parameters that don't feed back on the dynamics (REBX_OPERATOR_RECORDER)
     */

    drag->operator_type = REBX_OPERATOR_UPDATER;
    drag->step_function = simple_drag;  // set function pointer to what we wrote above
    rebx_add_operator(rebx, drag);      // Now it's initialized, add to REBOUNDx

    /* Now we set the parameters that our custom functions above need
     * Before setting them, we need to register them with their corresponding types
     */

    rebx_register_param(rebx, "starkconst", REBX_TYPE_DOUBLE);
    rebx_register_param(rebx, "dragconst", REBX_TYPE_DOUBLE);

    rebx_set_param_double(rebx, &stark->ap, "starkconst", 1.e-5);
    rebx_set_param_double(rebx, &drag->ap, "dragconst", 1.e-5);

    /* To simplify things for other users, we can always add the new effects and register parameters in the REBOUNDx
     * repo itself. Feel free to send a pull request or contact me (tamayo.daniel@gmail.com) about adding it
     */

    double tmax = 5.e4;
    reb_simulation_integrate(sim, tmax);
    rebx_free(rebx);                            // Free all the memory allocated by rebx
}

This example is located in the directory examples/custom_effects

5.4. Gas dynamical friction

This example shows how to add the gas_df force.

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

int main(int argc, char* argv[]){
    struct reb_simulation* sim = reb_simulation_create();
    sim->integrator = REB_INTEGRATOR_BS;

    struct reb_particle bh = {0};
    bh.m     = 4e6;
    reb_simulation_add(sim, bh);

    double m = 1;
    double a = 206000;
    double e = 0.01;
    double inc = 0.17;
    double Omega = 0.;
    double omega = 0.;
    double f = 0.;

    struct reb_particle star = reb_particle_from_orbit(sim->G, bh, m, a, e, inc, Omega, omega, f);
    reb_simulation_add(sim, star);
    reb_simulation_move_to_com(sim);

    struct rebx_extras* rebx = rebx_attach(sim);
    struct rebx_force* gdf = rebx_load_force(rebx, "gas_dynamical_friction");
    rebx_add_force(rebx, gdf);

    rebx_set_param_double(rebx, &gdf->ap, "gas_df_rhog", 0.20);
    rebx_set_param_double(rebx, &gdf->ap, "gas_df_alpha_rhog", -1.5);
    rebx_set_param_double(rebx, &gdf->ap, "gas_df_cs", 20);
    rebx_set_param_double(rebx, &gdf->ap, "gas_df_alpha_cs", -0.5);
    rebx_set_param_double(rebx, &gdf->ap, "gas_df_xmin", 0.045);
    rebx_set_param_double(rebx, &gdf->ap, "gas_df_hr", 0.01);
    rebx_set_param_double(rebx, &gdf->ap, "gas_df_Qd", 5.0);

    double delta_t = 6.28e5;
    for (int i = 0; i < 100; i++){
        reb_simulation_integrate(sim, sim->t + delta_t);
        struct reb_orbit o = reb_orbit_from_particle(sim->G, sim->particles[1], sim->particles[0]);
        printf("%f %f %f %e\n", sim->t, o.a, o.e, o.inc);
    }

    rebx_free(rebx);    // this explicitly frees all the memory allocated by REBOUNDx
    reb_simulation_free(sim);
}

This example is located in the directory examples/gas_dynamical_friction

5.5. Inner disk edge.

This example shows how to add an inner disk edge when using modify_orbits_forces or modify_orbits_direct. See detailed annotations and explanations for the parameters in the InnerDiskEdge ipython example and Implemented Effects documentation for REBOUNDx.

#include <stdio.h>
#include <string.h>
#include "rebound.h"
#include "reboundx.h"
#include "core.h"

int main(int argc, char* argv[]){
    struct reb_simulation* sim = reb_simulation_create();
    /*sim units are ('yr', 'AU', 'Msun')*/
    sim->G = 4*M_PI*M_PI;

    struct reb_particle star = {0};
    star.m     = 1.;
    reb_simulation_add(sim, star);

    double m = 0.00001;
    double a = 1;
    double e = 0;
    double inc = 0.;
    double Omega = 0.;
    double omega = 0.;
    double f = 0.;

    struct reb_particle p1 = reb_particle_from_orbit(sim->G, star, m, a, e, inc, Omega, omega, f);
    reb_simulation_add(sim, p1);
    reb_simulation_move_to_com(sim);

    sim->dt = 0.002;  //The period at the inner disk edge divided by 20, for a disk edge location at 0.1 AU
    sim->integrator = REB_INTEGRATOR_WHFAST;

    struct rebx_extras* rebx = rebx_attach(sim);

    struct rebx_force* mof = rebx_load_force(rebx, "modify_orbits_forces");
    rebx_set_param_double(rebx, &mof->ap, "ide_position", 0.1);
    rebx_set_param_double(rebx, &mof->ap, "ide_width", 0.02); // Planet will stop within 0.02 AU of the inner disk edge
    rebx_add_force(rebx, mof);

    double tmax = 1.e4;
    rebx_set_param_double(rebx, &sim->particles[1].ap, "tau_a", -tmax);
    rebx_set_param_double(rebx, &sim->particles[1].ap, "tau_e", -tmax/100.);

    reb_simulation_integrate(sim, tmax);
    rebx_free(rebx);
    reb_simulation_free(sim);
}

This example is located in the directory examples/inner_disk_edge

5.6. Exponential Migration

This example shows how to add exponential migration to a REBOUND simulation. Both the Reboundx force (exponential_migration.c) and this example are based on ``modify_orbits_forces’’

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

void heartbeat(struct reb_simulation* sim);

int main(int argc, char* argv[]){
    struct reb_simulation* sim = reb_simulation_create();
    // Setup constants
    sim->dt             = 0.012;        // initial timestep.
    sim->heartbeat = heartbeat;

    struct reb_particle p = {0};
    p.m     = 1.;
    reb_simulation_add(sim, p);

    double m = 0.0001;
    double a1 = 24.0;
    double e = 0.01;
    double inc = 0.;
    double Omega = 0.;
    double omega = 0.;
    double f = 0.;

    struct reb_particle p1 = reb_particle_from_orbit(sim->G, p, m, a1, e, inc, Omega, omega, f);
    reb_simulation_add(sim,p1);
    reb_simulation_move_to_com(sim);

    struct rebx_extras* rebx = rebx_attach(sim);

    struct rebx_force* em = rebx_load_force(rebx, "exponential_migration");  // add force that adds velocity kicks to give exponential migration
    rebx_add_force(rebx, em);

    // Set the timescales for each particle.
    double tmax = 4.e4;

    rebx_set_param_double(rebx, &sim->particles[1].ap, "em_tau_a", 2.e3);         // add migration e-folding timescale
    rebx_set_param_double(rebx, &sim->particles[1].ap, "em_aini", 24.0);         // add initial semimajor axis
    rebx_set_param_double(rebx, &sim->particles[1].ap, "em_afin", 30.0);         // add final semimajor axis

     rebx_set_param_int(rebx, &em->ap, "coordinates", REBX_COORDINATES_PARTICLE);
     rebx_set_param_int(rebx, &sim->particles[0].ap, "primary", 1);

    reb_simulation_integrate(sim, tmax);
    rebx_free(rebx);    // Free all the memory allocated by rebx
    reb_simulation_free(sim);
}

void heartbeat(struct reb_simulation* sim){
    // output a e body
    if(reb_simulation_output_check(sim, 1.e2)){
        const struct reb_particle sun = sim->particles[0];
        const struct reb_orbit orbit = reb_orbit_from_particle(sim->G, sim->particles[1], sun); // calculate orbit of particles[1]
        printf("%f\t%f\t%f\t%f\n",sim->t,orbit.a, orbit.e);
    }
}

This example is located in the directory examples/exponential_migration

5.7. Post-Newtonian correction from general relativity

This example shows how to add post-newtonian corrections to REBOUND simulations with reboundx. If you have GLUT installed for the visualization, press ‘w’ and/or ‘c’ for a clearer view of the whole orbit.

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

int main(int argc, char* argv[]){
    struct reb_simulation* sim = reb_simulation_create();
    sim->dt = 1.e-8;

    struct reb_particle star = {0};
    star.m     = 1.;
    star.hash  = reb_hash("star");
    reb_simulation_add(sim, star);

    double m = 1.e-5;
    double a = 1.e-4; // put planet close to enhance precession so it's visible in visualization (this would put planet inside the Sun!)
    double e = 0.2;
    double inc = 0.;
    double Omega = 0.;
    double omega = 0.;
    double f = 0.;

    struct reb_particle planet = reb_particle_from_orbit(sim->G, star, m, a, e, inc, Omega, omega, f);
    planet.hash = reb_hash("planet");
    reb_simulation_add(sim, planet);
    reb_simulation_move_to_com(sim);

    struct rebx_extras* rebx = rebx_attach(sim);
    // Could also add "gr" or "gr_full" here.  See documentation for details.
    struct rebx_force* gr = rebx_load_force(rebx, "gr");
    rebx_add_force(rebx, gr);
    // Have to set speed of light in right units (set by G & initial conditions).  Here we use default units of AU/(yr/2pi)
    rebx_set_param_double(rebx, &gr->ap, "c", 10065.32);

    double tmax = 5.e-1;
    reb_simulation_integrate(sim, tmax);
    rebx_free(rebx);    // this explicitly frees all the memory allocated by REBOUNDx
    reb_simulation_free(sim);
}

This example is located in the directory examples/gr

5.8. Adding Post-Newtonian correction from general relativity as an operator.

This example shows how to add post-newtonian corrections to REBOUND simulations as an operator (see REBOUNDx paper and corresponding IntegrateForce.ipynb jupyter notebook for more details). If you have GLUT installed for the visualization, press ‘w’ and/or ‘c’ for a clearer view of the whole orbit.

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

double E0;
void heartbeat(struct reb_simulation* sim);

int main(int argc, char* argv[]){
    // We first set up a system similar to the EPIC system discussed in the REBOUNDx paper
    struct reb_simulation* sim = reb_simulation_create();
    sim->G = 4*M_PI*M_PI;

    struct reb_particle star = {0};
    star.m     = 0.93;
    reb_simulation_add(sim, star);

    double m = 1.35e-6;
    double a = 0.013; // put planet close to enhance precession so it's visible in visualization (this would put planet inside the Sun!)
    double e = 0.01;
    double inc = 0.;
    double Omega = 0.;
    double omega = 0.;
    double f = 0.;

    struct reb_particle planet1 = reb_particle_from_orbit(sim->G, star, m, a, e, inc, Omega, omega, f);

    m = 1.2e-5;
    a = 0.107; // put planet close to enhance precession so it's visible in visualization (this would put planet inside the Sun!)
    e = 0.01;

    struct reb_particle planet2= reb_particle_from_orbit(sim->G, star, m, a, e, inc, Omega, omega, f);
    reb_simulation_add(sim, planet1);
    reb_simulation_add(sim, planet2);
    reb_simulation_move_to_com(sim);

    sim->dt = 1.e-4;
    sim->heartbeat = heartbeat;
    sim->integrator = REB_INTEGRATOR_WHFAST;

    // Now we add GR. This is a velocity dependent force. With WHFast this would cause errors on long timescales, so we integrate the force in a separate step
    // See the REBOUNDx paper and the corresponding example in ipython_examples for more details
    // By adding a separate force step for GR and integrating across it we keep an oscillatory energy error

    struct rebx_extras* rebx = rebx_attach(sim);
    struct rebx_force* gr = rebx_load_force(rebx, "gr");
    rebx_set_param_double(rebx, &gr->ap, "c", 63197.8); // AU/yr

    struct rebx_operator* integrateforce = rebx_load_operator(rebx, "integrate_force");
    rebx_set_param_pointer(rebx, &integrateforce->ap, "force", gr);
    rebx_set_param_int(rebx, &integrateforce->ap, "integrator", REBX_INTEGRATOR_RK2);
    rebx_add_operator(rebx, integrateforce);

    double tmax = 10.;
    E0 = rebx_gr_hamiltonian(rebx, gr);
    reb_simulation_integrate(sim, tmax);
    rebx_free(rebx);    // this explicitly frees all the memory allocated by REBOUNDx
    reb_simulation_free(sim);
}

void heartbeat(struct reb_simulation* sim){
    if(reb_simulation_output_check(sim, 0.1)){
        struct rebx_force* gr = rebx_get_force(sim->extras, "gr");
        double E = rebx_gr_hamiltonian(sim->extras, gr);
        printf("t=%f\tEnergy Error=%e\n", sim->t, fabs((E-E0)/E0));
    }
}

This example is located in the directory examples/integrate_force

5.9. Tracking a particle’s minimum distance from the central star.

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

int main(int argc, char* argv[]){
    struct reb_simulation* sim = reb_simulation_create();
    struct rebx_extras* rebx = rebx_attach(sim);

    struct reb_particle star = {0};
    star.m     = 1.;
    star.hash  = reb_hash("star");
    reb_simulation_add(sim, star);

    double m = 0.;
    double a = 1.;
    double e = 0.9;
    double inc = 0.;
    double Omega = 0.;
    double omega = 0.;
    double f = M_PI;

    struct reb_particle planet = reb_particle_from_orbit(sim->G, star, m, a, e, inc, Omega, omega, f);
    reb_simulation_add(sim, planet);
    reb_simulation_move_to_com(sim);

    struct rebx_operator* track_min_distance = rebx_load_operator(rebx, "track_min_distance");
    rebx_add_operator(rebx, track_min_distance);

    // Wee add a min_distance parameter to the particle whose distance we want to track, and set it
    // to a particular value. In any timestep that the distance drops below this value, the value is updated.
    rebx_set_param_double(rebx, &sim->particles[1].ap, "min_distance", 5.);

    // By default distance is measured from sim->particles[0].  We can specify a different particle by a hash (unnecessary here):

    rebx_set_param_uint32(rebx, &sim->particles[1].ap, "min_distance_from", sim->particles[0].hash);


    struct reb_orbit orbit = {0};
    rebx_set_param_pointer(rebx, &sim->particles[1].ap, "min_distance_orbit", &orbit);

    double tmax = 10.;
    reb_simulation_integrate(sim, tmax);

    // At any point in the integration, we can check the `min_distance` parameter and output it as needed.
    double* min_distance = rebx_get_param(rebx, sim->particles[1].ap, "min_distance");
    struct reb_orbit* orbitptr = rebx_get_param(rebx, sim->particles[1].ap, "min_distance_orbit");

    printf("Particle's minimum distance from the star over the integration = %e\n", *min_distance);
    printf("Semimajor axis and eccentricity at closest approach: a=%.2f, e=%.2f\n", orbitptr->a, orbitptr->e);

    rebx_free(rebx);    // this explicitly frees all the memory allocated by REBOUNDx
    reb_simulation_free(sim);
}

This example is located in the directory examples/track_min_distance

5.10. Adding parameters to objects in REBOUNDx

This example walks through adding parameters to particles, forces and operators in REBOUNDx

#include <stdio.h>
#include <math.h>
#include "rebound.h"
#include "reboundx.h"

int main(int argc, char* argv[]){
    // We start by making a simulation and adding GR with REBOUNDx
    struct reb_simulation* sim = reb_simulation_create();

    struct reb_particle star = {0};
    star.m     = 1.;
    reb_simulation_add(sim, star);

    struct reb_particle planet = {0};
    planet.x = 1.;
    planet.vy = 1.;

    reb_simulation_add(sim, planet);

    struct rebx_extras* rebx = rebx_attach(sim);
    struct rebx_force* gr = rebx_load_force(rebx, "gr");
    rebx_add_force(rebx, gr);

    /* The documentation page https://reboundx.readthedocs.io/en/latest/effects.html lists the various required and optional parameters that need to be set for each effect in REBOUNDx.
     *
     * Parameters are stored in particles, forces and operators through the linked list 'ap', which is common to all three objects.
     *
     * For simple types (int and double) we add them through their corresponding setters.
     * In both cases, we pass the rebx struct, a pointer to the head of the linked list that we want to modify, and the value we want to set.
     */

    double c = 10064.915; // speed of light in default units of AU, Msun and yr/2pi
    rebx_set_param_double(rebx, &gr->ap, "c", c);

    // After setting the parameters we want to set, we would integrate as usual.

    double tmax = 10.;
    reb_simulation_integrate(sim, tmax);

    /* At any point, we can access the parameters we set (e.g., some effects could update these values as the simulation progresses). We get all parameter types back with rebx_get_param, which returns a void pointer that we are responsible for casting to the correct type. Since we are not modifying the linked list, we don't pass a reference to ap like above*/

    double* new_c = rebx_get_param(rebx, gr->ap, "c");

    /* It is important to check returned parameters!
     * If the parameter is not found, it will return a NULL pointer.
     * If this happens and you dereference it, you will get a segmentation fault.
     * If you get a seg fault, the first thing you should check are returned pointers.
     */

    if (new_c != NULL){
        printf("c=%f\n", *new_c);
    }

    /* The above functionality is probably enough for most users.
     * If you find you need to do more complicated things, read below!
     *
     * More complicated types can be added with set_param_pointer.
     * Instead of passing a value, we now pass a pointer to the variable.
     * As a simple example adding the gr pointer to particles[1] (with no meaning whatsoever):
     */

    rebx_set_param_pointer(rebx, &sim->particles[1].ap, "force", gr);

    /* In order to be able to write binary files, and for interoperability with Python, REBOUNDx keeps a registered list of parameter names and their corresponding types.
     * Above we have used parameter names that have been registered by various effects.
     * If you try to use a name that's not on the list, REBOUNDx will print an error and continue execution without adding the parameter.
     */

    rebx_set_param_int(rebx, &gr->ap, "q", 7);

    int* q = rebx_get_param(rebx, gr->ap, "q");

    if (q != NULL){
        printf("q=%d\n", *q);
    }
    else{
        printf("q is NULL because we didn't register the name first.\n");
    }

    /* You could register the name permanently in rebx_register_params in core.c, or you can do it manually in problem.c, passing a name and rebx_param_type enum (defined in reboundx.h):
     */

    rebx_register_param(rebx, "q", REBX_TYPE_INT);

    rebx_set_param_int(rebx, &gr->ap, "q", 7);
    q = rebx_get_param(rebx, gr->ap, "q");

    if (q != NULL){
        printf("q=%d\n", *q);
    }
    else{
        printf("q is NULL because we didn't register the name first.\n");
    }

    /* Finally, you might want to add your own custom structs (e.g. from another library).
     * This can be done straightforwardly by adding it as a pointer.
     */

    struct SPH_sim{
        double dt;
        int Nparticles;
    };

    struct SPH_sim my_sph_sim = {.dt = 0.1, .Nparticles=10000};

    rebx_register_param(rebx, "sph", REBX_TYPE_POINTER);
    rebx_set_param_pointer(rebx, &gr->ap, "sph", &my_sph_sim);

    // If you need to get it back:

    struct SPH_sim* sph = rebx_get_param(rebx, gr->ap, "sph");

    if (sph != NULL){
        printf("dt=%f\tNparticles=%d\n", sph->dt, sph->Nparticles);
    }

    /* One caveat is that since REBOUNDx does not know the structure definition, custom parameters will not be written or read from REBOUNDx binary files*/

    rebx_free(rebx);    // this explicitly frees all the memory allocated by REBOUNDx
    reb_simulation_free(sim);
}

This example is located in the directory examples/parameters

5.11. Yarkovsky effect on a small body

This example simulates a single asteroid at .5 AU orbiting around the Sun for 100,000 years to demonstartate how the Yarkovsky effect can change the semi-major axis of an orbiting body. Changing the value of the ‘yark_flag’ parameter between -1, 0, and 1 switches which version of the effect is being used. For more information on the different versions of the effect and what they’re good for, please visit the ipython example and the documentation for this effect.

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

int main(int argc, char* argv[]) {

    struct reb_simulation* sim = reb_simulation_create(); //creates simulation

    sim->G = 4*M_PI*M_PI;  // use units of AU, yr and solar masses
    sim->dt = .05;         //timestep for simulation in yrs
    sim->integrator = REB_INTEGRATOR_WHFAST; //integrator for sim

    //following adds star with mass of Sun to sim
    struct reb_particle star = {0};
    star.m = 1.;
    reb_simulation_add(sim, star);

    //following variables are the orbital elements of only asteroid in sim
    double m = 0;
    double a = .5;
    double e = 0;
    double inc = 0;
    double Omega = 0;
    double omega = 0;
    double f = 0;

    //creates asteroid
    struct reb_particle asteroid_1 = reb_particle_from_orbit(sim->G, star, m, a, e, inc, Omega, omega, f);

    reb_simulation_add(sim,asteroid_1);

    struct reb_particle* const particles = sim->particles; //pointer for the particles in the sim

    struct rebx_extras* rebx = rebx_attach(sim);
    struct rebx_force* yark = rebx_load_force(rebx, "yarkovsky_effect");

    double au_conv = 1.495978707e11;
    double msun_conv = 1.9885e30;
    double yr_conv = 31557600.0;

    double density = (3000.0*au_conv*au_conv*au_conv)/msun_conv;
    double c = (2.998e8*yr_conv)/au_conv;
    double lstar = (3.828e26*yr_conv*yr_conv*yr_conv)/(msun_conv*au_conv*au_conv);
    double albedo = .017;
    double stef_boltz = ((5.670e-8)*yr_conv*yr_conv*yr_conv)/(msun_conv);
    double emissivity = .9;
    double k = .25;
    double Gamma = (310*sqrt(yr_conv)*yr_conv*yr_conv)/msun_conv;
    double rotation_period = 15470.9/yr_conv;
    double sx = 0.0;
    double sy = 0.0;
    double sz = 1.0;

    rebx_set_param_double(rebx, &sim->particles[1].ap, "ye_body_density", density);
    rebx_set_param_int(rebx, &sim->particles[1].ap, "ye_flag", 0);
    rebx_set_param_double(rebx, &yark->ap, "ye_lstar", lstar);
    rebx_set_param_double(rebx, &yark->ap, "ye_c", c);
    particles[1].r = 1000/au_conv;

    rebx_set_param_double(rebx, &sim->particles[1].ap, "ye_albedo", albedo);
    rebx_set_param_double(rebx, &yark->ap, "ye_stef_boltz", stef_boltz);
    rebx_set_param_double(rebx, &sim->particles[1].ap, "ye_emissivity", emissivity);
    rebx_set_param_double(rebx, &sim->particles[1].ap, "ye_k", k);
    rebx_set_param_double(rebx, &sim->particles[1].ap, "ye_thermal_inertia", Gamma);
    rebx_set_param_double(rebx, &sim->particles[1].ap, "ye_rotation_period", rotation_period);
    rebx_set_param_double(rebx, &sim->particles[1].ap, "ye_spin_axis_x", sx);
    rebx_set_param_double(rebx, &sim->particles[1].ap, "ye_spin_axis_y", sy);
    rebx_set_param_double(rebx, &sim->particles[1].ap, "ye_spin_axis_z", sz);

    rebx_add_force(rebx, yark);

    double tmax = 100000;

    reb_simulation_integrate(sim, tmax); //integrates system for tmax years

    struct reb_orbit o= reb_orbit_from_particle(sim->G, sim->particles[1], sim->particles[0]); //gives orbital parameters for asteroid after sim

    double final_a = o.a; //final semi-major axis of asteroid after sim
    double final_e = o.e;

    printf("CHANGE IN SEMI-MAJOR AXIS: %1.30f\n", (final_a-a)); //prints difference between the intitial and final semi-major axes of asteroid

    printf("CHANGE IN ECCENTRICITY: %1.30f\n", (final_e-e));

}

This example is located in the directory examples/yarkovsky_effect

5.12. Kozai cycles

This example uses the IAS15 integrator to simulate a Lidov Kozai cycle of a planet perturbed by a distant star. The integrator automatically adjusts the timestep so that even very high eccentricity encounters are resolved with high accuracy.

This is the same Kozai example implemented in Lu et. al (2023) Also, see the ipython examples prefixed TidesSpin for in-depth exploration of the parameters that can be set in this simulation.

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

void heartbeat(struct reb_simulation* r);
double tmax = 1e5; // kept short to run quickly.
                   // set to 3e5 for a full cycle
                   // or 7e6 * 2 * M_PI to reproduce the paper plot

int main(int argc, char* argv[]){
    struct reb_simulation* sim = reb_simulation_create();
    // Initial conditions
    // Setup constants
    sim->dt                 = M_PI*1e-1;     // initial timestep
    sim->integrator         = REB_INTEGRATOR_IAS15; // IAS15 is used for its adaptive timestep:
                                                   // in a Kozai cycle the planet experiences close encounters during the high-eccentricity epochs.
                                                   // A fixed-time integrator (for example, WHFast) would need to apply the worst-case timestep to the whole simulation
    sim->heartbeat          = heartbeat;

    // Initial conditions
    struct reb_particle star = {0};
    star.m  = 1;
    star.r = 0.00465;
    reb_simulation_add(sim, star);

    // struct reb_particle planet = {0};
    double planet_m  = 0.054 * 9.55e-4; // A Jupiter-like planet
    double planet_r = 0.3 * 4.676e-4;
    double planet_a = 2.;
    double planet_e = 0.001;
    reb_simulation_add_fmt(sim, "m r a e", planet_m, planet_r, planet_a, planet_e);

    // The perturber - treated as a point particle
    double perturber_m  = 1;
    double perturber_a = 50.;
    double perturber_e = 0.7 * M_PI / 180.;
    double perturber_inc = 80. * M_PI / 180.;
    reb_simulation_add_fmt(sim, "m a e inc", perturber_m, perturber_a, perturber_e, perturber_inc);

    // Add REBOUNDx effects
    // First tides_spin
    struct rebx_extras* rebx = rebx_attach(sim);

    struct rebx_force* effect = rebx_load_force(rebx, "tides_spin");
    rebx_add_force(rebx, effect);


    // Sun
    const double solar_k2 = 0.07;
    rebx_set_param_double(rebx, &sim->particles[0].ap, "k2", solar_k2);
    rebx_set_param_double(rebx, &sim->particles[0].ap, "I", 0.07 * star.m * star.r * star.r);

    const double solar_spin_period = 4.6 * 2. * M_PI / 365.;
    const double solar_spin = (2 * M_PI) / solar_spin_period;
    rebx_set_param_vec3d(rebx, &sim->particles[0].ap, "Omega", (struct reb_vec3d){.z=solar_spin}); // Omega_x = Omega_y = 0 by default

    const double solar_Q = 1e6;
    struct reb_orbit orb = reb_orbit_from_particle(sim->G, sim->particles[1], sim->particles[0]);
    // In the case of a spin that is synchronous with a circular orbit, tau is related to the tidal quality factor Q through the orbital mean motion n (see Lu et al. 2023 for discussion). Clearly that's not the case here, but gives us a reasonable starting point to start turning this knob
    double solar_tau = 1 / (2 * solar_Q * orb.n);
    rebx_set_param_double(rebx, &sim->particles[0].ap, "tau", solar_tau);

    // Planet
    const double planet_k2 = 0.4;
    rebx_set_param_double(rebx, &sim->particles[1].ap, "k2", planet_k2);
    rebx_set_param_double(rebx, &sim->particles[1].ap, "I", 0.25 * planet_m * planet_r * planet_r);

    const double spin_period_p = 1. * 2. * M_PI / 365.; // days to reb years
    const double spin_p = (2. * M_PI) / spin_period_p;
    const double theta_p = 0. * M_PI / 180.;
    const double phi_p = 0. * M_PI / 180;
    struct reb_vec3d Omega_sv = reb_tools_spherical_to_xyz(spin_p, theta_p, phi_p);
    rebx_set_param_vec3d(rebx, &sim->particles[1].ap, "Omega", Omega_sv);

    const double planet_Q = 3e5;
    rebx_set_param_double(rebx, &sim->particles[1].ap, "tau", 1./(2.*planet_Q*orb.n));


    // add GR precession:
    struct rebx_force* gr = rebx_load_force(rebx, "gr_potential");
    rebx_add_force(rebx, gr);
    rebx_set_param_double(rebx, &gr->ap, "c", 10065.32); // in default units

    reb_simulation_move_to_com(sim);

    // Let's create a reb_rotation object that rotates to new axes with newz pointing along the total ang. momentum, and x along the line of
    // nodes with the invariable plane (along z cross newz)
    struct reb_vec3d newz = reb_vec3d_add(reb_simulation_angular_momentum(sim), rebx_tools_spin_angular_momentum(rebx));
    struct reb_vec3d newx = reb_vec3d_cross((struct reb_vec3d){.z =1}, newz);
    struct reb_rotation rot = reb_rotation_init_to_new_axes(newz, newx);
    rebx_simulation_irotate(rebx, rot); // This rotates our simulation into the invariable plane aligned with the total ang. momentum (including spin)

    rebx_spin_initialize_ode(rebx, effect);

    system("rm -v output.txt");        // delete previous output file
    reb_simulation_integrate(sim, tmax);

    rebx_free(rebx);
    reb_simulation_free(sim);
}

void heartbeat(struct reb_simulation* sim){
    // Output spin and orbital information to file
    if(reb_simulation_output_check(sim, 10)){        // outputs every 10 REBOUND years
      struct rebx_extras* const rebx = sim->extras;
      FILE* of = fopen("output.txt", "a");
      if (of==NULL){
          reb_simulation_error(sim, "Can not open file.");
          return;
      }

      struct reb_particle* sun = &sim->particles[0];
      struct reb_particle* p1 = &sim->particles[1];
      struct reb_particle* pert = &sim->particles[2];

      // orbits
      struct reb_orbit o1 = reb_orbit_from_particle(sim->G, *p1, *sun);
      double a1 = o1.a;
      double e1 = o1.e;
      double i1 = o1.inc;
      double Om1 = o1.Omega;
      double pom1 = o1.pomega;
      struct reb_vec3d n1 = o1.hvec;

      struct reb_particle com = reb_particle_com_of_pair(sim->particles[0],sim->particles[1]);
      struct reb_orbit o2 = reb_orbit_from_particle(sim->G, *pert, com);
      double a2 = o2.a;
      double e2 = o2.e;
      double i2 = o2.inc;
      double Om2 = o2.Omega;
      double pom2 = o2.pomega;

      struct reb_vec3d* Omega_sun = rebx_get_param(rebx, sun->ap, "Omega");

      // Interpret planet spin in the rotating planet frame
      struct reb_vec3d* Omega_p_inv = rebx_get_param(rebx, p1->ap, "Omega");

      // Transform spin vector into planet frame, w/ z-axis aligned with orbit normal and x-axis aligned with line of nodes
      struct reb_vec3d line_of_nodes = reb_vec3d_cross((struct reb_vec3d){.z =1}, n1);
      struct reb_rotation rot = reb_rotation_init_to_new_axes(n1, line_of_nodes); // Arguments to this function are the new z and x axes
      struct reb_vec3d srot = reb_vec3d_rotate(*Omega_p_inv, rot); // spin vector in the planet's frame

      // Interpret the spin axis in the more natural spherical coordinates
      double mag_p;
      double theta_p;
      double phi_p;
      reb_tools_xyz_to_spherical(srot, &mag_p, &theta_p, &phi_p);

      fprintf(of, "%e,%e,%e,%e,%e,%e,%e,%e,%e,%e,%e,%e,%e,%e,%e,%e,%e\n", sim->t, Omega_sun->x, Omega_sun->y, Omega_sun->z, mag_p, theta_p, phi_p, a1, e1, i1, Om1, pom1, a2, e2, i2, Om2, pom2); // print spins and orbits

      fclose(of);
    }

    if(reb_simulation_output_check(sim, 20.*M_PI)){        // outputs to the screen
        reb_simulation_output_timing(sim, tmax);
    }
}

This example is located in the directory examples/tides_spin_kozai

5.13. Obliquity Sculpting of Kepler Multis (Millholland & Laughlin 2019)

In particular, this simulates the obliquity sculpting of the Kepler multis due to convergent migration over 4 Myr. Two planets are initialized just wide of the 3:2 MMR (0.173 and 0.233 AU) are migrated inward. After 2 Myr, migration is turned off. After a period of chaotic obliquity evolution, the inner planet is excited to high obliquity and maintained indefinitely Result is based on Figure 3 in Millholland & Laughlin (2019), and reproduces Figure 3 in Lu et. al (2023) For a more in-depth description of the various parameters that can be set in this simulation, please see the ipython examples for consistent tides & spin (any notebook with the prefix TidesSpin) and migration (Migration.ipynb)

integration time is artificially shortened to run quickly. Full run in paper with tmax = 4e6 orbits takes ~ 2 days

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

void heartbeat(struct reb_simulation* sim);
double tmax = 100 * 2 * M_PI; // set short to run quickly. Set to 4e6 * 2 * M_PI in paper

int main(int argc, char* argv[]){
    struct reb_simulation* sim = reb_simulation_create();
    // Exact parameters from Millholland & Laughlin (2019)
    const double solar_mass = 1.;
    const double solar_rad = 0.00465;
    reb_simulation_add_fmt(sim, "m r", solar_mass, solar_rad);// Central object

    const double p1_mass = 5. * 3.0e-6; // in Earth masses * 1 Earth Mass / 1 Solar Mass
    const double p1_rad = 2.5 * 4.26e-5; // in Earth rad * 1 Earth rad / 1 AU
    reb_simulation_add_fmt(sim, "m a e r inc Omega pomega M", p1_mass, 0.17308688, 0.01, p1_rad, 0.5 * (M_PI / 180.), 0.0 * (M_PI / 180.), 0.0 * (M_PI / 180.), 0.0 * (M_PI / 180.)); // Planet 1

    const double p2_mass = 5. * 3.0e-6;
    const double p2_rad = 2.5 * 4.26e-5;
    reb_simulation_add_fmt(sim, "m a e r inc Omega pomega M", p2_mass, 0.23290608, 0.01, p2_rad, -0.431 * (M_PI / 180.), 0.0 * (M_PI / 180.), 0.0 * (M_PI / 180.), 0.0 * (M_PI / 180.)); // Planet 2
    sim->N_active = 3;
    sim->integrator = REB_INTEGRATOR_WHFAST;
    sim->dt = 1e-3;
    sim->heartbeat = heartbeat;

    // Add REBOUNDx Additional effects
    // First Spin
    struct rebx_extras* rebx = rebx_attach(sim);

    struct rebx_force* effect = rebx_load_force(rebx, "tides_spin");
    rebx_add_force(rebx, effect);
    // Exact parameters from Millholland & Laughlin (2019)
    // Sun
    const double solar_spin_period = 20 * 2 * M_PI / 365;
    const double solar_spin = (2 * M_PI) / solar_spin_period;
    const double solar_Q = 1000000.;
    rebx_set_param_double(rebx, &sim->particles[0].ap, "k2", 0.07);
    rebx_set_param_double(rebx, &sim->particles[0].ap, "I", 0.07 * solar_mass * solar_rad * solar_rad);
    rebx_set_param_vec3d(rebx, &sim->particles[0].ap, "Omega", (struct reb_vec3d){.z = solar_spin}); // Omega_x = Omega_y = 0 by default

    // We assume tau = 1/(2*n*Q) with n the mean motion, even though the spin is not synchronized with the orbit (see Lu et al. (2023))
    struct reb_orbit orb = reb_orbit_from_particle(sim->G, sim->particles[1], sim->particles[0]);
    rebx_set_param_double(rebx, &sim->particles[0].ap, "tau", 1./(2.*orb.n*solar_Q));

    // P1
    const double spin_period_1 = 5. * 2. * M_PI / 365.; // 5 days in REBOUND time units
    const double spin_1 = (2. * M_PI) / spin_period_1;
    const double planet_Q = 10000.;
    rebx_set_param_double(rebx, &sim->particles[1].ap, "k2", 0.4);
    rebx_set_param_double(rebx, &sim->particles[1].ap, "I", 0.25 * p1_mass * p1_rad * p1_rad);
    rebx_set_param_vec3d(rebx, &sim->particles[1].ap, "Omega", (struct reb_vec3d){.y=spin_1 * -0.0261769, .z=spin_1 * 0.99965732});

    rebx_set_param_double(rebx, &sim->particles[1].ap, "tau", 1./(2.*orb.n*planet_Q));

    // P2
    double spin_period_2 = 3. * 2. * M_PI / 365.; // 3 days in REBOUND time units
    double spin_2 = (2. * M_PI) / spin_period_2;
    rebx_set_param_double(rebx, &sim->particles[2].ap, "k2", 0.4);
    rebx_set_param_double(rebx, &sim->particles[2].ap, "I", 0.25 * p2_mass * p2_rad * p2_rad);
    rebx_set_param_vec3d(rebx, &sim->particles[2].ap, "Omega", (struct reb_vec3d){.y=spin_2 * 0.0249736, .z=spin_2 * 0.99968811});

    struct reb_orbit orb2 = reb_orbit_from_particle(sim->G, sim->particles[2], sim->particles[0]);
    rebx_set_param_double(rebx, &sim->particles[2].ap, "tau", 1./(2.*orb2.n*planet_Q));

    // And migration
    struct rebx_force* mo = rebx_load_force(rebx, "modify_orbits_forces");
    rebx_add_force(rebx, mo);

    // Set migration parameters
    rebx_set_param_double(rebx, &sim->particles[1].ap, "tau_a", -5e6 * 2 * M_PI);
    rebx_set_param_double(rebx, &sim->particles[2].ap, "tau_a", (-5e6 * 2 * M_PI) / 1.1);

    reb_simulation_move_to_com(sim);

    // Let's create a reb_rotation object that rotates to new axes with newz pointing along the total ang. momentum, and x along the line of
    // nodes with the invariable plane (along z cross newz)
    struct reb_vec3d newz = reb_vec3d_add(reb_simulation_angular_momentum(sim), rebx_tools_spin_angular_momentum(rebx));
    struct reb_vec3d newx = reb_vec3d_cross((struct reb_vec3d){.z =1}, newz);
    struct reb_rotation rot = reb_rotation_init_to_new_axes(newz, newx);
    rebx_simulation_irotate(rebx, rot); // This rotates our simulation into the invariable plane aligned with the total ang. momentum (including spin)
    rebx_spin_initialize_ode(rebx, effect);

    // Run simulation
    system("rm -v output_orbits.txt"); // remove previous output files
    system("rm -v output_spins.txt");

    reb_simulation_integrate(sim, tmax/2);

    printf("Migration Switching Off\n");
    rebx_set_param_double(rebx, &sim->particles[1].ap, "tau_a", INFINITY);
    rebx_set_param_double(rebx, &sim->particles[2].ap, "tau_a", INFINITY);

    reb_simulation_integrate(sim, tmax);

    rebx_free(rebx);
    reb_simulation_free(sim);
}

void heartbeat(struct reb_simulation* sim){
  if(reb_simulation_output_check(sim, tmax/100000)){        // outputs every 100 REBOUND years
    struct rebx_extras* const rebx = sim->extras;
    FILE* of_orb = fopen("output_orbits.txt", "a");
    FILE* of_spins = fopen("output_spins.txt", "a");
    if (of_orb == NULL || of_spins == NULL){
        reb_simulation_error(sim, "Can not open file.");
        return;
    }

    struct reb_particle* sun = &sim->particles[0];
    struct reb_particle* p1 = &sim->particles[1];
    struct reb_particle* p2 = &sim->particles[2];

    // Orbit information
    struct reb_orbit o1 = reb_orbit_from_particle(sim->G, *p1, *sun);
    double a1 = o1.a;
    double e1 = o1.e;
    double i1 = o1.inc;
    double Om1 = o1.Omega;
    double pom1 = o1.pomega;
    struct reb_vec3d norm1 = o1.hvec;

    struct reb_orbit o2 = reb_orbit_from_particle(sim->G, *p2, *sun);
    double a2 = o2.a;
    double e2 = o2.e;
    double i2 = o2.inc;
    double Om2 = o2.Omega;
    double pom2 = o2.pomega;
    struct reb_vec3d norm2 = o2.hvec;

    // Spin vectors - all initially in invariant plane
    struct reb_vec3d* Omega_sun = rebx_get_param(rebx, sun->ap, "Omega");

    // Interpret both planet spin vectors in the rotating planet frame in spherical coordinates
    struct reb_vec3d* Omega_p1 = rebx_get_param(rebx, p1->ap, "Omega");

    struct reb_vec3d lon1 = reb_vec3d_cross((struct reb_vec3d){.z =1}, norm1);  // Line of nodes is the new x-axis
    struct reb_rotation rot1 = reb_rotation_init_to_new_axes(norm1, lon1);      // Arguments to this function are the new z and x axes
    struct reb_vec3d sv1 = reb_vec3d_rotate(*Omega_p1, rot1);

    double mag1;
    double theta1;
    double phi1;
    reb_tools_xyz_to_spherical(sv1, &mag1, &theta1, &phi1);

    struct reb_vec3d* Omega_p2 = rebx_get_param(rebx, p2->ap, "Omega");

    struct reb_vec3d lon2 = reb_vec3d_cross((struct reb_vec3d){.z =1}, norm2); // Line of nodes is the new x-axis
    struct reb_rotation rot2 = reb_rotation_init_to_new_axes(norm2, lon2); // Arguments to this function are the new z and x axes
    struct reb_vec3d sv2 = reb_vec3d_rotate(*Omega_p2, rot2);

    double mag2;
    double theta2;
    double phi2;
    reb_tools_xyz_to_spherical(sv2, &mag2, &theta2, &phi2);

    fprintf(of_orb, "%e,%e,%e,%e,%e,%e,%e,%e,%e,%e,%e,%e,%e,%e,%e,%e,%e\n", sim->t, a1, e1, i1, pom1, Om1, norm1.x, norm1.y, norm1.z, a2, e2, i2, pom2, Om2, norm2.x, norm2.y, norm2.z);  // prints the spins and orbits of all bodies
    fclose(of_orb);
    fprintf(of_spins, "%e,%e,%e,%e,%e,%e,%e,%e,%e,%e,%e,%e,%e,%e,%e,%e\n", sim->t, Omega_sun->x, Omega_sun->y, Omega_sun->z, mag1, theta1, phi1, mag2, theta2, phi2, Omega_p1->x, Omega_p1->y, Omega_p1->z, Omega_p2->x, Omega_p2->y, Omega_p2->z);
    fclose(of_spins);
  }

  if(reb_simulation_output_check(sim, 100.*M_PI)){        // outputs to the screen
      reb_simulation_output_timing(sim, tmax);
  }
}

This example is located in the directory examples/tides_spin_migration_driven_obliquity_tides

5.14. Saving and loading simulations

This example demonstrates how to restart a simulation with all REBOUNDx effects and parameters.

#include <stdio.h>
#include <string.h>
#include "rebound.h"
#include "reboundx.h"
#include "core.h"

int main(int argc, char* argv[]){
    struct reb_simulation* sim = reb_simulation_create(); // make a simple sim with star and 1 planet
    struct reb_particle p = {0};
    p.m     = 1.;
    reb_simulation_add(sim, p);
    double m = 0.;
    double a = 1.e-4;
    double e = 0.2;
    double inc = 0.;
    double Omega = 0.;
    double omega = 0.;
    double f = 0.;

    struct reb_particle p1 = reb_particle_from_orbit(sim->G, p, m, a, e, inc, Omega, omega, f);
    reb_simulation_add(sim, p1);
    sim->dt = 1.e-8;
    sim->integrator = REB_INTEGRATOR_WHFAST;

    struct rebx_extras* rebx = rebx_attach(sim);

    // Add some (arbitrary) parameters to effects and particles

    struct rebx_force* gr = rebx_load_force(rebx, "gr");
    rebx_set_param_double(rebx, &gr->ap, "c", 3e4);
    rebx_set_param_int(rebx, &gr->ap, "gr_source", 0);
    rebx_add_force(rebx, gr);

    struct rebx_operator* mm = rebx_load_operator(rebx, "modify_mass");
    rebx_set_param_double(rebx, &sim->particles[0].ap, "tau_mass", -1.e4);
    rebx_add_operator(rebx, mm);

    struct reb_particle* pptr = &sim->particles[1];
    rebx_set_param_double(rebx, &pptr->ap, "tau_e", 1.5);   // no meaning, just to test storage
    rebx_set_param_int(rebx, &pptr->ap, "max_iterations", 42);

    double* c = rebx_get_param(rebx, gr->ap, "c");
    int* gr_source = rebx_get_param(rebx, gr->ap, "gr_source");
    double* tau_mass = rebx_get_param(rebx, sim->particles[0].ap, "tau_mass");
    double* tau_e = rebx_get_param(rebx, sim->particles[1].ap, "tau_e");
    int* max_iterations = rebx_get_param(rebx, sim->particles[1].ap, "max_iterations");

    struct rebx_operator* integforce = rebx_load_operator(rebx, "integrate_force");
    rebx_set_param_pointer(rebx, &integforce->ap, "force", gr);
    rebx_add_operator(rebx, integforce);

    printf("c: Original = %f\n", *c);
    printf("gr_source: Original = %d\n", *gr_source);
    printf("tau_mass: Original = %f\n", *tau_mass);
    printf("tau_e: Original = %f\n", *tau_e);
    printf("max_iterations: Original = %d\n", *max_iterations);

    // We now have to save both a REBOUND binary (for the simulation) and a REBOUNDx one (for parameters and effects)
    reb_simulation_integrate(sim, 1.e-4);
    reb_simulation_save_to_file(sim, "reb.bin");
    rebx_output_binary(rebx, "rebx.bin");

    rebx_free(rebx);
    reb_simulation_free(sim);

    // We now reload the simulation and the rebx instance (which adds previously loaded effects to the simulation)
    sim = reb_simulation_create_from_file("reb.bin", 0);
    rebx = rebx_create_extras_from_binary(sim, "rebx.bin");

    gr = rebx_get_force(rebx, "gr");
    mm = rebx_get_operator(rebx, "modify_mass");

    c = rebx_get_param(rebx, gr->ap, "c");
    gr_source = rebx_get_param(rebx, gr->ap, "gr_source");
    tau_mass = rebx_get_param(rebx, sim->particles[0].ap, "tau_mass");
    tau_e = rebx_get_param(rebx, sim->particles[1].ap, "tau_e");
    max_iterations = rebx_get_param(rebx, sim->particles[1].ap, "max_iterations");

    printf("c: Loaded = %f\n", *c);
    printf("gr_source: Loaded = %d\n", *gr_source);
    printf("tau_mass: Loaded = %f\n", *tau_mass);
    printf("tau_e: Loaded = %f\n", *tau_e);
    printf("max_iterations: Loaded = %d\n", *max_iterations);

    // You would now integrate as usual
    double tmax = 1.e-4;
    reb_simulation_integrate(sim, tmax);
    rebx_free(rebx);
    reb_simulation_free(sim);
}

This example is located in the directory examples/saving_and_loading_simulations

5.15. Stellar evolution with interpolated mass data

This example shows how to change a particle’s mass by interpolating time-series data during a REBOUND simulation. If you have GLUT installed for visualization, press ‘w’ to see the orbits as wires. You can zoom out by holding shift, holding down the mouse and dragging. Press ‘c’ to better see migration/e-damping.

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

void heartbeat(struct reb_simulation* sim);
struct rebx_interpolator* stellarmass;
struct rebx_extras* rebx;
double tmax = 1e4;
int main(int argc, char* argv[]) {
     struct reb_simulation* sim = reb_simulation_create();
    sim->G = 4*M_PI*M_PI; // use units of AU, yr and solar masses
     sim->heartbeat = heartbeat;
     sim->integrator = REB_INTEGRATOR_WHFAST;

     struct reb_particle sun = {0};
     sun.m   = 1.;
     reb_simulation_add(sim, sun);
     // Initialize planets on circular orbits, each 2 times farther than last.
     struct reb_particle planet = {0};
     planet.x = 1.;
     planet.vy = 2.*M_PI;
     reb_simulation_add(sim, planet);
     planet.x *= 2.;
     planet.vy /= sqrt(2.);
     reb_simulation_add(sim, planet);
     planet.x *= 2.;
     planet.vy /= sqrt(2.);
     reb_simulation_add(sim, planet);

     reb_simulation_move_to_com(sim);
     rebx = rebx_attach(sim); // initialize reboundx

     // To set how the mass will change, we pass three equally-sized arrays
     // necessary to interpolate the time-mass values. Here we have six (6)
     // values that correspond to a star losing mass with an e-damping timescale
     // of tmax (1e4) up to 12,500 yr. The effect will use a cubic spline to
     // interpolate any intermediate values needed by the simulation.
    int n = 6; // size of arrays
     double times[] = {0, 2500, 5000, 7500, 10000, 12500}; // in yr
     double values[] = {1., 0.77880078307, 0.60653065971, 0.47236655274, 0.36787944117, 0.28650479686}; // in Msun
    stellarmass = rebx_create_interpolator(rebx, n, times, values, REBX_INTERPOLATION_SPLINE);

    reb_simulation_integrate(sim, tmax);
    rebx_free_interpolator(stellarmass);
     rebx_free(rebx); // explicitly free all the memory allocated by REBOUNDx
}

void heartbeat(struct reb_simulation* sim) {
     if (reb_simulation_output_check(sim, tmax/1000)) {
        sim->particles[0].m = rebx_interpolate(rebx, stellarmass, sim->t);
        reb_simulation_move_to_com(sim);
        struct reb_orbit o = reb_orbit_from_particle(sim->G, sim->particles[1], sim->particles[0]);
        printf("t=%e, Sun mass = %f, planet semimajor axis = %f\n", sim->t, sim->particles[0].m, o.a);
     }
}

This example is located in the directory examples/parameter_interpolation

5.16. Adding gravitational harmonics (J2, J4) to particles

This example shows how to add a J2 and J4 harmonic to particles. If you have GLUT installed for the visualization, press ‘w’ and/or ‘c’ for a clearer view of the whole orbit.

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

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

    struct reb_particle star = {0};
    star.m     = 1.;
    star.hash  = reb_hash("star");
    reb_simulation_add(sim, star);

    double m = 0.;
    double a = 1.;
    double e = 0.2;
    double inc = 0.;
    double Omega = 0.;
    double omega = 0.;
    double f = 0.;

    struct reb_particle planet = reb_particle_from_orbit(sim->G, star, m, a, e, inc, Omega, omega, f);
    planet.hash = reb_hash("planet");
    reb_simulation_add(sim, planet);
    reb_simulation_move_to_com(sim);

    struct rebx_extras* rebx = rebx_attach(sim);
    struct rebx_force* gh = rebx_load_force(rebx, "gravitational_harmonics");
    rebx_add_force(rebx, gh);

    rebx_set_param_double(rebx, &sim->particles[0].ap, "J2", 0.1);
    rebx_set_param_double(rebx, &sim->particles[0].ap, "J4", 0.01);
    rebx_set_param_double(rebx, &sim->particles[0].ap, "R_eq", 0.01);

    double tmax = 1.e5;
    reb_simulation_integrate(sim, tmax);
    rebx_free(rebx);    // this explicitly frees all the memory allocated by REBOUNDx
    reb_simulation_free(sim);
}

This example is located in the directory examples/J2

5.17. Constant time lag model for tides (Hut 1981)

In particular, this simulates post-main sequence tidal interactions between the Earth and Sun near its tip-RGB phase. Definitely see the corresponding ipython example, as well as the documentation, for more explanations along the way of the various parameters and assumptions.

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

int main(int argc, char* argv[]){
    struct reb_simulation* sim = reb_simulation_create();
    sim->G = 4*M_PI*M_PI;           // Units of AU, yr and Msun

    struct reb_particle sun = {0};
    sun.m = 0.86;                   // post-MS in Msun
    reb_simulation_add(sim, sun);

    struct reb_orbit eo = {0};      // for Earth
    eo.a = 1.0;                     // in AU
    eo.e = 0.03;
    double e_mass = 2.988e-6;       // if planets don't have mass, they don't raise any tides
    struct reb_particle ep = reb_particle_from_orbit(sim->G, sun, e_mass, eo.a, eo.e, eo.inc, eo.Omega, eo.omega, eo.f);
    reb_simulation_add(sim, ep);
    reb_simulation_move_to_com(sim);

    // Add REBOUNDx Additional Effect
    struct rebx_extras* rebx = rebx_attach(sim);
    struct rebx_force* tides = rebx_load_force(rebx, "tides_constant_time_lag");
    rebx_add_force(rebx, tides);

    // We first have to give the body being tidally distorted a physical radius, or nothing will happen
    sim->particles[0].r = 0.85;     // AU

    // At a minimum, have to set tctl_k2 (potential Love number of degree 2) parameter, which will add the conservative potential of the equilibrium tidal distortion
    rebx_set_param_double(rebx, &sim->particles[0].ap, "tctl_k2", 0.03);

    // We add dissipation by adding a constant time lag tctl_tau
    rebx_set_param_double(rebx, &sim->particles[0].ap, "tctl_tau", 0.04);

    // We also can set the angular rotation rate of bodies OmegaMag. Implementation assumes Omega is along z axis
    // If you need a more general implementation, see the tides_spin effect. Implementation will assume zero if not specified.

    rebx_set_param_double(rebx, &sim->particles[0].ap, "OmegaMag", 2*M_PI/(30./365.)); // 30 day rotation rate
    //exit(1);
    // Run simulation
    double tmax = 2.5e5; // years
    reb_simulation_integrate(sim, tmax);
    rebx_free(rebx);
    reb_simulation_free(sim);
}

This example is located in the directory examples/tides_constant_time_lag

5.18. Stochastic forces on a single planet

This example shows how to add a stochastic force to a single planet. As a result, the planet will undergo a random walk in its orbital parameters. See also [Rein (2010)](https://ui.adsabs.harvard.edu/abs/2010PhDT…….196R/abstract) and [Rein and Papaloizou (2009)](https://ui.adsabs.harvard.edu/abs/2009A%26A…497..595R/abstract)

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

void heartbeat(struct reb_simulation* r);

double tmax = 1e4*M_PI*2.0; // 1e4 orbits

int main(int argc, char* argv[]){
    // Let's first setup the simulation.
    // Note that we are using the WHFast integrator with a fixed timestep.
    // It's important to point out that the IAS15 integrator is not well
    // suited for stochastic forces because it automatically reduces the
    // timestep if it doesn't achieve an accuracy near machine precision.
    // Because the stochastic forces are random, it might never converge.

    struct reb_simulation* sim = reb_simulation_create();
    sim->integrator     = REB_INTEGRATOR_WHFAST;
    sim->dt             = 1e-2;          // At ~100 AU, orbital periods are ~1000 yrs, so here we use ~1% of that, in sec
    sim->heartbeat      = heartbeat;

    reb_simulation_add_fmt(sim, "m", 1.); // Sun
    reb_simulation_add_fmt(sim, "m a", 1e-3, 1.0); // Jupiter mass planet at 1AU
    reb_simulation_move_to_com(sim);

    // Next, we add the `stochastic_forces` module in REBOUNDx
    struct rebx_extras* rebx = rebx_attach(sim);
    struct rebx_force* sto = rebx_load_force(rebx, "stochastic_forces");
    rebx_add_force(rebx, sto);

    // We can now turn on stochastic forces for a particle by setting
    // the field `kappa` to a finite value. This parameter determines
    // the strength of the stochastic forces and is relative to the
    // gravitational force that the particle experiences from the
    // center of mass interior to its orbit.
    rebx_set_param_double(rebx, &sim->particles[1].ap, "kappa", 1e-5);

    // The auto-correlation time of the stochastic forces is by default
    // the orbital period. We can set it to a fraction or multiple of
    // the orbital period by chanhging the `tau_kappa` parameter.

    // rebx_set_param_double(rebx, &sim->particles[1].ap, "tau_kappa", 1.); // Not needed. Already the default.


    // Let's integrate the system.
    reb_simulation_integrate(sim, tmax);

    rebx_free(rebx);
    reb_simulation_free(sim);
}

void heartbeat(struct reb_simulation* sim){
    // Periodically output the semi-major axis of the planet.
    if(reb_simulation_output_check(sim, 1.e2*M_PI*2.0)){
        reb_simulation_output_timing(sim, tmax);
        FILE* f = fopen("orbit.txt", "a+");
        struct reb_orbit o =  reb_orbit_from_particle(sim->G, sim->particles[1], sim->particles[0]);
        fprintf(f, "%.8e\t%.8e\n", sim->t, o.a);
        fclose(f);
    }
}

This example is located in the directory examples/stochastic_forces

5.19. Adding Lense-Thirring effect

This example shows how to add the Lense-Thirring effect to a simulation. If you have GLUT installed for the visualization, press ‘w’ and/or ‘c’ for a clearer view of the whole orbit.

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

int main(int argc, char* argv[]){
    struct reb_simulation* sim = reb_simulation_create();
    sim->G = 4*M_PI*M_PI; // units of AU, yr, Msun
    struct reb_particle star = {0};
    star.m     = 1.;
    reb_simulation_add(sim, star);
    double omega = 90.361036076; //solar rotation rate in rad/year
    double C_I = 0.06884; //solar moment of inertia prefactor
    double R_eq = 0.00465247264; //solar equatorial radius in AU

    struct reb_particle planet = {0};  // add a planet on a circular orbit (with default units where G=1)
    const double mp = 1.7e-7;   // approximate values for mercury in units of Msun and AU
    const double a = 0.39;
    const double e = 0.21;
    reb_simulation_add_fmt(sim, "m a e", mp, a, e);

    struct rebx_extras* rebx = rebx_attach(sim);  // first initialize rebx
    struct rebx_force* lense  = rebx_load_force(rebx, "lense_thirring"); // add our new force
    rebx_add_force(rebx, lense);
    rebx_set_param_vec3d(rebx, &sim->particles[0].ap, "Omega", (struct reb_vec3d){.x=0, .y=0, .z=omega});
    rebx_set_param_double(rebx, &sim->particles[0].ap, "I", C_I*star.m*R_eq*R_eq);

    // Have to set speed of light in right units (set by G & initial conditions).  Here we use default units of AU/(yr/2pi)
    rebx_set_param_double(rebx, &lense->ap, "lt_c", 63241.077); // speed of light in AU/yr

    double tmax = 1000.;
    reb_simulation_integrate(sim, tmax);
    struct reb_orbit orb = reb_orbit_from_particle(sim->G, sim->particles[1], sim->particles[0]);
    printf("Pericenter precession rate = %.3f arcsec/century\n", orb.pomega*206265/10.);
}

This example is located in the directory examples/lense_thirring

5.20. Pseudo-Synchronization (Hut 1981)

This example shows how to add quadrupole and tidal distortions to bodies with structure, letting us consistently track the spin-axis and dynamical evolution of the system. In particular, this simulates the pseudo-synchronization of a fiducial hot Jupiter. The hot Jupiter is initialized with a slight eccentricity, nontrivial obliquity and fast rotation. Under the influence of tidal dissipation, we see the following rapidly occur: circularization of the orbit, obliquity damping down to 0, and rotation settling to the pseudo-synchronous value predicted by Hut (1981) Definitely see the corresponding ipython example, as well as the documentation, for more in-depth explanations regarding the various parameters that may be set in this simulation. Also see Lu et. al (in review), Eggleton et. al (1998), Hut (1981).

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

void heartbeat(struct reb_simulation* sim);
double tmax = 1000 * 2 * M_PI;

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

    // Star
    const double solar_mass = 1.;
    const double solar_rad = 0.00465; // Bodies with structure require radius! This is one solar radius
    reb_simulation_add_fmt(sim, "m r", solar_mass, solar_rad);// Central object

    // Fiducial hot Jupiter
    const double p1_mass = 1. * 9.55e-4; // in Jupiter masses * 1 Jupiter Mass / 1 Solar Mass
    const double p1_rad = 1. * 4.676e-4; // in Jupiter rad * 1 jupiter rad / 1 AU
    const double p1_e = 0.01;
    const double p1_inc = 0.01;
    reb_simulation_add_fmt(sim, "m a e inc r", p1_mass, 0.04072, p1_e, p1_inc, p1_rad); // Planet 1

    sim->N_active = 2;
    sim->integrator = REB_INTEGRATOR_WHFAST;
    sim->dt = 1e-3;
    sim->heartbeat = heartbeat;

    // Add tides_spin as a REBOUNDx additional effect
    struct rebx_extras* rebx = rebx_attach(sim);
    struct rebx_force* effect = rebx_load_force(rebx, "tides_spin");
    rebx_add_force(rebx, effect);
    // Star
    // The following parameters are mandatory to set for the body to have structure:
    // The Love number (k2)
    // The three components of the anguarl spin frequency (Omega_x, Omega_y, Omega_z)
    const double solar_k2 = 0.07;
    rebx_set_param_double(rebx, &sim->particles[0].ap, "k2", solar_k2);
    const double solar_spin_period = 27 * 2 * M_PI / 365; // 27 Days in REBOUND time units
    const double solar_spin = (2 * M_PI) / solar_spin_period;
    rebx_set_param_vec3d(rebx, &sim->particles[0].ap, "Omega", (struct reb_vec3d){.z=solar_spin}); // Omega.x = Omega.y = 0 by default

    // While not technically necessary, the fully dimensional moment of inertia (I) should also be set
    // This is required to evolve the spin axis! Without setting this value, the spin axis will remain stationary.
    rebx_set_param_double(rebx, &sim->particles[0].ap, "I", 0.07 * solar_mass * solar_rad * solar_rad);

    // Finally, the last parameter is the constant time lag tau
    const double solar_Q = 1e6; // Often tidal dissipation is expressed in terms of a tidal quality factor Q
    struct reb_orbit orb = reb_orbit_from_particle(sim->G, sim->particles[1], sim->particles[0]);
    // In the case of a spin that is synchronous with a circular orbit, tau is related to the tidal quality factor Q through the orbital mean motion n
    double solar_tau = 1 / (2 * solar_Q * orb.n);
    // See Lu et. al (2023) for discussion of the cases in which this approximation relating Q and tau is valid

    rebx_set_param_double(rebx, &sim->particles[0].ap, "tau", solar_tau);

    // Planet - all of the above applies here too
    const double spin_period_1 = 0.5 * 2. * M_PI / 365.; // 0.5 days in reb years
    const double spin_1 = (2. * M_PI) / spin_period_1;
    const double planet_Q = 10000.;
    const double theta_1 = 30. * (M_PI / 180.);
    const double phi_1 = 0 * (M_PI / 180);
    rebx_set_param_double(rebx, &sim->particles[1].ap, "k2", 0.3);
    rebx_set_param_double(rebx, &sim->particles[1].ap, "I", 0.25 * p1_mass * p1_rad * p1_rad);

    // Let's consider a tilted planetary spin axis. The spin frequency vector is expressed in the same inertial reference frame as the simulation.
    // Given a polar angle theta (from the z axis) and azimuthal angle phi (from the x axis), we could either manually set spin axis components:
    struct reb_vec3d Omega_1;
    Omega_1.x = spin_1 * sin(theta_1) * cos(phi_1);
    Omega_1.y = spin_1 * sin(theta_1) * sin(phi_1);
    Omega_1.z = spin_1 * cos(theta_1);

    // Or use the built-in convenience function, which returns the Cartesian coordinates of the spin vector given:
    // magnitude, obliquity, and phase angle
    Omega_1 = reb_tools_spherical_to_xyz(spin_1, theta_1, phi_1);

    // For your own use case, you would just use whichever of the above two methods is most convenient (here we've done it three ways)
    // whichever of the two methods we use above, we need to remember to set the spin axis values
    rebx_set_param_vec3d(rebx, &sim->particles[1].ap, "Omega", Omega_1);
    rebx_set_param_double(rebx, &sim->particles[1].ap, "tau", 1./(2*planet_Q*orb.n));

    reb_simulation_move_to_com(sim);

    // Let's create a reb_rotation object that rotates our simulation to new axes with newz pointing along the total ang. momentum, and x along the line of
    // nodes with the invariable plane (along z cross newz)
    struct reb_vec3d newz = reb_vec3d_add(reb_simulation_angular_momentum(sim), rebx_tools_spin_angular_momentum(rebx));
    struct reb_vec3d newx = reb_vec3d_cross((struct reb_vec3d){.z =1}, newz);
    struct reb_rotation rot = reb_rotation_init_to_new_axes(newz, newx);
    rebx_simulation_irotate(rebx, rot); // This rotates our simulation into the invariable plane aligned with the total ang. momentum (including spin)
    rebx_spin_initialize_ode(rebx, effect); // We must call this function before starting our integration to track the spins

    system("rm -v output.txt"); // remove previous output file
    reb_simulation_integrate(sim, tmax);
    rebx_free(rebx);
    reb_simulation_free(sim);
}

void heartbeat(struct reb_simulation* sim){
    // Output spin and orbital information to file
    if(reb_simulation_output_check(sim, 10)){        // outputs every 10 REBOUND years
      struct rebx_extras* const rebx = sim->extras;
      FILE* of = fopen("output.txt", "a");
      if (of==NULL){
          reb_simulation_error(sim, "Can not open file.");
          return;
      }

      struct reb_particle* star = &sim->particles[0];
      struct reb_particle* p = &sim->particles[1];

      struct reb_orbit orb = reb_orbit_from_particle(sim->G, *p, *star);
      double a = orb.a;
      double Om = orb.Omega;
      double inc = orb.inc;
      double pom = orb.pomega;
      double e = orb.e;

      // The spin vector in the inertial (in this case, invariant frame)
      struct reb_vec3d* Omega_inv = rebx_get_param(rebx, p->ap, "Omega");

      // Transform spin vector into planet frame, w/ z-axis aligned with orbit normal and x-axis aligned with line of nodes
      struct reb_vec3d orbit_normal = orb.hvec;
      struct reb_vec3d line_of_nodes = reb_vec3d_cross((struct reb_vec3d){.z =1}, orbit_normal);
      struct reb_rotation rot = reb_rotation_init_to_new_axes(orbit_normal, line_of_nodes); // Arguments to this function are the new z and x axes
      struct reb_vec3d srot = reb_vec3d_rotate(*Omega_inv, rot); // spin vector in the planet's frame

      // Interpret the spin axis in the more natural spherical coordinates
      double mag;
      double theta;
      double phi;
      reb_tools_xyz_to_spherical(srot, &mag, &theta, &phi);

      fprintf(of, "%e,%e,%e,%e,%e,%e,%e,%e,%e\n", sim->t, a, inc, e, pom, Om, mag, theta, phi);
      fclose(of);
    }

    if(reb_simulation_output_check(sim, 20.*M_PI)){        // outputs to the screen
        reb_simulation_output_timing(sim, tmax);
    }
}

This example is located in the directory examples/tides_spin_pseudo_synchronization

5.21. Exponential mass loss/gain

This example shows how to add exponential mass loss to bodies in the integration. The primary’s mass loss causes planets’ orbits to expand. For more flexible variations in mass or other parameters see the parameter interpolation example.

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

void heartbeat(struct reb_simulation* sim);
double tmax = 1.e4;

int main(int argc, char* argv[]){
     struct reb_simulation* sim = reb_simulation_create();
    sim->G = 4*M_PI*M_PI;               // use units of AU, yr and solar masses
     sim->heartbeat = heartbeat;

     struct reb_particle sun = {0};
     sun.m   = 1.;
     reb_simulation_add(sim, sun);

     struct reb_particle planet = {0};   // Initialize planets on circular orbits, each 2 times farther than last.
     planet.m = 1.e-3;
     planet.x = 1.;
     planet.vy = 2.*M_PI;
     reb_simulation_add(sim,planet);
     planet.x *= 2.;
     planet.vy /= sqrt(2.);
     reb_simulation_add(sim, planet);
     planet.x *= 2.;
     planet.vy /= sqrt(2.);
     reb_simulation_add(sim, planet);

     reb_simulation_move_to_com(sim);

     struct rebx_extras* rebx = rebx_attach(sim); // initialize reboundx
    struct rebx_operator* modify_mass = rebx_load_operator(rebx, "modify_mass");

    /* The function rebx_add_operator will choose how to add the operator to the integration
     * scheme based on the integrator being used and the properties of the operator.
     * This is typically a half operator timestep before the main REBOUND timestep, and half afterward.
     */

     rebx_add_operator(rebx, modify_mass);

    /* If you wanted to make your own choices, you can add individual operator steps.
     * In this case you would pass additional parameters. Say we wanted to add a full operaator timestep after the main REBOUND timestep;
     *
     * dt_fraction = 1. // Fraction of a REBOUND timestep (sim->dt) operator should act
     * timing = REBX_TIMING_POST; // Should happen POST timestep
     * name = "modify_mass_post"; // Name identifier
     * rebx_add_operator_step(rebx, modify_mass, dt_fraction, timing, name);
     */

     // To set an exponential mass loss rate, we set the e-folding timescale (positive for growth, negative for loss)
    // Here have the star lose mass with e-damping timescale = tmax
    rebx_set_param_double(rebx, &sim->particles[0].ap, "tau_mass", -tmax);

     // We can approximate a linear mass loss/growth rate if the rate is small by taking tau_mass = M_initial / mass_loss_rate (or growth)
     double M_dot = 1.e-12; // mass growth rate for the planet (in simulation units--here Msun/yr)
    double tau_mass = sim->particles[1].m / M_dot; // first planet gains mass at linear rate M_dot
    rebx_set_param_double(rebx, &sim->particles[1].ap, "tau_mass", tau_mass);

     reb_simulation_integrate(sim, tmax);
     rebx_free(rebx);        // this explicitly frees all the memory allocated by REBOUNDx
}

void heartbeat(struct reb_simulation* const sim){
     // Output masses and semimajor of the inner planet 100 times over the time of the simulation
    if(reb_simulation_output_check(sim, tmax/100.)){
             struct reb_orbit o = reb_orbit_from_particle(sim->G, sim->particles[1], sim->particles[0]);
             printf("t=%e, Sun mass = %f, planet mass = %e, planet semimajor axis = %f\n", sim->t, sim->particles[0].m, sim->particles[1].m, o.a);
     }
 }

This example is located in the directory examples/modify_mass

5.22. General central force.

This example shows how to add a general central force. If you have GLUT installed for the visualization, press ‘w’ and/or ‘c’ for a clearer view of the whole orbit.

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

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

    struct reb_particle star = {0};
    star.m     = 1.;
    reb_simulation_add(sim, star);

    double m = 0.;
    double a = 1.; // put planet close to enhance precession so it's visible in visualization (this would put planet inside the Sun!)
    double e = 0.2;
    double inc = 0.;
    double Omega = 0.;
    double omega = 0.;
    double f = 0.;

    struct reb_particle planet = reb_particle_from_orbit(sim->G, star, m, a, e, inc, Omega, omega, f);
    reb_simulation_add(sim, planet);
    reb_simulation_move_to_com(sim);

    struct rebx_extras* rebx = rebx_attach(sim);
    struct rebx_force* force = rebx_load_force(rebx, "central_force");
    rebx_add_force(rebx, force);

    /* We first choose a power for our central force (here F goes as r^-1).
     * We then need to add it to the particle(s) that will act as central sources for this force.*/

    struct reb_particle* ps = sim->particles;
    double gammacentral = -1.;
    rebx_set_param_double(rebx, &ps[0].ap, "gammacentral", gammacentral);

    // The other parameter to set is the normalization Acentral (F=Acentral*r^gammacentral). E.g.,

    rebx_set_param_double(rebx, &ps[0].ap, "Acentral", 1.e-4);

    /* We can also use the function rebx_central_force_Acentral to calculate the Acentral required
     * for particles[1] (around primary particles[0]) to have a pericenter precession rate of
     * pomegadot, given a gammacentral value: */

    double pomegadot = 1.e-3;
    double Acentral = rebx_central_force_Acentral(ps[1], ps[0], pomegadot, gammacentral);
    rebx_set_param_double(rebx, &ps[0].ap, "Acentral", Acentral);

    double tmax = 3.e4;
    reb_simulation_integrate(sim, tmax);
    rebx_free(rebx);    // this explicitly frees all the memory allocated by REBOUNDx
    reb_simulation_free(sim);
}

This example is located in the directory examples/central_force

5.23. Migration and other orbit modifications

This example shows how to add migration, eccentricity damping and pericenter precession to a REBOUND simulation. If you have GLUT installed for visualization, press ‘w’ to see the orbits as wires. You can zoom out by holding shift, holding down the mouse and dragging. Press ‘c’ to better see migration/e-damping.

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

void heartbeat(struct reb_simulation* sim);

int main(int argc, char* argv[]){
    struct reb_simulation* sim = reb_simulation_create();
    // Setup constants
    sim->integrator = REB_INTEGRATOR_WHFAST;
    sim->dt             = 0.012;        // initial timestep.
    sim->heartbeat = heartbeat;

    struct reb_particle p = {0};
    p.m     = 1.;
    reb_simulation_add(sim, p);

    double m = 0.;
    double a1 = 1.;
    double a2 = 2.;
    double e = 0.4;
    double inc = 0.;
    double Omega = 0.;
    double omega = 0.;
    double f = 0.;

    struct reb_particle p1 = reb_particle_from_orbit(sim->G, p, m, a1, e, inc, Omega, omega, f);
    struct reb_particle p2 = reb_particle_from_orbit(sim->G, p, m, a2, e, inc, Omega, omega, f);
    reb_simulation_add(sim,p1);
    reb_simulation_add(sim,p2);
    reb_simulation_move_to_com(sim);

    struct rebx_extras* rebx = rebx_attach(sim);

    // There are two options for how to modify orbits.  You would only choose one (comment the other out).
    // You can't set precession separately with modify_orbits_forces (eccentricity and inclination damping induce pericenter and nodal precession).

    struct rebx_operator* mo = rebx_load_operator(rebx, "modify_orbits_direct");                                     // directly update particles' orbital elements each timestep
    rebx_add_operator(rebx, mo);
    //struct rebx_force* mo = rebx_load_force(rebx, "modify_orbits_forces");                                         // add forces that orbit-average to give exponential a and e damping
    //rebx_add_force(rebx, mo);

    // Set the timescales for each particle.
    double tmax = 5.e4;

     rebx_set_param_double(rebx, &sim->particles[1].ap, "tau_a", -tmax);         // add semimajor axis damping on inner planet (e-folding timescale)
    rebx_set_param_double(rebx, &sim->particles[1].ap, "tau_omega", -tmax/10.); // add linear precession (set precession period). Won't do anything for modify_orbits_forces
    rebx_set_param_double(rebx, &sim->particles[2].ap, "tau_e", -tmax/10.);     // add eccentricity damping on particles[2] (e-folding timescale)

    /* One can also adjust a coupling parameter between eccentricity and semimajor axis damping.  We use the parameter p
     * as defined by Deck & Batygin (2015).  The default p=0 corresponds to no coupling, while p=1 corresponds to e-damping
     * at constant angular momentum.  This is only implemented for modify_orbits_direct.
      * modify_orbits_forces damps eccentricity at constant angular momentum.
     *
     * Additionally, the damping by default is done in Jacobi coordinates.  If you'd prefer to use barycentric
     * coordinates, or coordinates referenced to a particular particle, set a coordinates parameter in the effect
     * parameters returned by rebx_add to REBX_COORDINATES_BARYCENTRIC or REBX_COORDINATES_PARTICLE.
     * If the latter, add a 'primary' flag to the reference particle (not neccesary for barycentric):
     */

     rebx_set_param_int(rebx, &mo->ap, "coordinates", REBX_COORDINATES_PARTICLE);
     rebx_set_param_double(rebx, &mo->ap, "p", 1.); // doesn't do anything for modify_orbits_forces
     rebx_set_param_int(rebx, &sim->particles[0].ap, "primary", 1);

    reb_simulation_integrate(sim, tmax);
    rebx_free(rebx);    // Free all the memory allocated by rebx
    reb_simulation_free(sim);
}

void heartbeat(struct reb_simulation* sim){
    // output a e and pomega (Omega + omega) of inner body
    if(reb_simulation_output_check(sim, 5.e2)){
        const struct reb_particle sun = sim->particles[0];
        const struct reb_orbit orbit = reb_orbit_from_particle(sim->G, sim->particles[1], sun); // calculate orbit of particles[1]
        printf("%f\t%f\t%f\t%f\n",sim->t,orbit.a, orbit.e, orbit.pomega);
    }
}

This example is located in the directory examples/modify_orbits

5.24. Radiation forces on circumplanetary dust

This example shows how to integrate circumplanetary dust particles under the action of radiation forces using IAS15. We use Saturn’s Phoebe ring as an example, a distant ring of debris, The output is custom, outputting the semi-major axis of every dust particle relative to the planet.

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

void heartbeat(struct reb_simulation* r);

double tmax = 1e10;

int main(int argc, char* argv[]){
    struct reb_simulation* sim = reb_simulation_create();
    // Setup constants
    sim->integrator     = REB_INTEGRATOR_IAS15;
    sim->G              = 6.674e-11;    // Use SI units
    sim->dt             = 1e4;          // Initial timestep in sec
    sim->N_active       = 2;            // Only the sun and the planet affect other particles gravitationally
    sim->heartbeat      = heartbeat;
    sim->usleep     = 1000;             // Slow down integration (for visualization only)

    // sun
    struct reb_particle sun = {0};
    sun.m  = 1.99e30;                   // mass of Sun in kg
    reb_simulation_add(sim, sun);

    // Saturn (simulation set up in Saturn's orbital plane, i.e., inc=0, so only need 4 orbital elements
    double mass_sat = 5.68e26;          // Mass of Saturn
    double a_sat = 1.43e12;             // Semimajor axis of Saturn in m
    double e_sat = 0.056;               // Eccentricity of Saturn
    double pomega_sat = 0.;             // Angle from x axis to pericenter
    double f_sat = 0.;                  // True anomaly of Saturn
    double inc = 0.;
    double Omega = 0.;
    struct reb_particle saturn = reb_particle_from_orbit(sim->G, sun, mass_sat, a_sat, e_sat, inc, Omega, pomega_sat, f_sat);

    reb_simulation_add(sim, saturn);

    // Add REBOUNDx
    struct rebx_extras* rebx = rebx_attach(sim);
    struct rebx_force* rad = rebx_load_force(rebx, "radiation_forces");
    double c = 3.e8;                    // speed of light in SI units
    rebx_set_param_double(rebx, &rad->ap, "c", c);

    // Will assume particles[0] is the radiation source by default. You can also add a flag to a particle explicitly
    rebx_set_param_int(rebx, &sim->particles[0].ap, "radiation_source", 1);

    /* Dust particles
     Here we imagine particles launched from Saturn's irregular Satellite Phoebe.
     Such grains will inherit the moon's orbital elements (e.g. Tamayo et al. 2011)

     In order for a particle to feel radiation forces, we have to set their beta parameter,
     the ratio of the radiation pressure force to the gravitional force from the star (Burns et al. 1979).
     We do this in two ways below.*/

    double a_dust = 1.30e10;            // semimajor axis of satellite Phoebe, in m
    double e_dust = 0.16;               // eccentricity of Phoebe
    double inc_dust = 175.*M_PI/180.;   // inclination of Phoebe to Saturn's orbital plane
    double Omega_dust = 0.;             // longitude of ascending node
    double omega_dust = 0.;             // argument of pericenter
    double f_dust = 0.;                 // true anomaly

    // We first set up the orbit and add the particles
    double m_dust = 0.;                 // treat dust particles as massless
    struct reb_particle p = reb_particle_from_orbit(sim->G, sim->particles[1], m_dust, a_dust, e_dust, inc_dust, Omega_dust, omega_dust, f_dust);
    reb_simulation_add(sim, p);

    // For the first particle we simply specify beta directly.
    rebx_set_param_double(rebx, &sim->particles[2].ap, "beta", 0.1);

    // We now add a 2nd particle on the same orbit, but set its beta using physical parameters.
    struct reb_particle p2 = reb_particle_from_orbit(sim->G, sim->particles[1], 0., a_dust, e_dust, inc_dust, Omega_dust, omega_dust, f_dust);
    reb_simulation_add(sim, p2);

    /* REBOUNDx has a convenience function to calculate beta given the gravitational constant G, the star's luminosity and mass, and the grain's physical radius, density and radiation pressure coefficient Q_pr (Burns et al. 1979). */

    // Particle parameters
    double radius = 1.e-5;              // in meters
    double density = 1.e3;              // kg/m3 = 1g/cc
    double Q_pr = 1.;                   // Equals 1 in limit where particle radius >> wavelength of radiation
    double L = 3.85e26;                 // Luminosity of the sun in Watts

    double beta = rebx_rad_calc_beta(sim->G, c, sim->particles[0].m, L, radius, density, Q_pr);
    rebx_set_param_double(rebx, &sim->particles[3].ap, "beta", beta);

    printf("Particle 2 has beta = %f\n", beta);

    reb_simulation_move_to_com(sim);

    printf("Time\t\tEcc (p)\t\tEcc (p2)\n");
    reb_simulation_integrate(sim, tmax);
    rebx_free(rebx);                /* free memory allocated by REBOUNDx */
}

void heartbeat(struct reb_simulation* sim){
    if(reb_simulation_output_check(sim, 1.e8)){
        const struct reb_particle* particles = sim->particles;
        const struct reb_particle saturn = particles[1];
        const double t = sim->t;
        struct reb_orbit orbit = reb_orbit_from_particle(sim->G, sim->particles[2], saturn); /* calculate orbit of particles[2] around Saturn */
        double e2 = orbit.e;
        orbit = reb_orbit_from_particle(sim->G, sim->particles[3], saturn);
        double e3 = orbit.e;
        printf("%e\t%f\t%f\n", t, e2, e3);
    }
}

This example is located in the directory examples/rad_forces_circumplanetary