PROGRAM Radioactive_Decay
!------------------------------------------------------------------
! 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)
!------------------------------------------------------------------
IMPLICIT NONE
INTEGER, PARAMETER :: K = 8
REAL(KIND=K), PARAMETER :: Lambda = 0.1   ! decay constant
REAL(KIND=K), PARAMETER :: dt     = 0.1   ! time step
INTEGER, PARAMETER      :: N0     = 1000  ! initial number of parents
INTEGER                 :: i, N, Ntheory, NN
REAL(KIND=8)            :: t, Tmax, R, Error

 N    = N0
 NN   = N0
 t    = 0.0
 Tmax = 50.0

 ! initiate random number generator
 CALL RANDOM_SEED

 DO

   t = t + dt

   ! Theoratical (expeced) number of undecayed nuclei at time t 
   Ntheory = NINT( N0*EXP(-Lambda*t) )

   ! Loop over each remaining parent nucleus
   DO i=1,NN
     CALL RANDOM_NUMBER(R)

     ! Decide if the nucleus decays:
     IF( R<Lambda*dt ) N=N-1
   END DO
   NN = N

   ! Error from counting operation
   Error   = SQRT(REAL(N))

   PRINT '(F6.2,1x,2I8,1x,F8.3)',t, Ntheory, N, Error

   IF(t>=Tmax) EXIT

 END DO

END PROGRAM