PROGRAM cca
PUBLIC site(50,50),x(1000),y(1000),L,N
LIBRARY "csgraphics"
DIM next_particle(1000),first_particle(1000)
CALL initial(first_particle(),next_particle(),ncl)
DO
   CALL move(next_particle(),first_particle(),ncl)
LOOP until ncl = 1
END

SUB initial(first_particle(),next_particle(),ncl)
    PUBLIC nx(4),ny(4)
    DECLARE PUBLIC site(,),x(),y()
    DECLARE PUBLIC L,N
    RANDOMIZE
    DATA 1,0,-1,0,0,1,0,-1
    FOR nn = 1 to 4
        READ nx(nn),ny(nn)
    NEXT nn
    INPUT prompt "L = ": L
    INPUT prompt "N = ": N
    CALL compute_aspect_ratio(L+1,xwin,ywin)
    SET WINDOW 0,xwin,0,ywin
    CLEAR
    BOX LINES 1,L+1,1,L+1
    SET COLOR "blue"
    LET ncl = 0                   ! number of clusters
    FOR i = 1 to N
        DO
           LET x(i) = int(L*rnd) + 1
           LET y(i) = int(L*rnd) + 1
        LOOP until site(x(i),y(i)) = 0
        LET ncl = ncl + 1
        LET site(x(i),y(i)) = ncl
        BOX AREA x(i),x(i)+1,y(i),y(i)+1
        LET first_particle(ncl) = i
        LET next_particle(i) = i
        LET xi = x(i)
        LET yi = y(i)
        CALL neighbor(xi,yi,next_particle(),first_particle(),ncl)
    NEXT i
END SUB

SUB neighbor(xi,yi,next_particle(),first_particle(),ncl)
    DECLARE PUBLIC site(,),x(),y()
    DECLARE PUBLIC L,nx(),ny()
    DECLARE DEF pbc
    FOR nn = 1 to 4
        LET px = pbc(xi + nx(nn),L)
        LET py = pbc(yi + ny(nn),L)
        LET perim = site(px,py)
        LET part = site(xi,yi)
        IF perim <> 0 and perim <> part then
           CALL merge(perim,part,first_particle(),next_particle(),ncl)
        END IF
    NEXT nn
END SUB

SUB merge(c1,c2,first_particle(),next_particle(),ncl)
    DECLARE PUBLIC site(,), x(), y()
    LET p1 = first_particle(c1)
    LET p2 = first_particle(c2)
    LET p1next = next_particle(p1)
    LET p2next = next_particle(p2)
    LET next_particle(p1) = p2next
    LET next_particle(p2) = p1next
    DO
       LET site(x(p2next),y(p2next)) = c1
       LET p2next = next_particle(p2next)
    LOOP until p2next = p1next
    LET plast = first_particle(ncl)
    IF c2 <> ncl then
       LET p = plast
       DO
          LET site(x(p),y(p)) = c2
          LET p = next_particle(p)
       LOOP until p = plast
       LET first_particle(c2) = plast
    END IF
    LET ncl = ncl - 1
END SUB

SUB move(next_particle(),first_particle(),ncl)
    DECLARE PUBLIC site(,),x(),y()
    DECLARE PUBLIC L,nx(),ny()
    DECLARE DEF pbc
    LET c = int(ncl*rnd) + 1
    LET direction = int(4*rnd) + 1
    LET p1 = first_particle(c)
    LET i = p1
    LET dx = nx(direction)
    LET dy = ny(direction)
    DO
       LET site(x(i),y(i)) = 0
       SET COLOR "white"
       BOX AREA x(i),x(i)+1,y(i),y(i)+1
       LET x(i) = pbc(x(i) + dx,L)
       LET y(i) = pbc(y(i) + dy,L)
       LET i = next_particle(i)
    LOOP until i = p1
    DO
       SET COLOR "blue"
       BOX AREA x(i),x(i) + 1, y(i),y(i) + 1
       LET site(x(i),y(i)) = c
       LET i = next_particle(i)
    LOOP until i = p1
    DO
       LET xi = x(i)
       LET yi = y(i)
       CALL neighbor(xi,yi,next_particle(),first_particle(),ncl)
       LET i = next_particle(i)
    LOOP until i = p1
    SET COLOR "black"
    BOX LINES 1,L + 1,1,L + 1
END SUB

FUNCTION pbc(s,L)
    IF s > L then
       LET pbc = 1
    ELSE IF s < 1 then
       LET pbc = L
    ELSE
       LET pbc = s
    END IF
END DEF