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***

---

***RDLsolver.cpp***

---

***RDLsolver_func.cpp***