# Thread: Problem with Dawson's integral

1. ## Problem with Dawson's integral

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.

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
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!='\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:
Distance to homocentric point: 100mm
Initial current: 10nA

r_w: 0.0769mm
d_w: 127.33mm

Any help would be appreciated.

Thanks,
lordbubonicus 2. I get this, compiled with gcc
Code:
```\$ gcc foo.c
\$ ./a.exe
Please enter the value of the initial beam radius in mm: 0.5
Initial beam radius is 0.0005 m
Please enter the value of the distance to the homocentric point in mm: 100
The distance to the homocentric point is 0.1 m
Please enter the value of the total beam current in nA: 10
The initial beam current is 1e-08 A
The mass of the charged particles is 2.20852e-25 kg
The kinetic energy of the particles is 8.01088e-17 J
The axial velocity is 26934.2 m/s
k has a value of: 0.000103275
x has a value of 1.36858
initial r_w: 7.68291e-05
The value of Dawson's integral is: 0.254825
initial d_w: 0.06975
The radius of the beam waist is 0.0768291 mm
The distance from the aperture to the waist is 69.75 mm```
So how exactly does that differ from your expected result?

Remember that floats and doubles are approximations, so don't expect mathematical operations (say (1/3)*3, which is of course 1), to have the same computational result.

> value=(float)strtod(buffer, &ptr);
> value=(double)value;
I can't see what purpose these casts perform, except to perhaps remove some intermediate precision.

> j=1;
> r_i=1.0e-3*Validation(j);
Save a variable, and just do
r_i=1.0e-3*Validation(1); 3. The difference from the expected result is in the value of d_w - the value I obtain of 69.75 is just over half of the expected value of 127.33. I'm not concerned about the differenc in the value of r_w, as I'm aware of the nature of the approximations involved. But thank you for pointing it out. 4. Well I added some code to print the count array, and each value of sum, and got this.
Code:
```\$ ./a.exe
Please enter the value of the initial beam radius in mm: 0.5
Initial beam radius is 0.0005 m
Please enter the value of the distance to the homocentric point in mm: 100
The distance to the homocentric point is 0.1 m
Please enter the value of the total beam current in nA: 10
The initial beam current is 1e-08 A
The mass of the charged particles is 2.20852e-25 kg
The kinetic energy of the particles is 8.01088e-17 J
The axial velocity is 26934.2 m/s
k has a value of: 0.000103275
x has a value of 1.36858
initial r_w: 7.68291e-05
count=0.852144
count=0.236928
count=0.018316
count=0.000394
count=0.000002
count=0.000000
The value of Dawson's integral is: 0.254825
initial d_w: 0.06975
The radius of the beam waist is 0.0768291 mm
The distance from the aperture to the waist is 69.75 mm```
From which I deduce that the calculation is really the result of this expression.
ans=x*(1.0-(2.0/3.0)*x2*(1.0-(2.0/5.0)*x2*(1.0-(2.0/7.0)*x2)));
One thing I notice is that it doesn't look like a traditional power series (but I'm not a mathematician) as each term is x2 (not x2, x4, x6 etc).

> 0.5641895835
I also notice that the short-form calculation doesn't apply this "approximately a half" adjustment to the answer, but again, maybe that's OK too. 5. ## Problem solved

Thanks very much for the help Salem.

In the end I used a different algortihm for the integral, but still had problems with the short form. So I forced the program to compute the answer using the long form every time, which seemed to work. 6. Originally Posted by lordbubonicus ...I used a different algortihm...
Isn't it better to get to the bottom of things?

1. The following leads to a very poor approximation for your value of x. Shouldn't it be something like if (fabs(x) < 0.1) or (fabs(x) < 0.2)?
Code:
```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)));
}```
(You do know what your value of x is, right?)

2. The following is incorrect for two reasons:
Code:
```	for(i=1; i<=NMAX; i++)
{
d1=+2.0;
d2-=2.0;
e1*=e2;
sum+=count[i]*(e1/d1+1.0/(d2*e1));
}```
Obviously it should be d1 += 2.0, not d1=+2.0
Also, the changes to e1, d1, and d2 should be made after the term is added into sum

With these changes, I get 127.326 for your benchmark comparison.

Finally, if you want to get better precision, instead of just changing floats to doubles, shouldn't you do some error analysis to see whether your approximations require more terms? (Truncation error versus roundoff error.)

D. Popular pages Recent additions 