New Particle Trajectory Code
Jump to navigation
Jump to search
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;
}
*** NAME ***
*** NAME ***
*** NAME ***
*** NAME ***
*** NAME ***
*** NAME ***
*** NAME ***
*** NAME ***