Oborn at 18:06, 18 June 2009

2009-06-18T18:06:22Z

Oborn: 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 ...

2009-06-18T18:04:32Z

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

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

<pre><nowiki>

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

</nowiki></pre>

----

'''***RDLvector.cpp***'''

<pre><nowiki>

/* 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;
}

</nowiki></pre>

----

'''***RDLparticle.cpp***'''

<pre><nowiki>

</nowiki></pre>

----

'''***RDLsolver.cpp***'''

<pre><nowiki>

</nowiki></pre>

----

'''***RDLsolver_func.cpp***'''

<pre><nowiki>

</nowiki></pre>

RDYNLAB Code - Revision history

Oborn at 18:06, 18 June 2009

Oborn: 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 ...