Radiation forces on circumplanetary dust (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 shows how to integrate circumplanetary dust particles using the IAS15 integrator. The example sets the function pointer additional_forces
to a function that describes the radiation forces. The example uses a beta parameter of 0.01. The output is custom too, outputting the semi-major axis of every dust particle relative to the planet.
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "rebound.h"
void force_radiation(struct reb_simulation* r);
void heartbeat(struct reb_simulation* r);
double betaparticles = 0.01; // Beta parameter, defined as the ratio of radiation pressure over gravity.
double tmax = 1e6;
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-4; // Initial timestep.
r->N_active = 2; // Only the star and the planet are massive.
r->additional_forces = force_radiation;
r->heartbeat = heartbeat;
r->usleep = 5000; // Slow down integration (for visualization only)
// Star
struct reb_particle star = {0};
star.m = 1.;
reb_simulation_add(r, star);
// planet
struct reb_particle planet = {0};
planet.m = 1e-3;
planet.x = 1;
planet.vy = sqrt(r->G*(star.m+planet.m)/planet.x);
reb_simulation_add(r, planet);
// Dust particles
while(r->N<3){ // Three particles in total (star, planet, dust particle)
struct reb_particle p = {0};
p.m = 0; // massless
double _r = 0.001; // distance from planet planet
double v = sqrt(r->G*planet.m/_r);
p.x = _r;
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");
reb_simulation_integrate(r, tmax);
}
void force_radiation(struct reb_simulation* r){
struct reb_particle* particles = r->particles;
const struct reb_particle star = particles[0]; // cache
const int N = r->N;
const double G = r->G;
#pragma omp parallel for
for (int i=0;i<N;i++){
const struct reb_particle p = particles[i]; // cache
if (p.m!=0.) continue; // Only dust particles feel radiation forces
const double prx = p.x-star.x;
const double pry = p.y-star.y;
const double prz = p.z-star.z;
const double pr = sqrt(prx*prx + pry*pry + prz*prz); // distance relative to star
const double prvx = p.vx-star.vx;
const double prvy = p.vy-star.vy;
const double prvz = p.vz-star.vz;
const double c = 1.006491504759635e+04; // speed of light.
const double rdot = (prvx*prx + prvy*pry + prvz*prz)/pr; // radial velocity relative to star
const double F_r = betaparticles*G*star.m/(pr*pr);
// Equation (5) of Burns, Lamy, Soter (1979)
particles[i].ax += F_r*((1.-rdot/c)*prx/pr - prvx/c);
particles[i].ay += F_r*((1.-rdot/c)*pry/pr - prvy/c);
particles[i].az += F_r*((1.-rdot/c)*prz/pr - prvz/c);
}
}
void heartbeat(struct reb_simulation* r){
if(reb_simulation_output_check(r, M_PI*2.)){
reb_simulation_output_timing(r, tmax);
}
if(reb_simulation_output_check(r, M_PI*2.)){ // output every year
FILE* f = fopen("a.txt","ab");
const struct reb_particle* particles = r->particles;
const struct reb_particle planet = particles[1];
const double G = r->G;
const double t = r->t;
const int N = r->N;
for (int i=2;i<N;i++){
const struct reb_particle p = particles[i];
const double prx = p.x-planet.x;
const double pry = p.y-planet.y;
const double prz = p.z-planet.z;
const double pr = sqrt(prx*prx + pry*pry + prz*prz); // distance relative to star
const double pvx = p.vx-planet.vx;
const double pvy = p.vy-planet.vy;
const double pvz = p.vz-planet.vz;
const double pv = sqrt(pvx*pvx + pvy*pvy + pvz*pvz); // distance relative to star
const double a = -G*planet.m/( pv*pv - 2.*G*planet.m/pr ); // semi major axis
fprintf(f,"%e\t%e\n",t,a);
}
fclose(f);
}
}
This example is located in the directory examples/circumplanetarydust