A Monte Carlo Simulation of Radioactive Decay

This is a truly random proceses. The probability of decay is constant. The probabilty that a nucleus undergoes radioactive decay in time dt is p:

p = λ dt    (for λ dt <<1)

where λ (decay constant) is probability per unit time for the decay of each nucleus of a given nuclide.

Consider a system initially having N0 unstable nuclei. How does the number of parent nuclei, N, change with time?

Theoretically, the number of undecayed nuclei at time t is given by:

N = N 0 exp (- λ t )

where N0 is the number of parent(undecayed) nuclei at t = 0.

 Alogrithm: ``` Determine N0 (initial number of parents) Determine λ (decay constant) Determine Tmax (any time) Determine dt (time step) N = N0 LOOP from t=dt to Tmax, step dt LOOP over each remaining parent nucleus # Generate a random number R between [0,1] from a uniform distribution R = uniform_RandomNumber() # Decide if the nucleus decays: # reduce the number of parents by 1 IF( R < λ dt ) N = N-1 END LOOP over nuclei # calculate theoratical value Ntheory = N0 * exp(-λt) Record t, N and Ntheory END LOOP over time ```
rd-Mo99_to_Tc99.f90 