Code:
/*prj03
name: David Brown
student number: 0500830
usercode: phufai
This program calculates the beam waist radius, and the distance from the aperture to the beam waist for a
homocentric charged particle beam. The initial beam radius, distance from the aperture to the homocentric point, and
total beam current are input by the user and validated using a previously developed bulletproof input validation
program. The inputs are converted into SI units, and then used to calculate several variables. A series approximation
to Dawson's integral from Numerical Recipes in C is used to help find the distance to the beam waist.*/
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
/*Define the input parameters that are not specified by the user*/
#define PARTMASS 133
#define KE 500
/*Define the constants needed for the main function. Values taken from CODATA Recommended values of the Fundamental
Physical Constants, 2002; Review of Modern Physics 77, 1-107 (2005)*/
#define PI 4.0*atan(1.0) /*I thought M_PI was part of math.h, but got compiler errors when I tried to use it*/
#define EPSILON 8.854187817e-12
#define ELEC_CHARGE 1.60217653e-19
/*Define the constants and macros needed for the computation of Dawson's integral*/
#define NMAX 6
#define H 0.4
#define SQR(a) ((sqrarg=(a)) == 0.0?0.0 : sqrarg*sqrarg)
#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
/*function for computing dawson's algorithm is adapted from: W.H.Press, S.A.Teukolsky, W.T.Vetterling and B.P.Flannery,
Numerical Recipes in C: 2nd edition, Cambridge University Press(1992)*/
double DawsonIntegral(double x)
{
int i, n0;
double d1, d2, e1, e2, sum, x2, xp, xx, ans;
static double sqrarg;
static double count[NMAX+1];
static int init=0;
if(init==0)
{
init=1;
for(i=1; i<=NMAX; i++)
count[i]=exp(-SQR((2.0*i-1.0)*H));
}
if(fabs(x)<2.0)
{
x2=x*x;
ans=x*(1.0-(2.0/3.0)*x2*(1.0-(2.0/5.0)*x2*(1.0-(2.0/7.0)*x2)));
}
else
{
xx=fabs(x);
n0=2*(int)(0.5*xx/H+0.5);
xp=xx-n0*H;
e1=exp(2.0*xp*H);
e2=e1*e1;
d1=n0+1;
d2=d1-2.0;
sum=0.0;
for(i=1; i<=NMAX; i++)
{
d1=+2.0;
d2-=2.0;
e1*=e2;
sum+=count[i]*(e1/d1+1.0/(d2*e1));
}
ans=0.5641895835*SIGN(exp(-xp*xp), x)*sum;
}
return(ans);
}
/*I have used my own "bulletproof" input for this program. Code for the validation is adapted from code found at
http://faq.cprogramming.com/cgi-bin/smartfaq.cgi?answer=1043372399&id=1043284385*/
double Validation(int j)
{
int counter=0;
double value;
char buffer[BUFSIZ], *ptr;
if(j==1)
printf("Please enter the value of the initial beam radius in mm: ");
else if(j==2)
printf("Please enter the value of the distance to the homocentric point in mm: ");
else if(j==3)
printf("Please enter the value of the total beam current in nA: ");
while(counter==0)
{
if(fgets(buffer, sizeof(buffer), stdin) != NULL)
{
value=(float)strtod(buffer, &ptr);
if(buffer[0]!='\n' && (*ptr=='\n' || *ptr=='\0'))
{
/*In this implementation of the bulletproof input there is no need to check the magnitude of the
inputs against type limits, as the parameters all have specified limits already*/
if(j==1)
{
if(value<0.1)
printf("Value too small. Please re-enter the value: ");
else if(value>2.0)
printf("Value too large. Please re-enter the value: ");
else
counter=1;
}
if(j==2)
{
if(value<10.0)
printf("Value too small. Please re-enter the value: ");
else if(value>500.0)
printf("Value too large. Please re-enter the value: ");
else
counter=1;
}
if(j==3)
{
if(value<1)
printf("Value too small. Please re-enter the value: ");
else if(value>1.0e3)
printf("Value too large. Please re-enter the value: ");
else
counter=1;
}
}
else
printf("Invalid entry. Please re-enter the value: ");
}
}
/*recast value to double for use in calculations*/
value=(double)value;
return(value);
}
int main()
{
int j;
double r_i, d_h, I, m, E, r_i_dot;
double k, v, x, dawson;
double d_w, r_w;
/*Ask user for values for the unspecified parameters, validate them, and convert them to SI units one at a time*/
j=1;
r_i=1.0e-3*Validation(j);
printf("Initial beam radius is %g m\n", r_i); /*Limit the number of decimal places displayed.*/
j=2;
d_h=1.0e-3*Validation(j);
printf("The distance to the homocentric point is %g m\n", d_h);
j=3;
I=1.0e-9*Validation(j);
printf("The initial beam current is %g A\n", I);
/*Convert the defined input values into SI units. Conversion factors taken from
http://www1.bipm.org/en/si/si_brochure/chapter4/table7.html*/
m=1.66053886e-27*PARTMASS;
printf("The mass of the charged particles is %g kg\n", m);
E=1.60217653e-19*KE;
printf("The kinetic energy of the particles is %g J\n", E);
/*Now that all of the input values have been provided and converted, the calculations can be carried out*/
/*firstly, calculate the axial velocity*/
v=sqrt((2*E)/m);
printf("The axial velocity is %g m/s\n", v);
/*Now calculate the value of k*/
k=(PI*EPSILON*v*m)/(I*ELEC_CHARGE);
printf("k has a value of: %g\n", k);
/*Compute the radial velocity at the edge of the beam. Assume that both this and the axial velocity are
components of particle velocity in the direction of the beam (blue line in assignment brief). tan(theta)=r_i/d_h
v=velocity*cos(theta), and r_i_dot=velocity*sin(theta). Rearranging and substituting leads r_i_dot=v*tan(theta).
Therefore arrive at following formula for the radial velocity at the edge of the beam*/
r_i_dot=(r_i*v)/d_h;
/*Compute the value of x for use in the Dawson's integral function, and in calculating r_w*/
x=sqrt(k)*r_i_dot;
printf("x has a value of %g\n", x);
/*Can now compute the value of the beam waist radius*/
r_w=r_i*exp(-x*x);
printf("initial r_w: %g\n", r_w);
/*calculate the value of dawson's integral using the external function*/
dawson=DawsonIntegral(x);
printf("The value of Dawson's integral is: %g\n", dawson);
d_w=2*v*sqrt(k)*r_i*dawson;
printf("initial d_w: %g\n", d_w);
/*Reconvert the answers into mm before outputting the values*/
r_w=1e3*r_w;
d_w=1e3*d_w;
printf("The radius of the beam waist is %g mm\n", r_w);
printf("The distance from the aperture to the waist is %g mm\n", d_w);
return(0);
}
And here's the benchmark: