PROGRAM perc_cluster
! cluster generated by Hammersley, Leath, and Alexandrowicz algorithm
DIM xs(10000),ys(10000),status$(-1:2)
LIBRARY "csgraphics"
CALL initial(p,L,status$())
CALL initialize_arrays(xs(),ys())
CALL grow(p,L,N,xs(),ys(),status$())
CALL mass_dist(L,N,xs(),ys())
END

SUB initial(p,L,status$())
    RANDOMIZE
    INPUT prompt "value of L (odd) = ": L
    INPUT prompt "site occupation probability = ": p
    CALL compute_aspect_ratio(L,xwin,ywin)
    SET WINDOW 0,xwin,0,ywin
    SET COLOR "red"
    BOX CIRCLE 0,1,0,1
    FLOOD 0.5,0.5
    BOX KEEP 0,1,0,1 in status$(1)     ! occupied site
    CLEAR
    SET COLOR "blue"
    BOX CIRCLE 0,1,0,1
    FLOOD 0.5,0.5
    BOX KEEP 0,1,0,1 in status$(2)     ! perimeter site
    CLEAR
    SET COLOR "black"
    BOX CIRCLE 0,1,0,1
    FLOOD 0.5,0.5
    BOX KEEP 0,1,0,1 in status$(-1)    ! tested, unoccupied site
    CLEAR
    BOX LINES 0.5,L+0.5,0.5,L+0.5
    FOR y = 1 to L
        FOR x = 1 to L
            PLOT POINTS: x,y
        NEXT x
    NEXT y
    BOX SHOW status$(1) at 1,L+2
    PLOT TEXT, AT 1,L+2: "  occupied site"
    BOX SHOW status$(2) at 30,L+2
    PLOT TEXT, AT 30,L+2: "  perimeter site"
    BOX SHOW status$(-1) at 60,L+2
    PLOT TEXT, AT 60,L+2: "  tested site"
END SUB

SUB initialize_arrays(xs(),ys())
    MAT xs = 0
    MAT ys = 0
END SUB

SUB grow(p,L,N,xs(),ys(),status$())
    ! generate single percolation cluster          
    DIM perx(25000),pery(25000),site(0:131,0:131)
    DIM nx(4),ny(4)          ! set up direction vectors for lattice
    DATA 1,0,-1,0,0,1,0,-1
    ! set up boundary sites
    FOR i = 1 to L
        LET site(0,i) = -1
        LET site(L+1,i) = -1
        LET site(i,L+1) = -1
        LET site(i,0) = -1
    NEXT i
    ! seed at center of lattice
    LET xseed = int(L/2) + 1
    LET yseed = xseed
    LET site(xseed,yseed) = 1     ! seed site
    LET xs(1) = xseed
    LET ys(1) = yseed
    LET N = 1                ! number of sites in the cluster
    BOX SHOW status$(1) at xseed-0.5,yseed-0.5
    FOR i = 1 to 4
        ! nx,ny direction vectors for new perimeter sites
        READ nx(i),ny(i)
        ! perx,pery, positions of perimeter sites of seed
        LET perx(i) = xseed + nx(i)
        LET pery(i) = yseed + ny(i)
        ! perimeter sites labeled by 2
        LET site(perx(i),pery(i)) = 2  ! site placed on perimeter list
        BOX SHOW status$(2) at perx(i)-0.5,pery(i)-0.5
    NEXT i
    LET nper = 4             ! initial number of perimeter sites
    DO
       ! randomly choose perimeter site
       LET iper = int(rnd*nper) + 1
       LET x = perx(iper)    ! coordinate of a perimeter site
       LET y = pery(iper)
       ! relabel remaining perimeter sites so that
       ! last perimeter site in array replaces newly chosen site
       LET perx(iper) = perx(nper)
       LET pery(iper) = pery(nper)
       LET nper = nper - 1
       IF rnd < p then       ! site occupied
          LET site(x,y) = 1
          LET N = N + 1
          LET xs(N) = x      ! save position of occupied site
          LET ys(N) = y
          BOX SHOW status$(1) at x-0.5,y-0.5
          FOR nn = 1 to 4    ! find new perimeter sites
              LET xnew = x + nx(nn)
              LET ynew = y + ny(nn)
              IF site(xnew,ynew) = 0 then
                 LET nper = nper + 1
                 LET perx(nper) = xnew
                 LET pery(nper) = ynew
                 ! place site on perimeter list
                 LET site(xnew,ynew) = 2
                 BOX SHOW status$(2) at xnew-0.5,ynew-0.5
              END IF
          NEXT nn
       ELSE                  ! rnd >= p so site is not occupied        
          LET site(x,y) = -1
          BOX SHOW status$(-1) at x-0.5,y-0.5
       END IF
    LOOP until nper < 1      ! all perimeter sites tested
    BOX LINES 0.5,L+0.5,0.5,L+0.5      ! redraw box
END SUB

SUB mass_dist(L,N,xs(),ys())
    DIM mass(10000)
    PRINT "press any key or click mouse to see data"
    DO
       GET MOUSE xm,ym,s
    LOOP until key input or s <> 0
    FOR i = 1 to N
        LET xcm = xcm + xs(i)     ! compute center of mass
        LET ycm = ycm + ys(i)
    NEXT i
    LET xcm = xcm/N
    LET ycm = ycm/N
    FOR i = 1 to N
        LET dx = xs(i) - xcm
        LET dy = ys(i) - ycm
        LET r = int(sqr(dx*dx + dy*dy))
        ! distance from center of mass
        ! mass(r) = number of sites at distance r from center of mass
        IF r > 1 then LET mass(r) = mass(r) + 1
    NEXT i
    LET rprint = 2
    CLEAR
    PRINT " r "," m "," ln(r) "," ln(m) "
    FOR r = 2 to L/2
        LET masstotal = masstotal + mass(r)
        IF r = rprint then
           PRINT r, masstotal,log(r),log(masstotal)
           LET rprint = 2*rprint  ! use logarithmic scale for r
        END IF
    NEXT r
END SUB