Difference between revisions of "Particle Trajectory Code"

From New IAC Wiki
Jump to navigation Jump to search
Line 4: Line 4:
 
'''*** build.sh ***'''
 
'''*** build.sh ***'''
  
−
<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>
+
</nowiki></pre>
  
 
----
 
----
 
'''*** fields.c ***'''
 
'''*** fields.c ***'''
  
 +
<pre><nowiki>
 
#include <trajectory.h>
 
#include <trajectory.h>
 
#include <math.h>
 
#include <math.h>
Line 66: Line 67:
 
   return field;
 
   return field;
 
}
 
}
 +
</nowiki></pre>
  
 
----
 
----
 
'''*** func.c ***'''
 
'''*** func.c ***'''
  
 +
<pre><nowiki>
 
#include <trajectory.h>
 
#include <trajectory.h>
 
#include <gsl/gsl_errno.h>
 
#include <gsl/gsl_errno.h>
Line 103: Line 106:
 
   return GSL_SUCCESS;
 
   return GSL_SUCCESS;
 
}
 
}
 +
</nowiki></pre>
  
 
----
 
----
 
'''*** init.c ***'''
 
'''*** init.c ***'''
  
 +
<pre><nowiki>
 
#include <stdio.h>
 
#include <stdio.h>
 
#include <stdlib.h>
 
#include <stdlib.h>
Line 188: Line 193:
 
   return 0;
 
   return 0;
 
}
 
}
 +
</nowiki></pre>
  
 
----
 
----
 
'''*** readinput.c ***'''
 
'''*** readinput.c ***'''
  
 +
<pre><nowiki>
 
#include <stdio.h>
 
#include <stdio.h>
 
#include <stdlib.h>
 
#include <stdlib.h>
Line 242: Line 249:
 
   pars->eelectron=pars->eelectron*1.602E-13; /* Convert to Joules */
 
   pars->eelectron=pars->eelectron*1.602E-13; /* Convert to Joules */
 
}
 
}
 +
</nowiki></pre>
  
 
----
 
----
 
'''*** trajectory.c ***'''
 
'''*** trajectory.c ***'''
  
 +
<pre><nowiki>
 
/* trajectory.c  
 
/* trajectory.c  
 
  * A program to simulate the paths of an electron-positron pair
 
  * A program to simulate the paths of an electron-positron pair
Line 633: Line 642:
 
   return 0;
 
   return 0;
 
}
 
}
 +
</nowiki></pre>
  
 
----
 
----
 
'''*** trajectory.h ***'''
 
'''*** trajectory.h ***'''
  
 +
<pre><nowiki>
 
#ifndef _TRAJECTORY_H_
 
#ifndef _TRAJECTORY_H_
 
#define _TRAJECTORY_H_
 
#define _TRAJECTORY_H_
Line 688: Line 699:
  
 
#endif
 
#endif
 +
</nowiki></pre>

Revision as of 21:40, 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
Retrieved from "https://wiki.iac.isu.edu/index.php?title=Particle_Trajectory_Code&oldid=44235"