//
// Program: ndecay.cxx
//
// Description: solve the nuclear decay equation using the Euler
//              numerical method
//
// Author: A. Garfagnini
//
// Date: 1 October 2012
//

#include <iostream>
#include <fstream>
#include <unistd.h>
#include <libgen.h>


using namespace std;

void usage(const char * pname)
{
    cerr << "Usage: " << pname << " [-o file_name]" << endl;
}

int main(int argc, char ** argv)
{
   int opt;
   ofstream out_file;

    // Handle command line options
    while((opt = getopt(argc, argv, "o:h?")) != -1) {

        switch (opt) {

        case 'o':
            out_file.open(optarg);
            break;

        case 'h':
        case '?':
        default:
            usage( basename( argv[0] ) );
            return 1;
        }
    }


    // Handle user data input
    cout << "Number of Nuclei at t0: ";
    int n0;
    cin >> n0;
    if (out_file.is_open()) {
        out_file << "# Number of Nuclei at t0: " << n0 << endl;
    }

    cout << "Decay constant lambda: ";
    double lambda;
    cin >> lambda;
    if (out_file.is_open()) {
        out_file << "# Decay constant lambda: " << lambda << endl;
    }

    cout << "Time step h: ";
    double h;
    cin >> h;
    if (out_file.is_open()) {
        out_file << "Time step h: " << h << endl;
    }

    cout << "T_final: ";
    double t_f;
    cin >> t_f;
    if (out_file.is_open()) {
        out_file << "# T_final: " << t_f << endl;
    }


    // -- Evolve the system
    double t = 0.0;
    double ndecay = n0;

    if (out_file.is_open()) {
        out_file << t << " " << ndecay << endl;
    } else {
        cout << t << " " << ndecay << endl;
    }

    while (t < t_f) {

        t += h;
        ndecay = ndecay * (1.0 - lambda*h);

        if (out_file.is_open()) {
            out_file << t << " " << ndecay << endl;
        } else {
            cout << t << " " << ndecay << endl;
        }

    }

    if (out_file.is_open()) {
        out_file.close();
    }

    return 0;
}