Hi everybody,

I'm writing a program to evaluate Dawson's integral and use it to calculate some parameters for a theoretical particle beam. I'm using the function Dawson from Numerical recipes in C (2nd edition, with the float types chaged to doubles for better precision.

The problem I have is that the program is not giving me the expected answer. I have a benchmark model, but the values don't match. Any help would be appreciated.

And here's the benchmark: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); }

Initial radius: 0.5mm

Distance to homocentric point: 100mm

Initial current: 10nA

r_w: 0.0769mm

d_w: 127.33mm

Any help would be appreciated.

Thanks,

lordbubonicus