Radiation forces (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 provides an implementation of the Poynting-Robertson effect. The code is using the IAS15 integrator which is ideally suited for this velocity dependent force.
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "rebound.h"
void radiation_forces(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 = 1e5;
int main(int argc, char* argv[]){
struct reb_simulation* r = reb_simulation_create();
// Start the REBOUND visualization server. This
// allows you to visualize the simulation by pointing
// your web browser to http://localhost:1234
reb_simulation_start_server(r, 1234);
// setup constants
r->dt = 1e-3; // initial timestep
r->integrator = REB_INTEGRATOR_IAS15;
r->ri_ias15.epsilon = 1e-4; // accuracy parameter
r->N_active = 1; // the star is the only massive particle
r->force_is_velocity_dependent = 1;
r->additional_forces = radiation_forces; // setup callback function for velocity dependent forces
r->heartbeat = heartbeat;
// star is at rest at origin
struct reb_particle star = {0};
star.m = 1.;
reb_simulation_add(r, star);
// dust particles are initially on a circular orbit
while(r->N<2){
struct reb_particle p = {0};
p.m = 0; // massless
double a = 1.; // a = 1 AU
double v = sqrt(r->G*(star.m*(1.-betaparticles))/a);
double phi = reb_random_uniform(r, 0,2.*M_PI); // random phase
p.x = a*sin(phi); p.y = a*cos(phi);
p.vx = -v*cos(phi); p.vy = v*sin(phi);
reb_simulation_add(r, p);
}
remove("radius.txt"); // remove previous output
reb_simulation_integrate(r, tmax);
}
void radiation_forces(struct reb_simulation* r){
struct reb_particle* particles = r->particles;
const int N = r->N;
const struct reb_particle star = particles[0]; // cache
#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 in unit of G=1, M_sun=1, 1year=1
const double rdot = (prvx*prx + prvy*pry + prvz*prz)/pr; // radial velocity relative to star
const double F_r = betaparticles*r->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, 400.)){ // print some information to screen
reb_simulation_output_timing(r, tmax);;
}
if(reb_simulation_output_check(r, M_PI*2.*1000.)){ // output radial distance every 1000 years
FILE* f = fopen("radius.txt","ab");
struct reb_particle* particles = r->particles;
const struct reb_particle star = particles[0];
const int N = r->N;
for (int i=1;i<N;i++){
const struct reb_particle p = particles[i];
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
fprintf(f,"%e\t%e\n",r->t,pr);
}
fclose(f);
}
}
This example is located in the directory examples/prdrag