#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)
//------------------------------------------------------------------

#define dt     0.2      // time step
#define tmax   10       // maximum time
#define N      50       // number of elements (i.e. tmax/dt)
#define n0     1000     // initial number of parents
#define lambda 0.4      // decay constant

void radioactiveDecay(){

 Int_t      i, j, n;
 Double_t   nn[N], nTheory[N], t[N], error[N], zero[N], r;

 n = nn[0] = nTheory[0] = n0;
 error[0] = sqrt(n0);
 zero[0] = 0.0;

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

   for(j=1; j<N; j++)
   {
      t[j] = j*dt;
      zero[j] = 0.0;

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

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

         // Decide if the nucleus decays
         if( r<lambda*dt ) n--;
      }
      nn[j] = n;

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

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

   TCanvas *c1  = new TCanvas("c1","Radioactive Decay",10,10,800,500);
   TGraph  *gr1 = new TGraphErrors(N, t, nn, zero, error);
   TF1     *gr2 = new TF1("Theory","1000.0*exp(-0.4*x)",0,10);

   gr1->SetMarkerColor(1);
   gr1->SetMarkerStyle(20);
   gr1->SetMarkerSize(0.8);
   gr1->SetTitle("Radioactive Decay");
   gr1->GetXaxis()->SetTitle("time (s)");
   gr1->GetYaxis()->SetTitle("N");

   gr1->Draw("AP");
   gr2->Draw("SAME");
}