Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Bateman Equation in C++

  1. Jun 13, 2007 #1

    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 [Broken]

    Code (Text):

    #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
        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 cant be j
                    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!
    Last edited by a moderator: May 2, 2017
  2. jcsd
  3. Jun 13, 2007 #2


    User Avatar
    Staff Emeritus
    Science Advisor
    Gold Member

    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
  4. Jun 13, 2007 #3
    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.
  5. Nov 1, 2007 #4
    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?
  6. Nov 1, 2007 #5
    have you checked all of the values prior to use?
    If so, are you trying to calculate a null with an int/double?
  7. Dec 3, 2010 #6
    would it be possible to get a copy of the source code?
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook