PROGRAM dla
DIM site(-100 to 100,-100 to 100)
LIBRARY "csgraphics"
CALL initial(site(,),L)
CALL grow_cluster(site(,),L)
END

SUB initial(site(,),L)
    RANDOMIZE
    INPUT prompt "L = ": L
    CALL compute_aspect_ratio(L+2,xwin,ywin)
    SET WINDOW -xwin,xwin,-ywin,ywin
    BOX LINES -L-0.5,L+0.5,-L-0.5,L+0.5
    MAT site = 0             ! initialize sites
    FOR y = -L to L
        FOR x = -L to L
            PLOT POINTS: x,y
        NEXT x
    NEXT y
    SET COLOR "red"
END SUB

SUB grow_cluster(site(,),L)
    LET R0 = 3               ! start walker at distance R0 from origin
    LET site(0,0) = 1        ! seed site
    LET N = 1                ! number of particles in cluster
    BOX AREA -0.5,0.5,-0.5,0.5
    DO
       ! find random initial position of new walker
       LET theta = 2*pi*rnd
       LET x = int(R0*cos(theta))
       LET y = int(R0*sin(theta))
       CALL walk(site(,),L,x,y,R0,N)   ! random walk 
    LOOP until key input
END SUB

SUB walk(site(,),L,x,y,R0,N)
    ! walk until on a perimeter site of cluster  
    ! or walker strays too far from cluster
    DO
       LET onperimeter$ = "no"
       LET r = sqr(x*x + y*y)
       IF r < R0 + 1 then
          ! test if walker on perimeter
          CALL test(site(,),L,x,y,r,R0,N,onperimeter$)
       END IF
       LET step = int(r - R0) - 1      ! big step 
       IF step < 1 then LET step = 1
       IF onperimeter$ = "no" then
          LET random = rnd
          IF random < 0.25 then
             LET x = x + step
          ELSE IF random < 0.5 then
             LET x = x - step
          ELSE IF random < 0.75 then
             LET y = y + step
          ELSE
             LET y = y - step
          END IF
       END IF
    LOOP until onperimeter$ = "yes" or r > 2*R0
END SUB

SUB test(site(,),L,x,y,r,R0,N,onperimeter$)
    LET sum = site(x,y+1) + site(x,y-1) + site(x+1,y) + site(x-1,y)
    IF sum > 0 then          ! walker on perimeter site 
       LET site(x,y) = 1
       LET onperimeter$ = "yes"
       IF abs(x) <= L and abs(y) <= L then
          BOX AREA x-0.5,x+0.5,y-0.5,y+0.5
          LET N = N + 1
          SET CURSOR 2,1
          PRINT using "N = ####": N
       ELSE
          STOP               ! cluster outside box
       END IF
       IF r >= R0 then LET R0 = int(r+2)
    END IF
END SUB