Difference between revisions of "Particle Trajectory Code"
Jump to navigation
Jump to search
| Line 4: | Line 4: | ||
'''*** build.sh ***''' | '''*** build.sh ***''' | ||
| − | <pre><nowiki | + | <pre><nowiki> |
#!/bin/sh | #!/bin/sh | ||
| Line 12: | Line 12: | ||
gcc -g -O2 -Wall -I. -o check_field fields.c check_field.c -lm | gcc -g -O2 -Wall -I. -o check_field fields.c check_field.c -lm | ||
| − | + | </nowiki></pre> | |
---- | ---- | ||
Revision as of 21:42, 10 June 2009
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.
*** 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(ey1<=0.0 && ey2>=0.0 && ex1>ex2)
{
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(ey1>-0.05 && ey1<HMAGNET+0.05 && ex1>-0.05 && ex1<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(ey1>-0.05 && ey1<HMAGNET+0.05 && ex1>-0.05 && ex1<LMAGNET+0.05)
flag=6;
/* 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