Problem with Dawson's integral

This is a discussion on Problem with Dawson's integral within the C Programming forums, part of the General Programming Boards category; Hi everybody, I'm writing a program to evaluate Dawson's integral and use it to calculate some parameters for a theoretical ...

  1. #1
    Registered User
    Join Date
    Oct 2007
    Posts
    16

    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
    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:
    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

  2. #2
    and the hat of wrongness Salem's Avatar
    Join Date
    Aug 2001
    Location
    The edge of the known universe
    Posts
    32,451
    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);
    If you dance barefoot on the broken glass of undefined behaviour, you've got to expect the occasional cut.
    If at first you don't succeed, try writing your phone number on the exam paper.
    I support http://www.ukip.org/ as the first necessary step to a free Europe.

  3. #3
    Registered User
    Join Date
    Oct 2007
    Posts
    16
    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. #4
    and the hat of wrongness Salem's Avatar
    Join Date
    Aug 2001
    Location
    The edge of the known universe
    Posts
    32,451
    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[1]=0.852144
    count[2]=0.236928
    count[3]=0.018316
    count[4]=0.000394
    count[5]=0.000002
    count[6]=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.
    If you dance barefoot on the broken glass of undefined behaviour, you've got to expect the occasional cut.
    If at first you don't succeed, try writing your phone number on the exam paper.
    I support http://www.ukip.org/ as the first necessary step to a free Europe.

  5. #5
    Registered User
    Join Date
    Oct 2007
    Posts
    16

    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. #6
    Registered User
    Join Date
    Mar 2004
    Posts
    536
    Quote Originally Posted by lordbubonicus View Post
    ...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.
    Last edited by Dave Evans; 11-25-2007 at 03:05 PM.

Popular pages Recent additions subscribe to a feed

Similar Threads

  1. Memory problem with Borland C 3.1
    By AZ1699 in forum C Programming
    Replies: 16
    Last Post: 11-16-2007, 10:22 AM
  2. Problem with for loop calling external function
    By lordbubonicus in forum C Programming
    Replies: 2
    Last Post: 10-13-2007, 10:54 AM
  3. Someone having same problem with Code Block?
    By ofayto in forum C++ Programming
    Replies: 1
    Last Post: 07-12-2007, 08:38 AM
  4. A question related to strcmp
    By meili100 in forum C++ Programming
    Replies: 6
    Last Post: 07-07-2007, 02:51 PM
  5. WS_POPUP, continuation of old problem
    By blurrymadness in forum Windows Programming
    Replies: 1
    Last Post: 04-20-2007, 06:54 PM

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21