J2 precession (C)
Try it out this example!
REBOUND has been compiled with emscripten to WebAssembly. This lets you run this example interactively from within your browser at almost native speed. No installation is required. Click here to try it out.
This example presents an implementation of the J2 gravitational moment. The equation of motions are integrated with the 15th order IAS15 integrator. The parameters in this example have been chosen to represent those of Saturn, but one can easily change them or even include higher order terms in the multipole expansion.
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "rebound.h"
const double J2planet = 16298e-6; // J2 of Saturn (Murray and Dermott p 531)
const double Mplanet = 0.00028588598; // mass of Saturn in solar masses
const double Rplanet = 0.00038925688; // radius of Saturn in AU
const double ObliquityPlanet = 0.; // obliquity of the planet
const double tmax = 1e3; // Maximum integration time
void heartbeat(struct reb_simulation* r);
void force_J2(struct reb_simulation* r);
int main(int argc, char* argv[]){
struct reb_simulation* r = reb_simulation_create();
// Start the visualization web server.
// Point your browser to http://localhost:1234
reb_simulation_start_server(r, 1234);
// Setup constants
r->integrator = REB_INTEGRATOR_IAS15;
r->dt = 1e-6; // initial timestep
r->N_active = 2; // only the star and the planet are massive.
// Planet
struct reb_particle planet = {0};
planet.m = Mplanet;
reb_simulation_add(r, planet);
struct reb_particle p = {0}; // test particle
double a = Rplanet*3.; // small distance from planet (makes J2 important)
double e = 0.1;
double v = sqrt((1.+e)/(1.-e)*r->G*planet.m/a); // setup eccentric orbit (ignores J2)
p.x = (1.-e)*a;
p.vy = v;
p.x += planet.x; p.y += planet.y; p.z += planet.z;
p.vx += planet.vx; p.vy += planet.vy; p.vz += planet.vz;
reb_simulation_add(r, p);
reb_simulation_move_to_com(r);
remove("a.txt"); // delete previous output
// Setup callback functions
r->heartbeat = heartbeat;
r->additional_forces = force_J2;
reb_simulation_integrate(r, tmax);
reb_simulation_free(r);
}
void force_J2(struct reb_simulation* r){
if (J2planet==0) return;
// Star
const struct reb_particle planet = r->particles[0]; // cache
const int N = r->N;
#pragma omp parallel for
for (int i=1;i<N;i++){
const struct reb_particle p = r->particles[i]; // cache
const double sprx = p.x-planet.x;
const double spry = p.y-planet.y;
const double sprz = p.z-planet.z;
const double prx = sprx*cos(-ObliquityPlanet) + sprz*sin(-ObliquityPlanet);
const double pry = spry;
const double prz =-sprx*sin(-ObliquityPlanet) + sprz*cos(-ObliquityPlanet);
const double pr2 = prx*prx + pry*pry + prz*prz; // distance^2 relative to planet
const double fac = 3.*r->G*J2planet*planet.m*Rplanet*Rplanet/2./pow(pr2,3.5);
const double pax = fac*prx*(prx*prx + pry*pry - 4.*prz*prz);
const double pay = fac*pry*(prx*prx + pry*pry - 4.*prz*prz);
const double paz = fac*prz*(3.*(prx*prx + pry*pry) - 2.*prz*prz);
r->particles[i].ax += pax*cos(ObliquityPlanet) + paz*sin(ObliquityPlanet);
r->particles[i].ay += pay;
r->particles[i].az +=-pax*sin(ObliquityPlanet) + paz*cos(ObliquityPlanet);
}
}
void heartbeat(struct reb_simulation* r){
if(reb_simulation_output_check(r, 4000.*r->dt)){ // output something to screen
reb_simulation_output_timing(r, tmax);
}
if(reb_simulation_output_check(r,M_PI*2.*0.01)){ // output some orbital parameters to file
FILE* f = fopen("a.txt","ab");
const struct reb_particle planet = r->particles[0];
const int N = r->N;
for (int i=1;i<N;i++){
struct reb_orbit o = reb_orbit_from_particle(r->G, r->particles[i],planet);
fprintf(f,"%.15e\t%.15e\t%.15e\t%.15e\n",r->t,o.a,o.e,o.omega);
}
fclose(f);
}
}
This example is located in the directory examples/J2