DECLARE SUB GetPointSize ()
DECLARE SUB DrawPoint (xc!, yc!, colour!, s!)
DECLARE SUB PrintKeys ()
DECLARE SUB TrialSetup ()
DECLARE SUB TrialRun ()
DECLARE SUB InitAtoms ()
DECLARE SUB DrawAtomBox ()
DECLARE SUB DrawGraphBox ()
DECLARE SUB PlotPoint ()
DECLARE SUB PrintSummary ()
DECLARE SUB GetKey ()
DECLARE SUB GetRandomVars ()
DECLARE SUB GetHalfLife ()
DECLARE SUB GetTime ()
DECLARE SUB GetAtoms ()
DECLARE SUB PrintInstructions ()
DECLARE FUNCTION irand! (i1!, i2!)
SCREEN 12
RANDOMIZE TIMER

'''''constants

'limits
DIM SHARED maxatoms, maxseconds
maxatoms = 10000
maxseconds = 1000

'colours
DIM SHARED radioactive, decayed, white, half
radioactive = 4 'red
decayed = 2 'green
white = 15
half = 3 'blue

'screen locations
DIM SHARED symax
DIM SHARED pxmin, pxmax, pymin, pymax
DIM SHARED gxmin, gxmax, gymin, gymax, gxdiff, gydiff
DIM SHARED xloc, yloc, dxloc
symax = 480
pxmin = 0
pxmax = 440
pymin = 280
pymax = 480
gxmin = 0
gxmax = 440
gymin = 50
gymax = 250
gxdiff = gxmax - gxmin
gydiff = gymax - gymin
xloc = 57
yloc = 16
dxloc = 10

'''''variables

'user-changeable variables
DIM SHARED natoms, nseconds, halflife, time, atoms, pointsize
natoms = 100
nseconds = 500
halflife = 100
pointsize = 1

'program-changeable variables
DIM SHARED px(maxatoms), py(maxatoms), state(maxatoms)
DIM SHARED fraction
DIM SHARED a$

'''''main loop
DO

CALL TrialSetup
CALL GetKey
       
SELECT CASE a$
CASE CHR$(13) 'enter
        CALL TrialRun
        CALL GetKey
CASE "a"
        CALL GetAtoms
CASE "t"
        CALL GetTime
CASE "h"
        CALL GetHalfLife
CASE "s"
        CALL GetPointSize
CASE "r"
        CALL GetRandomVars
CASE CHR$(27) 'esc
        SYSTEM
END SELECT

LOOP

SUB DrawAtomBox

'print atom box title
LOCATE 1, 1
PRINT "Atoms in a Radioactive Substance"

'draw atom box
LINE (pxmin, symax - pymin)-(pxmax, symax - pymax), white, B

'draw atoms
FOR p = 1 TO natoms
        CALL DrawPoint(px(p), py(p), radioactive, pointsize)
NEXT p

END SUB

SUB DrawGraphBox

'print graph box title
LOCATE 14, 1
PRINT "Graph of ";
COLOR radioactive
PRINT "Radioactive";
COLOR white
PRINT " and ";
COLOR decayed
PRINT "Decayed";
COLOR white
PRINT " Atoms versus Time"

'draw graph box
LINE (gxmin, symax - gymin)-(gxmax, symax - gymax), white, B

'draw lines where halflife should occur
IF halflife <= nseconds THEN
        xhalf = gxmin + gxdiff * halflife / nseconds
        yhalf = symax - (gymin + gydiff / 2)
        LINE (xhalf, symax - gymin)-(xhalf, yhalf), half
        LINE (gxmin, yhalf)-(xhalf, yhalf), half
END IF

END SUB

SUB DrawPoint (xc, yc, colour, s)
        LINE (xc - s, symax - yc - s)-(xc + s, symax - yc + s), colour, BF
END SUB

SUB GetAtoms
        DO
                LOCATE yloc + 5, xloc + dxloc
                INPUT natoms
                IF natoms < 1 OR natoms > maxatoms THEN
                        PRINT "Number of atoms must be 1 to "; maxatoms
                END IF
        LOOP WHILE natoms < 1 OR natoms > maxatoms
END SUB

SUB GetHalfLife
        DO
                LOCATE yloc + 9, xloc + dxloc
                INPUT halflife
                IF halflife < 1 OR halflife > maxseconds THEN
                        PRINT "Half life must be 1 to "; maxseconds
                END IF
        LOOP WHILE halflife < 1 OR halflife > maxseconds
END SUB

SUB GetKey
a$ = ""
DO
        a$ = INKEY$
LOOP WHILE a$ = ""
END SUB

SUB GetPointSize
        DO
                LOCATE 14, xloc
                PRINT "Point Size"
                LOCATE 14, xloc + dxloc
                INPUT pointsize
                IF pointsize < 0 OR pointsize > 5 THEN
                        PRINT "Point size must be 0 to 5."
                END IF
        LOOP WHILE pointsize < 0 OR pointsize > 5
END SUB

SUB GetRandomVars
        texp = 10 ^ (INT(RND * 3))
        natoms = (INT(RND * 9) + 1) * 10 ^ (INT(RND * 4))
        nseconds = (INT(RND * 9) + 1) * texp
        halflife = (INT(RND * 9) + 1) * texp
END SUB

SUB GetTime
        DO
                LOCATE yloc + 11, xloc + dxloc
                INPUT nseconds
                IF nseconds < 1 OR nseconds > maxseconds THEN
                        PRINT "Time must be 1 to "; maxseconds
                END IF
        LOOP WHILE nseconds < 1 OR nseconds > maxseconds
END SUB

SUB InitAtoms
FOR p = 1 TO natoms
        px(p) = irand(pxmin + 1, pxmax - 1)
        py(p) = irand(pymin + 1, pymax - 1)
        state(p) = radioactive
NEXT p
END SUB

FUNCTION irand (i1, i2)
irand = INT(RND * (i2 - i1 + 1)) + i1
END FUNCTION

SUB PlotPoint

gx = gxmin + gxdiff * time / nseconds
gy = gymin + gydiff * atoms / natoms
CALL DrawPoint(gx, gy, radioactive, pointsize)
CALL DrawPoint(gx, gymin + gymax - gy, decayed, pointsize)

END SUB

SUB PrintKeys

LOCATE 1, xloc
PRINT "Half-Life Experiment"

LOCATE 3, xloc
PRINT "Key   Action"
LOCATE 4, xloc
PRINT "---   ------"
LOCATE 5, xloc
PRINT "A     Change Atoms"
LOCATE 6, xloc
PRINT "T     Change Time"
LOCATE 7, xloc
PRINT "H     Change Half-life"
LOCATE 8, xloc
PRINT "S     Change Point Size"
LOCATE 9, xloc
PRINT "R     Random Changes"
LOCATE 10, xloc
PRINT "Enter Start"
LOCATE 11, xloc
PRINT "P     Pause/Resume"
LOCATE 12, xloc
PRINT "Esc   Stop/Exit"

END SUB

SUB PrintSummary

LOCATE yloc + 1, xloc
PRINT "Atoms (y-axis)"
COLOR radioactive
LOCATE yloc + 2, xloc
PRINT "Radioactive"; atoms
COLOR half
LOCATE yloc + 3, xloc
PRINT "Half total "; natoms / 2
COLOR decayed
LOCATE yloc + 4, xloc
PRINT "Decayed    "; natoms - atoms
COLOR white
LOCATE yloc + 5, xloc
PRINT "Total      "; natoms

LOCATE yloc + 7, xloc
PRINT "Time (x-axis)"
LOCATE yloc + 8, xloc
PRINT "Remaining  "; nseconds - time
COLOR half
LOCATE yloc + 9, xloc
PRINT "Half-life  "; halflife
COLOR white
LOCATE yloc + 10, xloc
PRINT "Elapsed    "; time
LOCATE yloc + 11, xloc
PRINT "Total      "; nseconds

END SUB

SUB TrialRun

a$ = ""

FOR time = 1 TO nseconds

a$ = INKEY$
IF a$ = CHR$(27) THEN EXIT FOR 'esc
IF a$ = "p" THEN
        a$ = ""
        DO
                a$ = INKEY$
        LOOP WHILE a$ <> "p"
END IF

FOR p = 1 TO natoms
        IF state(p) <> decayed AND RND < fraction THEN
                state(p) = decayed
                atoms = atoms - 1
                CALL DrawPoint(px(p), py(p), decayed, pointsize)
        END IF
NEXT p

CALL PlotPoint
CALL PrintSummary

NEXT time

LOCATE 28, 1
PRINT "Press any key to continue."

END SUB

SUB TrialSetup

'calculate fraction of decays per second from halflife
fraction = 1 - 1 / (2 ^ (1 / halflife))

'select radioactive atoms with random positions
CALL InitAtoms

CLS

CALL PrintKeys
CALL DrawAtomBox
CALL DrawGraphBox

'initialise graph variables
atoms = natoms
time = 0

CALL PlotPoint
CALL PrintSummary

END SUB