PROGRAM nuclear_decay
! simulation of decay of unstable nuclei
DIM ncum(0 to 1000)
CALL initial(N0,p,ntrial,tmax,ncum())
CALL decay(N0,p,ntrial,tmax,ncum())
CALL output(ntrial,tmax,ncum())
END

SUB initial(N0,p,ntrial,tmax,ncum())
    RANDOMIZE
    INPUT prompt "initial number of unstable nuclei = ": N0
    INPUT prompt "decay probability for unstable nucleus = ": p
    INPUT prompt "number of time intervals per trial = ": tmax
    INPUT prompt "number of trials = ": ntrial
    FOR t = 1 to tmax
        LET ncum(t) = 0      ! accumulate number of unstable nuclei
    NEXT t
END SUB

SUB decay(N0,p,ntrial,tmax,ncum())
    DIM N(5000)
    FOR itrial = 1 to ntrial
        FOR i = 1 to N0
            LET N(i) = 1     ! nucleus = 1 if unstable
        NEXT i
        LET ncum(0) = ncum(0) + N0
        LET n_unstable = N0  ! # of unstable nuclei
        FOR t = 1 to tmax
            FOR i = 1 to N0
                IF N(i) = 1 then
                   IF rnd <= p then
                      LET N(i) = 0     ! nucleus decays
                      LET n_unstable = n_unstable - 1
                   END IF
                END IF
            NEXT i
            LET ncum(t) = ncum(t) + n_unstable   ! accumulate data
        NEXT t
    NEXT itrial
END SUB

SUB output(ntrial,tmax,ncum())
    ! print data to file to be read by separate program
    OPEN #1: name "decay.dat", create new, access output
    PRINT #1: "time", "mean number of unstable nuclei"
    FOR t = 0 to tmax
        PRINT #1: t, ncum(t)/ntrial
    NEXT t
    CLOSE #1
END SUB