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***
***RDLsolver.cpp***
***RDLsolver_func.cpp***