PROGRAM sandpile
! simple one-dimensional sandpile model
DIM height(51),move(100),N(0:51)
CALL initial(height(),N(),L,sand$,#1,#2)
LET grain = 0
DO
   LET height(1) = height(1) + 1  ! add grain to first site
   WINDOW #1
   BOX SHOW sand$ at 1,height(1) - 1
   LET topple = 0                 ! number of grains that topple
   LET grain = grain + 1          ! number of grains added
   DO
      CALL check(height(),L,move(),unstable)
      IF unstable = 1 then
         CALL slide(height(),L,move(),sand$,#1)
         LET topple = topple + 1
      END IF
   LOOP until unstable = 0
   LET N(topple) = N(topple) + 1
   CALL show_distribution(N(),L,grain,#2)
LOOP until key input
END

SUB initial(height(),N(),L,sand$,#1,#2)
    OPEN #1: screen 0,0.7,0,1
    INPUT prompt "lattice size L = ": L     ! suggest L = 20
    FOR i = 1 to L
        LET height(i) = 0
    NEXT i
    LET height(L+1) = 0           ! boundary condition
    SET WINDOW -1,L+1,-2,40
    SET COLOR "blue"
    BOX AREA 1.1,1.9,0.1,0.9
    BOX KEEP 1,2,0,1 in sand$
    CLEAR
    SET COLOR "black"
    PLOT LINES: 1,-0.05;L+1,-0.05
    ! initialize array
    FOR topple = 0 to 20
        LET N(topple) = 0
    NEXT topple
    OPEN #2: screen 0.72,1,0.4,1
    SET WINDOW -0.1,L+1,-0.05,1.01
    PLOT LINES: 0,-0.01;L+0.2,-0.01
    PLOT LINES: -0.05,-0.01;-0.05,1
END SUB

SUB check(height(),L,move(),unstable)
    LET unstable = 0
    FOR i = 1 to L
        IF height(i) - height(i+1) > 1 then
           LET move(i) = 1
           LET unstable = 1
        ELSE
           LET move(i) = 0
        END IF
    NEXT i
END SUB

SUB slide(height(),L,move(),sand$,#1)
    WINDOW #1
    FOR i = 1 to L
        IF move(i) = 1 then
           LET height(i) = height(i) - 1    ! sand topples
           BOX CLEAR i,i + 1,height(i),height(i) + 1
           IF i < L then
              ! next site receives grain
              LET height(i+1) = height(i+1) + 1
              BOX SHOW sand$ at i+1,height(i+1) - 1
           END IF
        END IF
    NEXT i
END SUB

SUB show_distribution(N(),L,grain,#2)
    WINDOW #2
    ! erase previous graph
    SET COLOR "white"
    BOX AREA 0,L+0.2,0,1
    SET COLOR "black"
    FOR topple = 0 to L
        IF N(topple) > 0 then
           BOX AREA topple,topple+0.2,0,N(topple)/grain
        END IF
    NEXT topple
END SUB