RDYNLAB Code
Jump to navigation
Jump to search
The following files constitute the source code for RDYNLAB. They can be compiled with any source code intended to use the library functions or (preferably) compiled as a shared object library using GNU autoconf, automake, and libtool. Compiling a shared object library is beyond the scope of this document. The package, with the appropriate configure scripts and makefiles, can be obtained by contacting the author (kelleyk@byui.edu).
It should be noted that the library is written in C++, and will probably only compile or work properly with other C++ codes.
***rdynlab.h***
/* RDYNLAB: Relativistic particle DYNamics in the LAB frame. */
#ifndef _RDYNLAB_H_
#define _RDYNLAB_H_
#include <iostream>
using namespace std;
/* Some global control parameters */
#define RDL_TINY 1.e+00 /* Required numerical accuracy */
#define RDL_TMAX 1.e+18 /* Maximum time we will ever be allowed to run to.
* This number can be rediculously large, since the
* stopping time will normally be provided by the
* user. */
/* A few constants that we will need */
#define RDL_CLIGHT 2.99792458E8 /* Speed of light in vacuum */
#define RDL_MU0 1.25663706144E-6 /* Permeability of free space, needed for power
* loss of charged particles */
/***************************** Vector Class **********************************/
class RDL_Vector
{
public:
double x,y,z; /* The three components of a vector */
/* Construction and destruction operators */
RDL_Vector() {};
~RDL_Vector() {};
RDL_Vector(double,double,double);
/* Operators */
RDL_Vector operator+(const RDL_Vector&); /* Add vectors */
RDL_Vector operator+(void); /* "Positive" of a vector */
RDL_Vector operator-(const RDL_Vector&); /* Subtract vectors */
RDL_Vector operator-(void); /* "Negative" of a vector */
RDL_Vector operator*(double); /* Multiply vector by scalar */
double operator*(const RDL_Vector&); /* Vector dot product */
RDL_Vector operator/(double); /* Divide vector by scalar */
RDL_Vector operator^(const RDL_Vector&); /* Vector cross product */
bool operator==(const RDL_Vector&); /* Are vectors the same? */
/* Friends to this class */
friend ostream &operator<<(ostream &out,const RDL_Vector &V);
friend RDL_Vector operator*(double a, const RDL_Vector &V);
};
/* Define << operator for vectors */
ostream &operator<<(ostream &out, const RDL_Vector &V);
/* Allows us to do scalar times vector */
RDL_Vector operator*(double a, const RDL_Vector &V);
/* These return the length of a vector and a unit vector */
double length(const RDL_Vector &V);
RDL_Vector unit(const RDL_Vector &V);
/***************************** Particle Class ********************************/
class RDL_Particle
{
public:
double q; /* Charge, in Coulombs */
double m; /* Mass, in kg */
RDL_Vector x; /* Position, in meters */
RDL_Vector v; /* Velocity, in m/s */
/* Construction and Destruction operators */
RDL_Particle() {};
~RDL_Particle() {};
RDL_Particle(double, double, const RDL_Vector&, const RDL_Vector&);
/* Init with q,m,x,v */
/* We need to define the = operator */
RDL_Particle operator=(const RDL_Particle&);
/* Allow the << operator access... */
friend ostream &operator<<(ostream &out, const RDL_Particle &p);
/* Return particle properties */
double gamma(void); /* Returns relativistic factor gamma */
double energy(void); /* Returns total energy */
double kinetic_energy(void); /* Returns kinetic energy */
double rest_energy(void); /* Returns rest energy */
RDL_Vector momentum(void); /* Returns momentum vector */
};
/* Overload << operator for particles */
ostream &operator<<(ostream &out, const RDL_Particle &p);
/******************* Struct for passing information to GSL ********************/
struct RDL_solver_parameters
{
RDL_Particle p; /* The particle */
double t; /* The time */
RDL_Vector (*f)(const RDL_Particle&,double); /* Function returning force */
};
/******************** Declaration of main solver routine **********************/
/* This is the main solver routine. Input parameters are the particle,
* a function that returns the force, the initial time, the initial time step,
* a function that determines if the simulation should stop, and an ostream for
* writing output.
*
* The force function takes arguments of a particle (which will include position
* and velocity) and time. The stopping condition function takes arguments of a
* particle (which has the particle's position and velocity), the time, and the
* number of iterations the solver has completed. */
int RDL_solver(RDL_Particle&,
RDL_Vector (*)(const RDL_Particle&,double),
double, double,
int (*)(const RDL_Particle&,double,int),
ostream&);
/* This function is required by the solver routine, and will probably never
* need to be called from user space. */
int RDL_solver_func(double t,
const double y[],
double f[],
void *params);
#endif
***RDLvector.cpp***
/* RDLvector.c: 3-D vectors and functions */
#include <rdynlab.h>
#include <math.h>
using namespace std;
/* Constructor */
RDL_Vector::RDL_Vector (double xin, double yin, double zin)
{
x = xin;
y = yin;
z = zin;
}
/* Operators */
RDL_Vector RDL_Vector::operator+(const RDL_Vector ¶m)
{
RDL_Vector temp;
temp.x = x + param.x;
temp.y = y + param.y;
temp.z = z + param.z;
return temp;
}
RDL_Vector RDL_Vector::operator+(void)
{
RDL_Vector temp;
temp.x = x;
temp.y = y;
temp.z = z;
return temp;
}
RDL_Vector RDL_Vector::operator-(const RDL_Vector ¶m)
{
RDL_Vector temp;
temp.x = x - param.x;
temp.y = y - param.y;
temp.z = z - param.z;
return temp;
}
RDL_Vector RDL_Vector::operator-(void)
{
RDL_Vector temp;
temp.x = -x;
temp.y = -y;
temp.z = -z;
return temp;
}
RDL_Vector RDL_Vector::operator*(double a)
{
RDL_Vector temp;
temp.x = x*a;
temp.y = y*a;
temp.z = z*a;
return temp;
}
double RDL_Vector::operator*(const RDL_Vector ¶m)
{
double temp;
temp = (x * param.x) + (y * param.y) + (z * param.z);
return temp;
}
RDL_Vector RDL_Vector::operator/(double a)
{
RDL_Vector temp;
temp.x = x/a;
temp.y = y/a;
temp.z = z/a;
return temp;
}
RDL_Vector RDL_Vector::operator^(const RDL_Vector ¶m)
{
RDL_Vector temp;
temp.x = (y * param.z) - (z * param.y);
temp.y = (z * param.x) - (x * param.z);
temp.z = (x * param.y) - (y * param.x);
return temp;
}
bool RDL_Vector::operator==(const RDL_Vector ¶m)
{
if(param.x==x && param.y==y && param.z==z) return true;
else return false;
}
/* << operator */
ostream &operator<<(ostream &out, const RDL_Vector &V)
{
out << "<" << V.x << ", " << V.y << ", " << V.z << ">";
return out;
}
/* This one allows double*RDL_Vector */
RDL_Vector operator*(double a, const RDL_Vector &V)
{
RDL_Vector temp;
temp.x = V.x*a;
temp.y = V.y*a;
temp.z = V.z*a;
return (temp);
}
/* Return the length of a vector, unit vector */
double length(const RDL_Vector &V)
{
RDL_Vector temp=V;
return sqrt(temp*temp);
}
RDL_Vector unit(const RDL_Vector &V)
{
RDL_Vector temp=V;
double s=length(temp);
/* Return zero if the length is zero */
if(s>0) return temp/s;
else return temp;
}
***RDLparticle.cpp***
/* RDLparticle.cpp: defines class for a particle */
#include <rdynlab.h>
#include <math.h>
RDL_Particle::RDL_Particle(double qin, double min, const RDL_Vector &xin, const RDL_Vector &vin)
{
q=qin;
m=min;
x=xin;
v=vin;
}
RDL_Particle RDL_Particle::operator=(const RDL_Particle ¶m)
{
q=param.q;
m=param.m;
x=param.x;
v=param.v;
return *this;
}
double RDL_Particle::gamma(void)
{
double s=length(v);
return 1./sqrt(1.-s*s/(RDL_CLIGHT*RDL_CLIGHT));
}
double RDL_Particle::energy(void)
{
return gamma()*m*RDL_CLIGHT*RDL_CLIGHT;
}
double RDL_Particle::kinetic_energy(void)
{
return (gamma()-1.0)*m*RDL_CLIGHT*RDL_CLIGHT;
}
double RDL_Particle::rest_energy(void)
{
return m*RDL_CLIGHT*RDL_CLIGHT;
}
RDL_Vector RDL_Particle::momentum(void)
{
return v*gamma()*m;
}
ostream &operator<<(ostream &out, const RDL_Particle &p)
{
out << "Current state of particle:\n";
out << " mass: " << p.m << "\n";
out << " charge: " << p.q << "\n";
out << " position: " << p.x << "\n";
out << " velocity: " << p.v << "\n";
return out;
}
***RDLsolver.cpp***
/* RDLsolver.cpp: The main solver routine for the RDYNLAB package */
#include <rdynlab.h>
#include <gsl/gsl_odeiv.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_diff.h>
#include <math.h>
#include <iostream>
#include <iomanip>
int RDL_solver(RDL_Particle &p,
RDL_Vector (*force) (const RDL_Particle&, double),
double ti,
double dt,
int (*stopping_conditions) (const RDL_Particle&, double, int),
ostream &out)
{
double y[6];
double t;
int status,iter,retval;
struct RDL_solver_parameters pars;
/* Initialize the parameters struct */
pars.p=p;
pars.t=ti;
pars.f=force;
/* Write information to ostream... */
out << scientific;
out << setprecision(3);
out << setfill(' ');
out << "# Beginning solver iterations. Output columns are as follows\n";
out << "# (1) Time, in seconds\n";
out << "# (2) X position of particle, in meters\n";
out << "# (3) Y position of particle, in meters\n";
out << "# (4) Z position of particle, in meters\n";
out << "# (5) X component of velocity, in m/s\n";
out << "# (6) Y component of velocity, in m/s\n";
out << "# (7) Z component of velocity, in m/s\n";
/* Set up GSL ODE solver system */
const gsl_odeiv_step_type *T=gsl_odeiv_step_rk8pd;
gsl_odeiv_step *s=gsl_odeiv_step_alloc(T,6);
gsl_odeiv_control *c=gsl_odeiv_control_y_new(RDL_TINY,0.0);
gsl_odeiv_evolve *e=gsl_odeiv_evolve_alloc(6);
gsl_odeiv_system sys={RDL_solver_func,NULL,6,&pars};
/* Initialize variable array */
y[0]=p.x.x;
y[1]=p.x.y;
y[2]=p.x.z;
y[3]=p.v.x;
y[4]=p.v.y;
y[5]=p.v.z;
/* Print out initial position */
out << setw(10) << ti << " ";
out << setw(10) << y[0] << " ";
out << setw(10) << y[1] << " ";
out << setw(10) << y[2] << " ";
out << setw(10) << y[3] << " ";
out << setw(10) << y[4] << " ";
out << setw(10) << y[5] << "\n";
/* Solver loop */
t=ti;
status=GSL_SUCCESS;
iter=1;
while (t<(RDL_TMAX+ti))
{
pars.p=p;
pars.t=t;
status=gsl_odeiv_evolve_apply(e,c,s,&sys,&t,(RDL_TMAX+ti),&dt,y);
out << setw(10) << t << " ";
out << setw(10) << y[0] << " ";
out << setw(10) << y[1] << " ";
out << setw(10) << y[2] << " ";
out << setw(10) << y[3] << " ";
out << setw(10) << y[4] << " ";
out << setw(10) << y[5] << "\n";
/* Update the particle and the struct that carries the info */
p.x.x=y[0];
p.x.y=y[1];
p.x.z=y[2];
p.v.x=y[3];
p.v.y=y[4];
p.v.z=y[5];
/* Check break conditions */
if(status!=GSL_SUCCESS) break;
retval=stopping_conditions(p,t,iter);
if(retval) break;
iter++;
}
/* Free memory allocated by GSL */
gsl_odeiv_evolve_free(e);
gsl_odeiv_control_free(c);
gsl_odeiv_step_free(s);
/* Return */
if(status!=GSL_SUCCESS) return 0;
else return retval;
}
***RDLsolver_func.cpp***
/* RDLsolver_func.cpp: This function is needed by the RDYNLAB main routine, and
* should net be called from user space */
#include <rdynlab.h>
#include <gsl/gsl_odeiv.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_diff.h>
#include <math.h>
#include <iostream>
#include <iomanip>
int RDL_solver_func(double t,const double y[],double f[], void *params)
{
/* Note: variables coming in array y are as follows:
* y[0]: x
* y[1]: y
* y[2]: z
* y[3]: vx
* y[4]: vy
* y[5]: vz */
/* Extract the particle and force from the input parameters */
RDL_solver_parameters *pars=(RDL_solver_parameters*)params;
RDL_Particle p=(pars->p);
RDL_Vector (*F)(const RDL_Particle&,const double)=pars->f;
/* Other variables we will need */
RDL_Vector Fpar,Fperp;
RDL_Vector a;
double gm,g3m;
double dgdt,g;
/* Update position and velocity of the particle from the y array. If I
* don't do this, the solver has no way of testing the accuracy of a step! */
p.x.x=y[0];
p.x.y=y[1];
p.x.z=y[2];
p.v.x=y[3];
p.v.y=y[4];
p.v.z=y[5];
/* Calculate gamma*m and gamma^3*m */
g=p.gamma();
gm=g*p.m;
g3m=g*g*gm;
/* Calculate the components of the force that are parallel and
* perpendicular to the velocity of the particle */
if(length(p.v)>0.)
{
Fpar=unit(p.v)*(F(p,t)*p.v)/length(p.v);
Fperp=F(p,t)-Fpar;
}
else /* Velocity is zero, no component of force is parallel to velocity! */
{
Fpar=p.v;
Fperp=F(p,t);
}
/* Calculate the acceleration w/o radiative losses (dv/dt) */
a=Fperp/gm+Fpar/g3m;
/* Use this acceleration to calculate dgamma/dt, derived from the
* Lienard formula: dgdt = 1/(mc^2) dE/dt */
dgdt=-RDL_MU0*p.q*p.q*g*g*g*g*g*g;
dgdt=dgdt/(6.*M_PI*p.m*RDL_CLIGHT*RDL_CLIGHT*RDL_CLIGHT);
dgdt=dgdt*(a*a-(p.v^a)*(p.v^a)/(RDL_CLIGHT*RDL_CLIGHT));
/* Update acceleration with radiative losses */
a=a+p.v*dgdt/g;
/* Put the time derivatives in the appropriate array */
f[0]=y[3];
f[1]=y[4];
f[2]=y[5];
f[3]=a.x;
f[4]=a.y;
f[5]=a.z;
/* return */
return GSL_SUCCESS;
}