Why is My C++ Program for the Bateman Equation Returning NAN?

[lapo3399]

Hello,

I am attempting to create a program that will evaluate the Bateman equation for radiactive decay series, given several decay constants stored in a text file. I am not sure where this is going wrong, but no matter what chain member I enter in as the desired nuclide to find the mass remaining after 1 million years, it outputs NAN - even when I choose the nuclide that I am certain to have a perceivable value. The equation I am emulating later in the program is shown below, please help to explain why this is not working! All help is appreciated!

http://img515.imageshack.us/img515/5097/screenshotml6.jpg

#include <iostream>
#include <math.h>
#include <fstream>
#include <string>

using namespace std;

int main(int argc, char** argv)
{
	string line; //for file input
	int iter; //stores which nuclide will be calculated for
	
	double dcs[23]; //array for decay constants
	int anum = 1; //counting var
	
	ifstream deconst;
	deconst.open ("decayconstants.txt", ios::in); //file with decay constants
	if (deconst.is_open())
	{
		while (! deconst.eof())
		{
			getline(deconst, line);
			dcs[anum]=strtod(line.c_str(), NULL); //store as double in decay constant array
			anum++;	
		}
		
	}
	
	std::cout << "What nuclide in the chain should be evaluated?";
	std::cin >> iter;
	
        //I separated the equation into multiple parts

	double consts=1; //part one: the multiplied decay constants at the beginning of the equation
	double init = 100; //the intial mass, N(O)
	double time = 1000000; //the variable t, one million years
	double sum=0; // the second half of the equation, the summation
	double totalmass; //the product of consts,init, and sum
		
	for (int i = 1; i < iter; i++)
	{
		consts=consts*dcs[i]; //multiply all the necessary decay constants 
	}
	
	for (int q = 1; q < iter; q++) //the summation
	{
		double expo = exp(-dcs[q]*time); //determine the numerator in the summation
		double denom = 1;
		 
		for (int w = 1; w < iter; w++) //the product
		{
			if (w==iter) //k can't be j
			{
				break;
			}
			else 
			{
				denom = denom * (dcs[w]-dcs[q]); //multiply together to get the product
			}
			
		}
		
		sum = sum + (expo/denom); //add up all of the summations
		
	}
	
	totalmass = sum*consts*init;

	std::cout << totalmass;
	
	return 0;
}

Thanks again!

 

[chroot]

Well, I can't run your program without the file it needs to read, but I'd suggest putting some cout statements throughout the code to print out the intermediate values.

- Warren

 

[lapo3399]

I've attempted that, and I finally realized what the problems where this morning. Firstly, the first for loop should read "i<(iter-1)", and I shouldn't have included the break statement in the other for loop.

 

[mehafroozs]

Hi,
I see that you said your program worked. I am working with a few students to develope a detector but i also needed them to study the decan chains of Radon. I think your program will be useful, do you think i could get a copy of it?

 

[KingMike]

have you checked all of the values prior to use?
If so, are you trying to calculate a null with an int/double?

 

[ustoehlker]

Hi,
would it be possible to get a copy of the source code?