# Radiation forces on circumplanetary dust (C)

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 <unistd.h>
#include <math.h>
#include "rebound.h"

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_create_simulation();
// 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->heartbeat         = heartbeat;
r->usleep        = 5000;            // Slow down integration (for visualization only)

// Star
struct reb_particle star = {0};
star.m  = 1.;

// planet
struct reb_particle planet = {0};
planet.m  = 1e-3;
planet.x  = 1;
planet.vy = sqrt(r->G*(star.m+planet.m)/planet.x);

// 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_move_to_com(r);

system("rm -v a.txt");

reb_integrate(r, tmax);
}

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_output_check(r, M_PI*2.)){
reb_output_timing(r, tmax);
}
if(reb_output_check(r, M_PI*2.)){ // output every year
FILE* f = fopen("a.txt","a");
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