Particle Trajectory Code
Jump to navigation
Jump to search
This is the old version of the code, which is now obsolete. The new version can be found here.
For some reason I am unable to load the source code files onto the wiki, so I have just copied and pasted them below. The boldface headers give the appropriate name for each file. You will need all of the files in order to compile the program.
Return to Simulating Particle Trajectories
*** build.sh ***
#!/bin/sh rm -f trajectory gcc -g -O2 -Wall -I. -o trajectory fields.c func.c init.c readinput.c trajectory.c -lm -lgsl -lgslcblas gcc -g -O2 -Wall -I. -o check_field fields.c check_field.c -lm
*** fields.c ***
#include <trajectory.h>
#include <math.h>
/* Functions to return the components of the E and B fields at position x,y,z */
double efield_x(double x, double y, double z)
{
return 0.0;
}
double efield_y(double x, double y, double z)
{
return 0.0;
}
double efield_z(double x, double y, double z)
{
return 0.0;
}
double bfield_x(double x, double y, double z)
{
return 0.0;
}
double bfield_y(double x, double y, double z)
{
return 0.0;
}
double bfield_z(double x, double y, double z)
{
/* Uniform B field: */
//return -0.0725;
/* Using the actual measured fields for the magnet. This functional form was
* determined using a seven parameter fit to the measured field */
double f,g;
double xcm,ycm;
double field;
double arg;
xcm=100.*(x-LMAGNET/2.);
ycm=100.*y;
arg=powf(fabs(xcm+0.064866),2.74157);
f=1.0019*exp(-fabs(arg)/(2.*23.4233*23.4233));
g=-0.0107702*ycm*ycm+0.174007*ycm+0.305011;
if(g<0.0) g=0.0;
field=-f*g*0.0896;
return field;
}
*** func.c ***
#include <trajectory.h>
#include <gsl/gsl_errno.h>
int func(double t, /* Unused, but solver requires this variable */
const double y[], /* The current position/momentum of the particle */
double f[], /* The time derivatives of the position/momentum */
void *params) /* Pointer to the parameters */
{
/* Note: array elements are x, y, z, px, py, pz!
* Start by reading the necessary parameters */
struct particle p=*(struct particle *)params;
double m=p.m;
double gamma=p.gamma;
double q=p.q;
/* Get the values for the components of the fields at this position */
double Ex=efield_x(y[0],y[1],y[2]);
double Ey=efield_y(y[0],y[1],y[2]);
double Ez=efield_z(y[0],y[1],y[2]);
double Bx=bfield_x(y[0],y[1],y[2]);
double By=bfield_y(y[0],y[1],y[2]);
double Bz=bfield_z(y[0],y[1],y[2]);
/* Calculate the derivatives */
f[0]=y[3]/(m*gamma);
f[1]=y[4]/(m*gamma);
f[2]=y[5]/(m*gamma);
f[3]=q*(Ex+(y[4]*Bz-y[5]*By)/(m*gamma));
f[4]=q*(Ey+(y[5]*Bx-y[3]*Bz)/(m*gamma));
f[5]=q*(Ez+(y[3]*By-y[4]*Bx)/(m*gamma));
return GSL_SUCCESS;
}
*** init.c ***
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <trajectory.h>
int init(struct particle *e, /* particle struct for electron */
struct particle *p, /* particle struct for positron */
struct parameters pars) /* parameters read from input file */
{
double ptot; /* Total momentum of a particle */
/* Set the initial position, energy, mass, charge, gamma, and components of
* momentum for the electron */
e->x=pars.dinset;
e->y=pars.hprime;
e->z=0.0;
e->e=pars.eelectron;
e->m=EMASS;
e->q=-ECHARGE;
e->gamma=e->e/(e->m*CLIGHT*CLIGHT);
if(e->gamma<1.00)
{
/* Total energy is less than rest energy-we need to toss this run*/
if(pars.type<1) fprintf(stdout,"Error: gamma for electron < 1.\n");
return 3;
}
ptot=sqrt(e->e*e->e/(CLIGHT*CLIGHT)-e->m*e->m*CLIGHT*CLIGHT);
e->px=ptot*cos(pars.theta);
e->py=ptot*sin(pars.theta);
e->pz=0.0;
/* Set the initial position, energy, mass, charge, gamma, and components of
* momentum for the positron */
p->x=e->x;
p->y=e->y;
p->z=e->z;
p->e=pars.egamma-pars.eelectron;
p->m=EMASS;
p->q=ECHARGE;
p->gamma=p->e/(p->m*CLIGHT*CLIGHT);
if(p->gamma<1.00)
{
/* Total energy is less than rest energy-we need to toss this run */
if(pars.type<1) fprintf(stdout,"Error: gamma for positron < 1.\n");
return 4;
}
ptot=sqrt(p->e*p->e/(CLIGHT*CLIGHT)-p->m*p->m*CLIGHT*CLIGHT);
p->px=ptot*cos(pars.theta);
p->py=ptot*sin(pars.theta);
p->pz=0.0;
/* Write out the initial status of each particle to the output file */
if(pars.type<1)
{
fprintf(stdout,"# Initial particle conditions:\n");
fprintf(stdout,"# Electron:\n");
fprintf(stdout,"# x: %g meters\n",e->x);
fprintf(stdout,"# y: %g meters\n",e->y);
fprintf(stdout,"# z: %g meters\n",e->z);
fprintf(stdout,"# m: %g kilograms\n",e->m);
fprintf(stdout,"# q: %g Coulombs\n",e->q);
fprintf(stdout,"# E: %g Joules\n",e->e);
fprintf(stdout,"# gamma: %g\n",e->gamma);
fprintf(stdout,"# p_x: %g kg m/s\n",e->px);
fprintf(stdout,"# p_y: %g kg m/s\n",e->py);
fprintf(stdout,"# p_z: %g kg m/s\n",e->pz);
fprintf(stdout,"# Positron:\n");
fprintf(stdout,"# x: %g meters\n",p->x);
fprintf(stdout,"# y: %g meters\n",p->y);
fprintf(stdout,"# z: %g meters\n",p->z);
fprintf(stdout,"# m: %g kilograms\n",p->m);
fprintf(stdout,"# q: %g Coulombs\n",p->q);
fprintf(stdout,"# E: %g Joules\n",p->e);
fprintf(stdout,"# gamma: %g\n",p->gamma);
fprintf(stdout,"# p_x: %g kg m/s\n",p->px);
fprintf(stdout,"# p_y: %g kg m/s\n",p->py);
fprintf(stdout,"# p_z: %g kg m/s\n",p->pz);
}
return 0;
}
*** readinput.c ***
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <trajectory.h>
void readinput(FILE *in, /* Pointer to input file */
struct parameters *pars) /* Pointer to parameters struct */
{
char buffer[BUFSIZ]; /* Read buffer */
fgets(buffer,BUFSIZ,in);
pars->type=atoi(buffer); /* Read in run type first (0=normal, >0=scan)*/
fgets(buffer,BUFSIZ,in);
pars->egamma=atof(buffer); /* Gamma ray energy, in MeV */
if(pars->type<1)
{
/* Read in for normal mode only. In scan mode, these parameters are
* chosen randomly */
fgets(buffer,BUFSIZ,in);
pars->eelectron=atof(buffer); /* Electron energy, in MeV */
fgets(buffer,BUFSIZ,in);
pars->theta=atof(buffer); /* Rotation of magnet relative to beam, in deg */
fgets(buffer,BUFSIZ,in);
pars->hprime=atof(buffer); /* Height of converter above yoke */
}
fgets(buffer,BUFSIZ,in);
pars->dinset=atof(buffer); /* Inset of converter relative to magnet */
/* This information always written to stdout (output file for normal runs,
* log file for scan runs */
fprintf(stdout,"# Parameters read from input file:\n");
fprintf(stdout,"# E_gamma (gamma ray energy): %g MeV\n",pars->egamma);
if(pars->type<1)
{
fprintf(stdout,"# E_electron (electron energy): %g MeV\n",
pars->eelectron);
fprintf(stdout,"# Theta (incident angle): %g deg.\n",pars->theta);
fprintf(stdout,"# H' (\"y\" location of pair production): %g m\n",
pars->hprime);
}
fprintf(stdout,"# D_inset (\"x\" location of pair production): %g m\n",
pars->dinset);
/* Convert to proper units */
pars->theta=pars->theta*M_PI/180.; /* Convert to radians */
pars->egamma=pars->egamma*1.602E-13; /* Convert to Joules */
pars->eelectron=pars->eelectron*1.602E-13; /* Convert to Joules */
}
*** trajectory.c ***
/* trajectory.c
* A program to simulate the paths of an electron-positron pair
* produced from the iteraction of a gamma ray with a converter
* in a non-uniform magnetic field. */
/* New in this version: instead of printing full output for each run in a scan,
* just print out a normal type run input file */
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>
#include <time.h>
#include <trajectory.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv.h>
int main(int argc, char *argv[])
{
struct parameters pars; /* Struct for holding operating parameters */
struct particle e,p; /* Particles structs for the electron and positron */
double t,t1,h,y[6]; /* Initial and final times, initial time step, and
array for passing six kinematic variables to
solver routine */
/* These variables are used to approximate the exit angles of the particles */
double ex1,ey1,ex2,ey2,px1,py1,px2,py2;
double slopee,slopep,phie,phip,anglee,anglep;
/* GSL solver types: Runge-Kutta 8th order with 9th order error estimate */
const gsl_odeiv_step_type *T1 = gsl_odeiv_step_rk8pd;
const gsl_odeiv_step_type *T2 = gsl_odeiv_step_rk8pd;
FILE *in,*out; /* I/O files */
int i,stat; /* miscellaneous counters and status variables */
double bx;
double slope,intercept;
int numiter; /* number of iterations */
char fname[BUFSIZ]; /* Character buffer for filenames */
char buffer[BUFSIZ]; /* Character buffer for system commands */
int flag;
/* This is a variable to flag the "goodness" of a set of parameters:
* 0: Electron hits return yoke.
* 1: Neither particle has an exit angle greater than GOODANGLE
* 2: Electron only has an exit angle greater than GOODANGLE
* 3: Positron only has an exit angle greater than GOODANGLE
* 4: Both particles have exit angles greater than GOODANGLE
* 5: Electron gets trapped in field (not good)
* 6: Positron gets trapped in field (not good) */
/* Check command line arguments and open input file */
if(argc<2)
{
fprintf(stderr,"You must specify an input file.\n");
exit(1);
}
fprintf(stdout,"# Using input file %s\n",argv[1]);
in=fopen(argv[1],"r");
if(!in)
{
fprintf(stderr,"Unable to open file %s\n",argv[1]);
exit(2);
}
/* Read input file */
readinput(in,&pars);
fclose(in);
/* Write some information to the output/log file - always to stdout */
fprintf(stdout,"# Magnet geometry:\n");
fprintf(stdout,"# Length: %g meters\n",LMAGNET);
fprintf(stdout,"# Height: %g meters\n",HMAGNET);
if(pars.type>0)
{
fprintf(stdout,"# Columns in table below are as follows:\n");
fprintf(stdout,"# (1) Iteration number\n");
fprintf(stdout,"# (2) Energy of electron (in MeV)\n");
fprintf(stdout,"# (3) Angle of magnet relative to beam (in deg.)\n");
fprintf(stdout,"# (4) Position of converter above return yoke (in m)\n");
fprintf(stdout,"# (5) \"Goodness\" flag:\n");
fprintf(stdout,"# 0: Electron doesn't clear return yoke.\n");
fprintf(stdout,"# 1: Neither particle has an exit angle > 90 deg.\n");
fprintf(stdout,"# 2: Electron has an exit angle > 90 deg.\n");
fprintf(stdout,"# 3: Positron has an exit angle > 90 deg.\n");
fprintf(stdout,"# 4: Both particles have exit angles > 90 deg.\n");
}
/* Begin scan loop */
if(pars.type<1) numiter=1;
else
{
numiter=pars.type;
srand(time(NULL));
}
for(i=1;i<=numiter;i++)
{
flag=-1;
/* Set random initial particle conditions and parameters, if this is a
* scan */
if(pars.type>0)
{
pars.eelectron=(double)rand()/(double)RAND_MAX*pars.egamma;
pars.theta=(double)rand()/(double)RAND_MAX*MAXANGLE*M_PI/180.;
pars.hprime=(double)rand()/(double)RAND_MAX*HMAGNET;
}
/* Open an output file, if this is a scan. Otherwise, we will direct
* output to stdout */
if(pars.type>0)
{
sprintf(fname,"%8.8d.pars",i);
out=fopen(fname,"w");
fprintf(out,"0 # Normal run type - generates verbose output\n");
fprintf(out,"%13.5E # Gamma ray energy (MeV)\n",
pars.egamma/1.602E-13);
fprintf(out,"%13.5E # Electron energy (MeV)\n",
pars.eelectron/1.602E-13);
fprintf(out,"%13.5E # Angle of photon relative to magnet (deg)\n",
pars.theta*180./M_PI);
fprintf(out,"%13.5E # \"Y\" location of converter above yoke (m)\n",
pars.hprime);
fprintf(out,"%13.5E # \"X\" location of converter from left side ",
pars.dinset);
fprintf(out,"of magnet (m)\n");
fclose(out);
}
/* Initialize the particles */
stat=init(&e,&p,pars);
if(stat==0)
{
/* This information is useful to have in the output */
if(pars.type<1)
{
fprintf(stdout,"# ELECTRON ENERGY (MEV): %g\n",
pars.eelectron/1.602E-13);
fprintf(stdout,"# RELATIVE ANGLE (DEG.): %g\n",pars.theta*180./M_PI);
fprintf(stdout,"# CONVERTER HEIGHT (M) : %g\n",pars.hprime);
}
/* Initialize the variables used to determine the exit angles */
ex1=e.x;
ey1=e.y;
ex2=ex1;
ey2=ey1;
px1=p.x;
py1=p.y;
px2=px1;
py2=py1;
/* Solve EOM. Continue process until both particles are out of the field.
* Do the electron first, then the positron! */
if(pars.type<1)
{
fprintf(stdout,"# Beginning solver iterations, electron first!\n");
fprintf(stdout,"# Columns in table below are as follows:\n");
fprintf(stdout,"# Time, in seconds.\n");
fprintf(stdout,"# X position of electron, in meters.\n");
fprintf(stdout,"# Y position of electron, in meters.\n");
fprintf(stdout,"# Z position of electron, in meters.\n");
}
/* Stepsize control: Keep accuracy to TINY. */
gsl_odeiv_step *S1 = gsl_odeiv_step_alloc (T1,6);
gsl_odeiv_control *C1 = gsl_odeiv_control_y_new (TINY,0.0);
gsl_odeiv_evolve *E1 = gsl_odeiv_evolve_alloc(6);
gsl_odeiv_system sys1 = {func, NULL, 6, &e};
t=0.0; /* Start at t=0 */
t1=LONGTIME; /* Iterate to time... We'll break long before this. */
h=FSTEP; /* Initial step size */
/* Initialize variable array for solver */
y[0]=e.x;
y[1]=e.y;
y[2]=e.z;
y[3]=e.px;
y[4]=e.py;
y[5]=e.pz;
/* Print out initial position of electron at time t=0 */
if(pars.type<1)
fprintf(stdout,"%13.5E %13.5E %13.5E %13.5E\n",t,y[0],y[1],y[2]);
/* Main solver loop */
while(t<t1)
{
int status=gsl_odeiv_evolve_apply(E1,C1,S1,&sys1,&t,t1,&h,y);
if(status!=GSL_SUCCESS) break;
/* Update particle struct and exit angle variables
* after each iteration */
e.x=y[0];
e.y=y[1];
e.z=y[2];
e.px=y[3];
e.py=y[4];
e.pz=y[5];
ex1=ex2;
ey1=ey2;
ex2=e.x;
ey2=e.y;
/* Print new position to output file */
if(pars.type<1)
fprintf(stdout,"%13.5E %13.5E %13.5E %13.5E\n",t,y[0],y[1],y[2]);
/* Break if electron hits the yoke */
if(ey2<=0.0 && ey1>0.0 && ex2>0.0 && ex2<LMAGNET)
{
flag=0;
break;
}
/* Break if particle is more than 40 cm away from the magnet */
if(y[0]>LMAGNET+0.40 || y[0]<-0.40 || y[1]<-0.40 || y[1]>HMAGNET+0.40)
break;
}
/* Check if the electron ever made it out of the magnet */
if(flag!=0)
{
if(ey2>-0.05 && ey2<HMAGNET+0.05 && ex2>-0.05 && ex2<LMAGNET+0.05)
flag=5;
}
/* Free the solver structures so they can be reused */
gsl_odeiv_evolve_free(E1);
gsl_odeiv_control_free(C1);
gsl_odeiv_step_free(S1);
/* Now repeat for positron, but start with two empty lines to output */
if(pars.type<1)
{
fprintf(stdout,"\n\n# Continuing iterations, for the positron\n");
fprintf(stdout,"# Columns in table below are as follows:\n");
fprintf(stdout,"# Time, in seconds.\n");
fprintf(stdout,"# X position of positron, in meters.\n");
fprintf(stdout,"# Y position of positron, in meters.\n");
fprintf(stdout,"# Z position of positron, in meters.\n");
}
/* Stepsize control: Keep accuracy to TINY. */
gsl_odeiv_step *S2 = gsl_odeiv_step_alloc (T2,6);
gsl_odeiv_control *C2 = gsl_odeiv_control_y_new (TINY,0.0);
gsl_odeiv_evolve *E2 = gsl_odeiv_evolve_alloc(6);
gsl_odeiv_system sys2 = {func, NULL, 6, &p};
t=0.0; /* Start at t=0 */
t1=LONGTIME; /* Iterate to time... We'll break long before this. */
h=FSTEP; /* Initial step size */
/* Initialize variable array for solver */
y[0]=p.x;
y[1]=p.y;
y[2]=p.z;
y[3]=p.px;
y[4]=p.py;
y[5]=p.pz;
/* Print out initial position of positron at time t=0 */
if(pars.type<1)
fprintf(stdout,"%13.5E %13.5E %13.5E %13.5E\n",t,y[0],y[1],y[2]);
/* Main solver loop */
while(t<t1)
{
int status=gsl_odeiv_evolve_apply(E2,C2,S2,&sys2,&t,t1,&h,y);
if(status!=GSL_SUCCESS) break;
/* Update particle struct and exit angle
* variables after each iteration */
p.x=y[0];
p.y=y[1];
p.z=y[2];
p.px=y[3];
p.py=y[4];
p.pz=y[5];
px1=px2;
py1=py2;
px2=p.x;
py2=p.y;
/* Print new position to output file */
if(pars.type<1)
fprintf(stdout,"%13.5E %13.5E %13.5E %13.5E\n",t,y[0],y[1],y[2]);
/* Break if particle is more than 40 cm away from the magnet */
if(y[0]>LMAGNET+0.40 || y[0]<-0.40 || y[1]<-0.40 || y[1]>HMAGNET+0.40)
break;
}
/* Check if the positron ever made it out of the magnet */
if(flag!=0)
{
if(py2>-0.05 && py2<HMAGNET+0.05 && px2>-0.05 && px2<LMAGNET+0.05)
{
if(flag!=5) flag=6;
else flag=7;
}
}
/* Free these variables so they can be reused later on */
gsl_odeiv_evolve_free(E2);
gsl_odeiv_control_free(C2);
gsl_odeiv_step_free(S2);
/* For reference/plotting purposes, write out points to define the
* beamline at 5 cm before and after the magnet. */
if(pars.type<1)
{
fprintf(stdout,"\n\n# These points define the beamline...\n");
for(bx=-0.05;bx<=LMAGNET+0.05;bx=bx+0.025)
{
slope=tan(pars.theta);
intercept=pars.hprime-pars.dinset*slope;
fprintf(stdout,"%13.5E %13.5E %13.5E\n",
bx, slope*(bx)+intercept,0.0);
}
}
/* Use ex1, ex2, etc... to determine exit angles relative to beam, and
* set the goodness flag */
if(flag<0)
{
flag=1;
slopee=(ey2-ey1)/(ex2-ex1);
phie=atan(slopee);
slopep=(py2-py1)/(px2-px1);
phip=atan(slopep);
if(phie<0) phie=phie+M_PI;
if(phip<0) phip=phip+M_PI;
anglee=M_PI-phie+pars.theta;
anglep=phip-pars.theta;
if(anglee>M_PI) anglee=anglee-M_PI;
if(anglep>M_PI) anglep=anglep-M_PI;
if(anglep>GOODANGLE*M_PI/180.) flag=3;
if(anglee>GOODANGLE*M_PI/180.) flag=flag+1;
}
/* Print results to stdout, if this is a scan */
if(pars.type>0)
{
fprintf(stdout,"%8d %13.5E %13.5E %13.5E %2d X\n",i,
pars.eelectron/1.602E-13,
pars.theta*180./M_PI,
pars.hprime,
flag);
/* Include this line to provide a counter on the term... */
fprintf(stderr,"Completed set %8d of %8d\n",i,numiter);
}
}
else if(pars.type>0) /* Don't count this iteration if scanning */
{
i--;
}
else /* In normal mode, this run was just no good. Exit */
{
return stat;
}
/* Delete the output file if this is a scan and flag <2 */
if(flag<2 && pars.type>0)
{
sprintf(buffer,"rm -f %s",fname);
system(buffer);
}
}
/* Should be done now. Exit gracefully. */
return 0;
}
*** trajectory.h ***
#ifndef _TRAJECTORY_H_
#define _TRAJECTORY_H_
#include <stdio.h>
#define CLIGHT 2.998E8 /* Speed of light, in m/s */
#define ECHARGE 1.602E-19 /* Electron charge, in C */
#define EMASS 9.11E-31 /* Electron mass, in kg */
#define LMAGNET 0.2032 /* Length of magnet */
#define HMAGNET 0.1524 /* Height of magnet */
#define MAXANGLE 75 /* Maximum angle of magnet relative to beam for scan */
#define GOODANGLE 90 /* Desired minimum exit angle for particles */
#define TINY 1.e-12 /* Required numerical accuracy */
#define FSTEP 1.e-18 /* First stepsize */
#define LONGTIME 3.e-8 /* Maximum iteration time */
struct particle
{
double x,y,z; /* specifies position of particle */
double px,py,pz; /* specifies (relativistic) momenta of the particle */
double e; /* total energy of particle */
double gamma; /* relativistic parameter for particle */
double m,q; /* mass and charge of the particle */
};
struct parameters
{
int type; /* Run type (0=normal, >0 = # of iterations for scan) */
double egamma; /* Energy of incident gamma ray */
double eelectron; /* Energy of the electron */
double theta; /* Angle of incident gamma relative to magnet axis */
double hprime; /* Height above yoke ("Y") where pair is produced */
double dinset; /* Inset distance ("X") where pair is produced */
};
/* Read the input file */
void readinput(FILE* ,struct parameters*);
/* Initialize particle structs */
int init(struct particle*, struct particle*, struct parameters);
/* Defines the six coupled ODE's */
int func(double, const double [], double [], void*);
/* These next six functions return the x, y, and z components of the E and B
* fields */
double efield_x(double, double, double);
double efield_y(double, double, double);
double efield_z(double, double, double);
double bfield_x(double, double, double);
double bfield_y(double, double, double);
double bfield_z(double, double, double);
#endif
Return to Simulating Particle Trajectories