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