#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>

//------------------------------------------------------------------
// Monte Carlo Simulation of Radioactive Decay
// 
// This program outputs time elaped vs 
// number of parent nuclei.
//
// lambda  : decay constant (1/s)
// n0      : initial number of undecayed nucleus
// nTheory : expected number of remaining nucleus at time t
// dt      : time step (s)
// tmax    : any time you want to stop the decay(s)
//------------------------------------------------------------------

int main(){

 const double   lambda = 0.1;   // decay constant
 const double   dt     = 0.1;   // time step
 const int      n0     = 1000;  // initial number of parents
 int            i, n, nn, nTheory;
 double         t, tmax, r, error;
 int            nint(double);

 n    = n0;
 nn   = n0;
 t    = 0.0;
 tmax = 50.0;

 // initiate random number generator
 srand(time(NULL));

 while(t<=tmax)
 {
    t += dt;

    // Theoratical (expeced) number of undecayed nuclei at time t 
    nTheory = nint( n0*exp(-lambda*t) );

    // Loop over each remaining parent nucleus
    for(i=1; i<=nn; i++) {
       r = double(rand())/RAND_MAX;

       // Decide if the nucleus decays
       if( r<lambda*dt ) n--;
    }

    nn = n;

    // Error from counting operation
    error = sqrt(n);

    printf("%6.2lf %8d %8d %8.3lf\n",t, nTheory, n, error);
 }

return 0;

} /* main */

// Returns nearest int of a double number x
int nint(double x)
{
  int    nearest_int = int(x);
  double real_part   = x - int(x);

  if(real_part >= 0.5) nearest_int++;

  return nearest_int;
}