Particle Trajectory Code

From New IAC Wiki
Revision as of 21:38, 10 June 2009 by Oborn (talk | contribs)
Jump to navigation Jump to search

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

  1. include <trajectory.h>
  2. 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 ***

  1. include <trajectory.h>
  2. 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 ***

  1. include <stdio.h>
  2. include <stdlib.h>
  3. include <math.h>
  4. 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 ***

  1. include <stdio.h>
  2. include <stdlib.h>
  3. include <math.h>
  4. 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 */
  1. include <stdio.h>
  2. include <stdlib.h>
  3. include <unistd.h>
  4. include <math.h>
  5. include <time.h>
  6. include <trajectory.h>
  7. include <gsl/gsl_errno.h>
  8. 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 ***

  1. ifndef _TRAJECTORY_H_
  2. define _TRAJECTORY_H_
  1. include <stdio.h>
  1. define CLIGHT 2.998E8 /* Speed of light, in m/s */
  2. define ECHARGE 1.602E-19 /* Electron charge, in C */
  3. define EMASS 9.11E-31 /* Electron mass, in kg */
  4. define LMAGNET 0.2032 /* Length of magnet */
  5. define HMAGNET 0.1524 /* Height of magnet */
  6. define MAXANGLE 75 /* Maximum angle of magnet relative to beam for scan */
  7. define GOODANGLE 90 /* Desired minimum exit angle for particles */
  8. define TINY 1.e-12 /* Required numerical accuracy */
  9. define FSTEP 1.e-18 /* First stepsize */
  10. 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);

  1. endif
Retrieved from "https://wiki.iac.isu.edu/index.php?title=Particle_Trajectory_Code&oldid=44234"