Difference between revisions of "New Particle Trajectory Code"

From New IAC Wiki
Jump to navigation Jump to search
Line 520: Line 520:
 
----
 
----
  
−
'''*** NAME ***'''
+
'''*** lforce.cpp ***'''
  
 
<pre><nowiki>
 
<pre><nowiki>
  
 +
#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);
 +
}
  
 
</nowiki></pre>
 
</nowiki></pre>
Line 530: Line 570:
 
----
 
----
  
−
'''*** NAME ***'''
+
'''*** main.cpp ***'''
  
 
<pre><nowiki>
 
<pre><nowiki>
  
 +
/* 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;
 +
}
  
 
</nowiki></pre>
 
</nowiki></pre>
Line 540: Line 697:
 
----
 
----
  
−
'''*** NAME ***'''
+
'''*** pairspec.h ***'''
  
 
<pre><nowiki>
 
<pre><nowiki>
  
 +
#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
  
 
</nowiki></pre>
 
</nowiki></pre>
Line 550: Line 739:
 
----
 
----
  
−
'''*** NAME ***'''
+
'''*** randomize.cpp ***'''
  
 
<pre><nowiki>
 
<pre><nowiki>
  
 +
#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;
 +
}
  
 
</nowiki></pre>
 
</nowiki></pre>

Revision as of 18:14, 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;
}


*** NAME ***


*** NAME ***


*** NAME ***


*** NAME ***
Retrieved from "https://wiki.iac.isu.edu/index.php?title=New_Particle_Trajectory_Code&oldid=44380"