      PROGRAM nuclear_decay
*     simulation of decay of unstable nuclei
      IMPLICIT NONE
      REAL*8 p
      INTEGER N0, ncum(0:1000), ntrial, tmax
      CALL initial(N0,p,ntrial,tmax,ncum)
      CALL decay(N0,p,ntrial,tmax,ncum)
      CALL output(ntrial,tmax,ncum)
      END
C
      SUBROUTINE initial(N0,p,ntrial,tmax,ncum)
      IMPLICIT NONE
      REAL*8 p
      REAL dmy, rnd
      INTEGER t, N0, ntrial, tmax, ncum(0:1000)
      dmy = rnd(-1)
      WRITE(*,'(A,$)') 'initial number of unstable nuclei = '
      READ(*,*) N0
      WRITE(*,'(A,$)') 'decay probability for unstable nucleus = '
      READ(*,*) p
      WRITE(*,'(A,$)') 'number of time intervals per trial = '
      READ(*,*) tmax
      WRITE(*,'(A,$)') 'number of trials = '
      READ(*,*) ntrial
      DO t = 0, tmax
         ncum(t) = 0                        ! accumulate number of unstable nuclei
      END DO
      END
C
      SUBROUTINE decay(N0,p,ntrial,tmax,ncum)
      IMPLICIT NONE
      REAL*8 p
      REAL rnd
      INTEGER N(5000), i, t, itrial, n_unstable, N0, ntrial, tmax, 
     +        ncum(0:1000)
      DO itrial = 1, ntrial
         DO i = 1, N0
            N(i) = 1                        ! nucleus = 1 if unstable
         END DO
         ncum(0) = ncum(0) + N0
         n_unstable = N0                    ! # of unstable nuclei
         DO t = 1, tmax
            DO i = 1, N0
               IF(N(i) .EQ. 1) THEN
                  IF(rnd(0) .LE. p) THEN
                     N(i) = 0               ! nucleus decays
                     n_unstable = n_unstable - 1
                  END IF
               END IF
            END DO
            ncum(t) = ncum(t) + n_unstable  ! accumulate data
         END DO
      END DO
      END
C
      SUBROUTINE output(ntrial,tmax,ncum)
*     print data to file to be read by separate program
      IMPLICIT NONE
      INTEGER IO1_, t, ntrial, tmax, ncum(0:1000)
      DATA IO1_/11/
      OPEN(IO1_, file = "decay.dat")
      WRITE(IO1_,*) '       time  mean number of unstable nuclei'
      DO t = 0, tmax
         WRITE(IO1_,'(I12,F13.6)') t, DBLE(ncum(t))/ntrial
      END DO
      CLOSE(IO1_)
      END
