Difference between revisions of "New Particle Trajectory Code"

From New IAC Wiki
Jump to navigation Jump to search
 
Line 799: Line 799:
 
----
 
----
  
−
'''*** NAME ***'''
+
'''*** readinput.cpp ***'''
  
 
<pre><nowiki>
 
<pre><nowiki>
  
 +
#include <pairspec.h>
 +
#include <stdio.h>
 +
#include <stdlib.h>
 +
#include <math.h>
 +
 +
/* Read control parameters from the input file */
 +
 +
void readinput(const char *fname, /* Filename (input) */
 +
              int *type,        /* Run type (return) */
 +
      double *eg,        /* Gamma ray energy (return) */
 +
      double *ee,        /* Electron energy (return) */
 +
      double *theta,    /* Angle of magnet axis (return) */
 +
      double *xi,        /* X position of production (return) */
 +
      double *yi)        /* Y position of production (return) */
 +
{
 +
  FILE *f;
 +
  char buffer[BUFSIZ];
 +
 +
  /* Open the file */
 +
 +
  f=fopen(fname,"r");
 +
 +
  /* Check that we were able to open file for reading */
 +
 +
  if(!f)
 +
  {
 +
    fprintf(stderr,"Cannot read file %s\n",fname);
 +
    exit(2);
 +
  }
 +
 +
  /* Read in run type */
 +
 +
  fgets(buffer,BUFSIZ,f);
 +
  *type=atoi(buffer);
 +
 +
  /* Read in the gamma ray energy, in MeV */
 +
 +
  fgets(buffer,BUFSIZ,f);
 +
  *eg=atof(buffer);
 +
 +
  /* If this is not a scan run, read in the electron energy (in MeV), the
 +
  * angle of the beam relative to the magnet axis (in degrees), and the
 +
  * Y position of the production site (in meters) */
 +
 
 +
  if(*type<1)
 +
  {
 +
    fgets(buffer,BUFSIZ,f);
 +
    *ee=atof(buffer);
 +
    fgets(buffer,BUFSIZ,f);
 +
    *theta=atof(buffer);
 +
    fgets(buffer,BUFSIZ,f);
 +
    *yi=atof(buffer);
 +
  }
 +
 +
  /* Read in the X position of the production site (in meters) */
 +
 +
  fgets(buffer,BUFSIZ,f);
 +
  *xi=atof(buffer);
 +
 +
  /* Write the data read in to stdout */
 +
 +
  fprintf(stdout,"# Parameters read from input file %s:\n",fname);
 +
 +
  if(*type>0) fprintf(stdout,"# Run type: scan with %d iterations\n",*type);
 +
  else if(*type==0) fprintf(stdout,"# Run type: single\n");
 +
  else fprintf(stdout,"# Run type: interactive\n");
 +
 +
  fprintf(stdout,"#  Gamma ray energy (MeV):        %g\n",*eg);
 +
  if(*type<1)
 +
  {
 +
    fprintf(stdout,"#  Electron energy (MeV):        %g\n",*ee);
 +
    fprintf(stdout,"#  Beam angle (degrees):          %g\n",*theta);
 +
    fprintf(stdout,"#  Y of production site (meters): %g\n",*yi);
 +
  }
 +
  fprintf(stdout,"#  X of production site (meters): %g\n",*xi);
 +
 +
  /* Convert all input quantities into SI units (RDYNLAB uses SI units) */
 +
 +
  *theta=*theta*M_PI/180.;
 +
  *eg=*eg*1.602E-13;
 +
  *ee=*ee*1.602E-13;
 +
 +
  /* Close the input file and return */
 +
 +
  fclose(f);
 +
  return;
 +
}
  
 
</nowiki></pre>
 
</nowiki></pre>
Line 809: Line 896:
 
----
 
----
  
−
'''*** NAME ***'''
+
'''*** scanheader.cpp ***'''
  
 
<pre><nowiki>
 
<pre><nowiki>
  
 +
#include <stdio.h>
 +
#include <pairspec.h>
 +
 +
/* Write out the header for the scan mode output file. This was just too bulky
 +
* in the main routine */
 +
 +
void scanheader(void)
 +
{
 +
  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: Both particles have exit angles > %g deg.\n",
 +
    (double)GOODANGLE);
 +
  fprintf(stdout,"#      1: Electron has an exit angle > %g deg.\n",
 +
    (double)GOODANGLE);
 +
  fprintf(stdout,"#      2: Positron has an exit angle > %g deg.\n",
 +
    (double)GOODANGLE);
 +
  fprintf(stdout,"#      3: Neither particle has an exit angle > %g deg.\n",
 +
    (double)GOODANGLE);
 +
  fprintf(stdout,"#      4: Electron hits yoke.\n");
 +
  fprintf(stdout,"#      5: Positron hits yoke.\n");
 +
  fprintf(stdout,"#      6: Both particles hit yoke.\n");
 +
  fprintf(stdout,"#      7: Electron is trapped in field.\n");
 +
  fprintf(stdout,"#      8: Positron is trapped in field.\n");
 +
  fprintf(stdout,"#      9: Both particles are trapped in field.\n");
 +
  return;
 +
}
  
 
</nowiki></pre>
 
</nowiki></pre>
Line 819: Line 936:
 
----
 
----
  
−
'''*** NAME ***'''
+
'''*** setflag.cpp ***'''
  
 
<pre><nowiki>
 
<pre><nowiki>
  
 +
#include <rdynlab.h>
 +
#include <pairspec.h>
 +
#include <math.h>
 +
 +
/* Set the output flag */
 +
 +
int setflag(int state,    /* State of electron, as returned by stopcond */
 +
            int statp,    /* State of positron, as returned by stopcond */
 +
    double theta,  /* Beam angle */
 +
    RDL_Vector ev, /* Final velocity of electron */
 +
    RDL_Vector pv) /* Final velocity of positron */
 +
{
 +
  int flag;
 +
  double phie,phip;
 +
  double garad;
 +
  RDL_Vector beam;
 +
 +
  flag=0;
 +
 +
  /* Notes: if return values state/statep are 1, particle is trapped
 +
  * in field. If return values are 2, particle hits yoke. If return
 +
  * values are 3, then both particles leave the field successfully. */
 +
 
 +
  while(2>1) /* Allows me to break once flag is set */
 +
  {
 +
    /* Check for yoke impacts */
 +
    if(state==2 && statp!=2) { flag=4; break; }
 +
    if(state!=2 && statp==2) { flag=5; break; }
 +
    if(state==2 && statp==2) { flag=6; break; }
 +
    /* Check for trapped particles */
 +
    if(state==1 && statp!=1) { flag=7; break; }
 +
    if(state!=1 && statp==1) { flag=8; break; }
 +
    if(state==1 && statp==1) { flag=9; break; }
 +
 +
    /* Calculate exit angles from final velocities */
 +
 +
    beam.x=cos(theta);
 +
    beam.y=sin(theta);
 +
    phie=acos((ev*beam)/(length(ev)*length(beam)));
 +
    phip=acos((pv*beam)/(length(pv)*length(beam)));
 +
 +
    garad=GOODANGLE*M_PI/180.;
 +
    if(phie>=garad && phip>=garad) { flag=0; break; }
 +
    if(phie>=garad && phip<garad)  { flag=1; break; }
 +
    if(phie<garad && phip>=garad)  { flag=2; break; }
 +
   
 +
    flag=3;
 +
    break;
 +
  }
 +
 +
  return flag;
 +
}
  
 
</nowiki></pre>
 
</nowiki></pre>
Line 829: Line 998:
 
----
 
----
  
−
'''*** NAME ***'''
+
'''*** stopcond.cpp ***'''
  
 
<pre><nowiki>
 
<pre><nowiki>
  
 +
#include <rdynlab.h>
 +
#include <pairspec.h>
 +
 +
/* Stopping conditions for calculating particle trajectories */
 +
 +
int stopcond(const RDL_Particle &p, double t,int iter)
 +
{
 +
  double m,b,x0;
 +
 +
  /* If the particle is still hanging around after this long, then it's
 +
  * trapped in the field */
 +
 +
  if(t>LONGTIME) return 1;
 +
 +
  /* Check to see if the particle has hit the yoke */
 +
 +
  if(p.v.y<0 && p.x.y<0) /* Only happens if y velocity is negative
 +
                          * and particle is below y=0 */
 +
  {
 +
    m=p.v.y/p.v.x;
 +
    b=p.x.y-m*p.x.x;
 +
    x0=-b/m;
 +
    if(x0>=0.0 && x0<=LMAGNET) return 2;
 +
  }
 +
 +
  /* Check to see if the particle has left some arbitrary bounds */
 +
 +
  if(p.x.x>2.0*LMAGNET || p.x.x<-LMAGNET || p.x.y>2.0*HMAGNET || p.x.y<-HMAGNET)
 +
    return 3;
 +
 +
  /* Stopping conditions have not been met. Return zero */
 +
 +
  return 0;
 +
}
  
 
</nowiki></pre>
 
</nowiki></pre>

Latest revision as of 18:16, 18 June 2009

I am unable to post code to the wiki, so I have copied and pasted it all here. You will need all of these files in order to compile the code. There are dependencies on GSL, PGPLOT, ncurses, and RDYNLAB.

Return to Simulating Particle Trajectories

*** build.sh ***

This is the script I use to build the program on inca. You will need to make the appropriate changes to locate the include files and libraries. 


#!/bin/sh

rm -f pairspec

g++ -g -O2 -Wall \
  -I. -I/home/kellkevi/local/include \
  -o pairspec \
  main.cpp \
  readinput.cpp \
  randomize.cpp \
  initialize.cpp \
  lforce.cpp \
  stopcond.cpp \
  scanheader.cpp \
  setflag.cpp \
  initplot.cpp \
  interactive.cpp \
  -L/usr/X11R6/lib -lX11 \
  -L/home/kellkevi/local/pgplot -lcpgplot -lpgplot \
  -L/home/kellkevi/local/lib -lrdynlab \
  -lm -lgsl -lgslcblas -lg2c -lncurses


*** initialize.cpp *** 


#include <pairspec.h>
#include <rdynlab.h>
#include <math.h>

/* Initialize the electron and positron particle structs based on the 
 * control parameters */

void initialize(double eg,        /* Gamma ray energy (input) */
                double ee,        /* Electron energy (input) */
		double theta,     /* Beam angle (input) */
		double yi,        /* Y of conversion (input) */
		double xi,        /* X of conversion (input) */
		RDL_Particle *e,  /* Electron (modify) */
		RDL_Particle *p)  /* RDL_Positron (modify) */
{
  double gamma,erest,ep,speed;
  RDL_Vector pos=RDL_Vector(xi,yi,0.0);

  /* Calculate positron energy and particle rest energy*/

  ep=eg-ee;
  erest=EMASS*RDL_CLIGHT*RDL_CLIGHT;
 
  /* Calculate initial velocity vector for electron */

  gamma=ee/erest;
  speed=RDL_CLIGHT*sqrt(1.-1./(gamma*gamma));
  RDL_Vector evi=RDL_Vector(speed*cos(theta),speed*sin(theta),0.0);

  /* Calculate initial velocity vector for positron */

  gamma=ep/erest;
  speed=RDL_CLIGHT*sqrt(1.-1./(gamma*gamma));
  RDL_Vector pvi=RDL_Vector(speed*cos(theta),speed*sin(theta),0.0);

  /* Set all necessary values for both particles */

  e->q=-ECHARGE;
  e->m=EMASS;
  e->x=pos;
  e->v=evi;

  p->q=ECHARGE;
  p->m=EMASS;
  p->x=pos;
  p->v=pvi;

  /* Return */

  return;
}


*** initplot.cpp *** 


#include <cpgplot.h>
#include <pairspec.h>
#include <ncurses.h>
#include <math.h>

/* Initialize PGPLOT and ncurses for interactive mode */

int initplot(int *r, 
             int *c, 
             double *lstepsize, 
	     double *estepsize, 
	     double *astepsize)
{
  int pgpid;
  int rows,cols;
  char title[25];
  int i;

  /* Initialize ncurses environment */

  initscr();
  noecho();
  cbreak();
  keypad(stdscr,TRUE);
  getmaxyx(stdscr,rows,cols);
  box(stdscr,'|','-');
  attron(A_BOLD | A_UNDERLINE);
  sprintf(title,"PAIRSPEC Interactive Mode");
  mvprintw(2,(cols-25)/2,title);
  attroff(A_BOLD | A_UNDERLINE);
  for(i=3;i<cols-3;i++) mvprintw(4,i,"-");
  attron(A_BOLD);
  mvprintw(5,(cols-9)/2,"Commands:");
  attroff(A_BOLD);
  mvprintw(6,(cols-73)/2,
   "(PgUp/PgDn) Incr/Decr gamma energy    (g) Manually enter gamma ray energy");
  mvprintw(7,(cols-73)/2,
   "(Home/End)  Incr/Decr electron energy (e) Manually enter electron energy ");
  mvprintw(8,(cols-73)/2,
   "(Ins/Del)   Incr/Decr beam angle      (t) Manually enter beam angle      ");
  mvprintw(9,(cols-73)/2,
   "(Lft/Rght)  Incr/Decr conversion X    (x) Manually enter conversion X    ");
  mvprintw(10,(cols-73)/2,
   "(Up/Down)   Incr/Decr conversion Y    (y) Manually enter conversion Y    ");
  mvprintw(12,(cols-73)/2,
   "(=) Larger stepsizes      (-) Smaller stepsizes    (0) Original stepsizes");
  mvprintw(13,(cols-73)/2,
   "                (q) Write a single mode input file and exit              ");
  for(i=3;i<cols-3;i++) mvprintw(14,i,"-");
  refresh();

  *lstepsize=LSTEPSIZE;
  *estepsize=ESTEPSIZE*1.602E-13;
  *astepsize=ASTEPSIZE*M_PI/180;

  /* Initialize PGPLOT device */

  pgpid=cpgopen("/XWINDOW");
  cpgslct(pgpid);
  cpgpap(6,HMAGNET/LMAGNET);
  cpgswin(-LMAGNET,2.*LMAGNET,-HMAGNET,2.*HMAGNET);
  cpgsvp(0.,1.,0.,1.);
  cpgask(0);
  cpgscr(0,1.,1.,1.);
  cpgscr(1,0.,0.,0.);
  cpgpage();

  *r=rows;
  *c=cols;

  return pgpid;
}


*** interactive.cpp *** 


#include <cpgplot.h>
#include <stdio.h>
#include <pairspec.h>
#include <math.h>
#include <ncurses.h>

/* User interaction with the program */

void go_interactive(double *eg, 
                    double *ee, 
		    double *theta,
		    double *yi,
		    double *xi,
		    int *i,
		    int rows,
		    int cols,
		    double *lstepsize,
		    double *estepsize,
		    double *astepsize)
{
  char buffer[BUFSIZ];
  double emin;
  FILE *out;
  int icte,ictp,j;
  double m,b,bx;
  double *ex=NULL,*ey=NULL,*px=NULL,*py=NULL;
  int ctrl;

  emin=EMASS*RDL_CLIGHT*RDL_CLIGHT;

  /* Read the X,Y data from the temporary file */

  out=fopen("pairspec.tmp","r");
  icte=0;
  ictp=0;
  buffer[0]='#';
  while(buffer[0]=='#') fgets(buffer,BUFSIZ,out);
  while(buffer[0]!='#')
  {
    ex=(double*)realloc(ex,(icte+1)*sizeof(double));
    ey=(double*)realloc(ey,(icte+1)*sizeof(double));
    sscanf(&buffer[11],"%lg%lg",&ex[icte],&ey[icte]);
    icte++;
    fgets(buffer,BUFSIZ,out);
  }
  while(buffer[0]=='#') fgets(buffer,BUFSIZ,out);
  while(buffer[0]!='#')
  {
    px=(double*)realloc(px,(ictp+1)*sizeof(double));
    py=(double*)realloc(py,(ictp+1)*sizeof(double));
    sscanf(&buffer[11],"%lg%lg",&px[ictp],&py[ictp]);
    ictp++;
    if(!(fgets(buffer,BUFSIZ,out))) break;
  }
  fclose(out);

  /***** Plot new results, using cpgplot *****/

  cpgeras();

  /* Start with Magnet limits */
  cpgsci(1);
  cpgsls(1);
  cpgmove(0,0);
  cpgdraw(LMAGNET,0);
  cpgdraw(LMAGNET,HMAGNET);
  cpgdraw(0,HMAGNET);
  cpgdraw(0,0);

  /* The beam line */
  cpgsci(4);
  cpgsls(2);
  m=tan(*theta);
  b=*yi-*xi*m;
  cpgmove(-LMAGNET,m*(-LMAGNET)+b);
  for(bx=-LMAGNET;bx<=2.*LMAGNET;bx=bx+0.05) cpgdraw(bx,m*bx+b);

  /* Now draw in the electron path */
  cpgsci(2);
  cpgsls(1);
  cpgmove(ex[0],ey[0]);
  for(j=1;j<icte;j++)
    cpgdraw(ex[j],ey[j]);

  /* And now the positron path */
  cpgsci(3);
  cpgmove(px[0],py[0]);
  for(j=1;j<ictp;j++)
    cpgdraw(px[j],py[j]);

  /* Mark the point where conversion takes place */
  cpgsci(1);
  cpgpt1(*xi,*yi,8);

  /* Print current status to terminal */

  attron(A_BOLD);
  mvprintw(15,(cols-15)/2,"Current Status:");
  attroff(A_BOLD);
  mvprintw(16,(cols-72)/2,
   "Gamma Ray Energy:  %9.4f MeV     Electron Energy: %9.4f MeV     ",
   *eg/1.602E-13,*ee/1.602E-13);
  mvprintw(17,(cols-72)/2,
   "Conversion site X: %9.4f meters  Beam angle:      %9.4f degrees ",
   *xi,*theta*180./M_PI);
  mvprintw(18,(cols-72)/2,
   "Conversion site Y: %9.4f meters                                 ",
   *yi);
  attron(A_BOLD);
  mvprintw(19,(cols-12)/2,"Step sizes:");
  attroff(A_BOLD);
  mvprintw(20,(cols-76)/2,
   "Energies: %9.4f MeV   Lengths: %9.4f meters   Angles: %9.4f deg.",
   *estepsize/1.602E-13,*lstepsize,*astepsize*180./M_PI);
  for(j=3;j<cols-3;j++) mvprintw(21,j,"-");
  move(rows,cols);
  refresh();

  /* Get control input from user */

  ctrl=getch();
  for(j=1;j<cols-1;j++) mvprintw(23,j," ");
  refresh();
  switch(ctrl)
  {
    case 113: /* Write input file and quit */
      out=fopen("pairspec.in","w");
      if(out)
      {
        fprintf(out,"%13d # Normal run type (verbose output)\n",0);
        fprintf(out,"%13.5E # Gamma ray energy (MeV)\n",*eg/1.602E-13);
        fprintf(out,"%13.5E # Electron energy (MeV)\n",*ee/1.602E-13);
        fprintf(out,"%13.5E # Beam angle (deg)\n",*theta*180./M_PI);
        fprintf(out,"%13.5E # \"Y\" of conversion site (m)\n",*yi);
        fprintf(out,"%13.5E # \"X\" of conversion site (m)\n",*xi);
        fclose(out);
      }
      else
        fprintf(stderr,"Warning: could not write to pairspec.in\n");
      
      *i=*i+1; /* This will cause us to leave the main loop */
      break;

    case 48: /* Revert stepsizes to original values */
      *lstepsize=LSTEPSIZE;
      *estepsize=ESTEPSIZE*1.602E-13;
      *astepsize=ASTEPSIZE*M_PI/180;
      break;

    case 61: /* Double the stepsizes */
      *lstepsize=*lstepsize*2.0;
      *estepsize=*estepsize*2.0;
      *astepsize=*astepsize*2.0;
      break;

    case 45: /* Halve the stepsizes */
      *lstepsize=*lstepsize/2.0;
      *estepsize=*estepsize/2.0;
      *astepsize=*astepsize/2.0;
      break;

    case 339: /* Increase gamma ray energy */
      *eg=*eg+*estepsize;
      break;

    case 338: /* Decrease gamma ray energy */
      *eg=*eg-*estepsize;
      if(*eg<*ee+emin)
      {
        attron(A_BOLD);
        mvprintw(23,(cols-43)/2,"Cannot reduce gamma ray energy any further!");
	attroff(A_BOLD);
	*eg=*ee+emin+0.0001*1.602E-13;
      }
      break;

    case 72: /* Increase electron energy */
      *ee=*ee+*estepsize;
      if(*ee>*eg-emin)
      {
        attron(A_BOLD);
        mvprintw(23,(cols-44)/2,"Cannot increase electron energy any further!");
	attroff(A_BOLD);
	*ee=*eg-emin-0.0001*1.602E-13;
      }
      break;

    case 70: /* Decrease electron energy */
      *ee=*ee-*estepsize;
      if(*ee<emin)
      {
        attron(A_BOLD);
        mvprintw(23,(cols-42)/2,"Cannot reduce electron energy any further!");
	attroff(A_BOLD);
	*ee=emin+0.0001*1.602E-13;
      }
      break;

    case 331: /* Increase beam angle */
      *theta=*theta+*astepsize;
      if(*theta>MAXANGLE*M_PI/180.)
      {
        attron(A_BOLD);
        mvprintw(23,(cols-39)/2,"Cannot increase beam angle any further!");
	attroff(A_BOLD);
	*theta=MAXANGLE*M_PI/180.;
      }
      break;

    case 126: /* Decrease beam angle */
      *theta=*theta-*astepsize;
      if(*theta<-(MAXANGLE*M_PI/180.))
      {
        attron(A_BOLD);
        mvprintw(23,(cols-39)/2,"Cannot decrease beam angle any further!");
	attroff(A_BOLD);
	*theta=-(MAXANGLE*M_PI/180.);
      }
      break;

    case 259: /* Increase Y */
      *yi=*yi+*lstepsize;
      break;

    case 258: /* Decrease Y */
      *yi=*yi-*lstepsize;
      break;

    case 260: /* Decrease X */
      *xi=*xi-*lstepsize;
      break;

    case 261: /* Increase X */
      *xi=*xi+*lstepsize;
      break;
      
    case 103: /* Manually enter gamma ray energy */
      mvprintw(23,3,"Enter new gamma ray energy (in MeV): ");
      echo();
      getstr(buffer);
      noecho();
      *eg=atof(buffer)*1.602E-13;
      if(*eg<*ee+emin)
      {
        for(j=1;j<cols-1;j++) mvprintw(23,j," ");
        attron(A_BOLD);
        mvprintw(23,(cols-43)/2,"Cannot reduce gamma ray energy any further!");
	attroff(A_BOLD);
	*eg=*ee+emin+0.0001*1.602E-13;
      }
      break;

    case 101: /* Manually enter electron energy */
      mvprintw(23,3,"Enter new electron energy (in MeV): ");
      echo();
      getstr(buffer);
      noecho();
      *ee=atof(buffer)*1.602E-13;
      if(*ee>*eg-emin)
      {
        for(j=1;j<cols-1;j++) mvprintw(23,j," ");
        attron(A_BOLD);
        mvprintw(23,(cols-44)/2,"Cannot increase electron energy any further!");
	attroff(A_BOLD);
	*ee=*eg-emin-0.0001*1.602E-13;
      }
      if(*ee<emin)
      {
        for(j=1;j<cols-1;j++) mvprintw(23,j," ");
        attron(A_BOLD);
        mvprintw(23,(cols-42)/2,"Cannot reduce electron energy any further!");
	attroff(A_BOLD);
	*ee=emin+0.0001*1.602E-13;
      }
      break;
      
    case 116: /* Manually enter beam angle */
      mvprintw(23,3,"Enter new beam angle (in degrees): ");
      echo();
      getstr(buffer);
      noecho();
      *theta=atof(buffer)*M_PI/180.;
      if(*theta>MAXANGLE*M_PI/180.)
      {
        for(j=1;j<cols-1;j++) mvprintw(23,j," ");
        attron(A_BOLD);
        mvprintw(23,(cols-39)/2,"Cannot increase beam angle any further!");
	attroff(A_BOLD);
	*theta=MAXANGLE*M_PI/180.;
      }
      if(*theta<-(MAXANGLE*M_PI/180.))
      {
        for(j=1;j<cols-1;j++) mvprintw(23,j," ");
        attron(A_BOLD);
        mvprintw(23,(cols-39)/2,"Cannot decrease beam angle any further!");
	attroff(A_BOLD);
	*theta=-(MAXANGLE*M_PI/180.);
      }
      break;

    case 120: /* Manually enter X */
      mvprintw(23,3,"Enter new X position for conversion (in meters): ");
      echo();
      getstr(buffer);
      noecho();
      *xi=atof(buffer);
      break;
      
    case 121: /* Manually enter Y */
      mvprintw(23,3,"Enter new Y position for conversion (in meters): ");
      echo();
      getstr(buffer);
      noecho();
      *yi=atof(buffer);
      break;

    default: /* Do nothing */
      break;
  }

  refresh();
  *i=*i-1; /* Keeps us from leaving the loop */
  free(ex);
  free(ey);
  free(px);
  free(py);

  return;
}


*** lforce.cpp *** 


#include <rdynlab.h>
#include <pairspec.h>
#include <math.h>

/* This function returns the magnetic field of the silver magnet as a function
 * of position */

double bfield_z(double x, double y, double z)
{
  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;
}

/* Lorentz force, assuming the field above */

RDL_Vector lforce(const RDL_Particle &p, double t)
{
  double field;
  
  field=bfield_z(p.x.x,p.x.y,p.x.z);

  RDL_Vector E=RDL_Vector(0,0,0);
  RDL_Vector B=RDL_Vector(0,0,field);
  RDL_Vector v=p.v;
  RDL_Vector F;

  F=(E*p.q) + ((v^B)*p.q);

  return(F);
}


*** main.cpp *** 


/* pairspec: a program to model particle trajectories in the pair spectrometer.
 * 
 * Built using the RDYNLAB package. Also requires gnuplot as a backend. */

#include <rdynlab.h>
#include <stdio.h>    
#include <pairspec.h> 
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <iostream>
#include <fstream>
#include <cpgplot.h>
#include <ncurses.h>

using namespace std;

int main(int argc, char *argv[])
{
  RDL_Particle e,p;          /* The electron and positron */
  double eg,ee,theta,xi,yi;  /* Control parameters */
  double m,b,bx;             /* Used to draw lines */
  int type;                  /* Run type (see input description) */
  int iter,i,flag;           /* Misc. counters and the output flag for scans */
  int state,statp;           /* Return values for trajectory calculations */
  ofstream output;           /* For directing rdynlab output */
  int pgpid;                 /* PGPLOT device identifier */
  int rows,cols;             /* Number of rows and columns on terminal */
  double lstep,estep,astep;  /* Parameter step sizes for interactive mode */

  /* Check command line arguments */

  if(argc<2)
  {
    fprintf(stderr,"You must specify an input file.\n");
    exit(1);
  }

  /* Read the input file */

  readinput(argv[1],&type,&eg,&ee,&theta,&xi,&yi);

  /* Open pgplot window if in interactive mode */
  
  if(type<0) pgpid=initplot(&rows,&cols,&lstep,&estep,&astep);

  /* Begin scan loop */

  if(type<1) iter=1;
  else
  {
    iter=type;
    srand(time(NULL));
    scanheader();
  }

  for(i=1;i<=iter;i++)
  {

    /* Set random parameters, if this is a scan */
    
    if(type>0) randomize(i,eg,&ee,&theta,&yi,&xi);

    /* Initialize the particles */

    initialize(eg,ee,theta,yi,xi,&e,&p);

    /* Let the particles run! */

    if(type==0)
    {
      fprintf(stdout,"##### ELECTRON #####\n");
      state=RDL_solver(e,lforce,0.0,FSTEP,stopcond,cout);
      fprintf(stdout,"\n\n##### POSITRON #####\n");
      statp=RDL_solver(p,lforce,0.0,FSTEP,stopcond,cout);
    }
    else
    {
      if(type<0) output.open("pairspec.tmp");
      else output.open("/dev/null");
      state=RDL_solver(e,lforce,0.0,FSTEP,stopcond,output);
      statp=RDL_solver(p,lforce,0.0,FSTEP,stopcond,output);
      output.close();
    }

    /* Set the flag value and write out run results, if this is a scan*/

    if(type>0)
    {
      flag=setflag(state,statp,theta,e.v,p.v);
      fprintf(stdout,"%8d %13.5E %13.5E %13.5E %2d X\n",i,
        ee/1.602E-13,theta*180./M_PI,yi,flag);
      fprintf(stderr,"Completed set %8d of %8d\n",i,iter);
    }

    /* If we are in interactive mode, plot the results to the screen */

    else if(type<0) 
      go_interactive(&eg,&ee,&theta,&yi,&xi,&i,rows,cols,&lstep,&estep,&astep);

    /* If we are in single mode, write the beam data to stdout */

    else
    {
      fprintf(stdout,"\n\n# These points define the beamline...\n");
      m=tan(theta);
      b=yi-xi*m;
      for(bx=-LMAGNET;bx<=2.*LMAGNET;bx=bx+0.05)
         fprintf(stdout,"%13.5E %13.5E %13.5E\n",bx,m*bx+b,0.0);
    }

  } /* End loop */

  if(type<0) { cpgclos(); endwin(); }

  return 0;
}


*** pairspec.h *** 


#ifndef _PAIRSPEC_H_
#define _PAIRSPEC_H_

#include <rdynlab.h>

/* Constants */

#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 LONGTIME 3.e-8 /* Maximum iteration time */
#define FSTEP 1.e-12   /* Initial step size */
#define LSTEPSIZE 0.01 /* Initial length step size for interactive mode */
#define ESTEPSIZE 1.00 /* Initial energy step size for interactive mode */
#define ASTEPSIZE 5.00 /* Initial angle step size for interactive mode */

/* Function definitions */

void readinput(const char*,int*,double*,double*,double*,double*,double*);
void randomize(int,double,double*,double*,double*,double*);
void initialize(double,double,double,double,double,RDL_Particle*,RDL_Particle*);
RDL_Vector lforce(const RDL_Particle&,double);
int stopcond(const RDL_Particle&,double,int);
void scanheader(void);
int setflag(int,int,double,RDL_Vector,RDL_Vector);
int initplot(int*,int*,double*,double*,double*);
void go_interactive(double*,double*,double*,double*,double*,int*,int,int,double*,double*,double*);

#endif


*** randomize.cpp *** 


#include <stdlib.h>
#include <pairspec.h>
#include <rdynlab.h>
#include <math.h>
#include <stdio.h>

/* Select random values for control parameters in scan mode */

void randomize(int iter,      /* Iteration number (input) */
               double eg,     /* Gamma ray energy (input) */
	       double *ee,    /* Electron energy (return) */
	       double *theta, /* Beam angle (return) */
	       double *yi,    /* Y position (return) */
	       double *xi)    /* X position (return) */
{
  double angle=MAXANGLE*M_PI/180.;       /* Maximum allowed angle, in radians */
  double ere=EMASS*RDL_CLIGHT*RDL_CLIGHT;/* Rest energy of electron */
  FILE *out;
  char fname[BUFSIZ];

  /* Randomize the electron energy, angle, and height within the appropriate
   * limits */

  *ee=0.0;
  while(*ee<ere || (eg-(*ee))<ere ) /* Kinematically allowed energies only! */
    *ee=(double)rand()/(double)RAND_MAX*eg;
  *theta=(double)rand()/(double)RAND_MAX*angle;
  *yi=(double)rand()/(double)RAND_MAX*HMAGNET;

  /* Write a single mode input file for this run */

  sprintf(fname,"%8.8d.pars",iter);
  out=fopen(fname,"w");
  if(out)
  {
    fprintf(out,"%13d # Normal run type - generates verbose output\n",0);
    fprintf(out,"%13.5E # Gamma ray energy (MeV)\n",eg/1.602E-13);
    fprintf(out,"%13.5E # Electron energy (MeV)\n",*ee/1.602E-13);
    fprintf(out,"%13.5E # Beam angle (deg)\n",*theta*180./M_PI);
    fprintf(out,"%13.5E # \"Y\" of conversion site (m)\n",*yi);
    fprintf(out,"%13.5E # \"X\" of conversion site (m)\n",*xi);
    fclose(out);
  }
  else
    fprintf(stderr,"Warning: could not write file to %s\n",fname);

  /* Return */

  return;
}


*** readinput.cpp *** 


#include <pairspec.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/* Read control parameters from the input file */

void readinput(const char *fname, /* Filename (input) */
               int *type,         /* Run type (return) */
	       double *eg,        /* Gamma ray energy (return) */
	       double *ee,        /* Electron energy (return) */
	       double *theta,     /* Angle of magnet axis (return) */
	       double *xi,        /* X position of production (return) */
	       double *yi)        /* Y position of production (return) */
{
  FILE *f;
  char buffer[BUFSIZ];

  /* Open the file */

  f=fopen(fname,"r");

  /* Check that we were able to open file for reading */

  if(!f)
  {
    fprintf(stderr,"Cannot read file %s\n",fname);
    exit(2);
  }

  /* Read in run type */

  fgets(buffer,BUFSIZ,f);
  *type=atoi(buffer);

  /* Read in the gamma ray energy, in MeV */

  fgets(buffer,BUFSIZ,f);
  *eg=atof(buffer);

  /* If this is not a scan run, read in the electron energy (in MeV), the
   * angle of the beam relative to the magnet axis (in degrees), and the
   * Y position of the production site (in meters) */
  
  if(*type<1)
  {
    fgets(buffer,BUFSIZ,f);
    *ee=atof(buffer);
    fgets(buffer,BUFSIZ,f);
    *theta=atof(buffer);
    fgets(buffer,BUFSIZ,f);
    *yi=atof(buffer);
  }

  /* Read in the X position of the production site (in meters) */

  fgets(buffer,BUFSIZ,f);
  *xi=atof(buffer);

  /* Write the data read in to stdout */

  fprintf(stdout,"# Parameters read from input file %s:\n",fname);

  if(*type>0) fprintf(stdout,"# Run type: scan with %d iterations\n",*type);
  else if(*type==0) fprintf(stdout,"# Run type: single\n");
  else fprintf(stdout,"# Run type: interactive\n");

  fprintf(stdout,"#   Gamma ray energy (MeV):        %g\n",*eg);
  if(*type<1)
  {
    fprintf(stdout,"#   Electron energy (MeV):         %g\n",*ee);
    fprintf(stdout,"#   Beam angle (degrees):          %g\n",*theta);
    fprintf(stdout,"#   Y of production site (meters): %g\n",*yi);
  }
  fprintf(stdout,"#   X of production site (meters): %g\n",*xi);

  /* Convert all input quantities into SI units (RDYNLAB uses SI units) */

  *theta=*theta*M_PI/180.;
  *eg=*eg*1.602E-13;
  *ee=*ee*1.602E-13;

  /* Close the input file and return */

  fclose(f);
  return;
}


*** scanheader.cpp *** 


#include <stdio.h>
#include <pairspec.h>

/* Write out the header for the scan mode output file. This was just too bulky
 * in the main routine */

void scanheader(void)
{
  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: Both particles have exit angles > %g deg.\n",
    (double)GOODANGLE);
  fprintf(stdout,"#       1: Electron has an exit angle > %g deg.\n",
    (double)GOODANGLE);
  fprintf(stdout,"#       2: Positron has an exit angle > %g deg.\n",
    (double)GOODANGLE);
  fprintf(stdout,"#       3: Neither particle has an exit angle > %g deg.\n",
    (double)GOODANGLE);
  fprintf(stdout,"#       4: Electron hits yoke.\n");
  fprintf(stdout,"#       5: Positron hits yoke.\n");
  fprintf(stdout,"#       6: Both particles hit yoke.\n");
  fprintf(stdout,"#       7: Electron is trapped in field.\n");
  fprintf(stdout,"#       8: Positron is trapped in field.\n");
  fprintf(stdout,"#       9: Both particles are trapped in field.\n");
  return;
}


*** setflag.cpp *** 


#include <rdynlab.h>
#include <pairspec.h>
#include <math.h>

/* Set the output flag */

int setflag(int state,     /* State of electron, as returned by stopcond */
            int statp,     /* State of positron, as returned by stopcond */
	    double theta,  /* Beam angle */
	    RDL_Vector ev, /* Final velocity of electron */
	    RDL_Vector pv) /* Final velocity of positron */
{
  int flag;
  double phie,phip;
  double garad;
  RDL_Vector beam;

  flag=0;

  /* Notes: if return values state/statep are 1, particle is trapped
   * in field. If return values are 2, particle hits yoke. If return
   * values are 3, then both particles leave the field successfully. */
  
  while(2>1) /* Allows me to break once flag is set */
  {
    /* Check for yoke impacts */
    if(state==2 && statp!=2) { flag=4; break; }
    if(state!=2 && statp==2) { flag=5; break; }
    if(state==2 && statp==2) { flag=6; break; }
    /* Check for trapped particles */
    if(state==1 && statp!=1) { flag=7; break; }
    if(state!=1 && statp==1) { flag=8; break; }
    if(state==1 && statp==1) { flag=9; break; }

    /* Calculate exit angles from final velocities */

    beam.x=cos(theta);
    beam.y=sin(theta);
    phie=acos((ev*beam)/(length(ev)*length(beam)));
    phip=acos((pv*beam)/(length(pv)*length(beam)));

    garad=GOODANGLE*M_PI/180.;
    if(phie>=garad && phip>=garad) { flag=0; break; }
    if(phie>=garad && phip<garad)  { flag=1; break; }
    if(phie<garad && phip>=garad)  { flag=2; break; }
    
    flag=3;
    break;
  }

  return flag;
}


*** stopcond.cpp *** 


#include <rdynlab.h>
#include <pairspec.h>

/* Stopping conditions for calculating particle trajectories */

int stopcond(const RDL_Particle &p, double t,int iter)
{
  double m,b,x0;

  /* If the particle is still hanging around after this long, then it's 
   * trapped in the field */

  if(t>LONGTIME) return 1;

  /* Check to see if the particle has hit the yoke */

  if(p.v.y<0 && p.x.y<0) /* Only happens if y velocity is negative 
                          * and particle is below y=0 */
  {
    m=p.v.y/p.v.x;
    b=p.x.y-m*p.x.x;
    x0=-b/m;
    if(x0>=0.0 && x0<=LMAGNET) return 2;
  }

  /* Check to see if the particle has left some arbitrary bounds */

  if(p.x.x>2.0*LMAGNET || p.x.x<-LMAGNET || p.x.y>2.0*HMAGNET || p.x.y<-HMAGNET)
    return 3;

  /* Stopping conditions have not been met. Return zero */

  return 0;
}
Retrieved from "https://wiki.iac.isu.edu/index.php?title=New_Particle_Trajectory_Code&oldid=44381"