RDYNLAB Code

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 &param) 
{
  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 &param)
{
  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 &param)
{
  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 &param)
{
  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 &param)
{
  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 &param)
{
  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;
}