Difference between revisions of "RDYNLAB Code"
Jump to navigation
Jump to search
(New page: 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 ...) |
(No difference)
|
Revision as of 18:04, 18 June 2009
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***