PROGRAM invasion
! generate invasion percolation cluster
DIM site(0 to 200,0 to 100),perx(5000),pery(5000)
LIBRARY "csgraphics"
CALL initial(Lx,Ly,water$,#1)
CALL assign(site(,),perx(),pery(),nper,Lx,Ly,water$,#1)
CALL invade(site(,),perx(),pery(),nper,Lx,Ly,water$,#1)
CALL average(site(,),Lx,Ly,#1)
END

SUB initial(Lx,Ly,water$,#1)
    RANDOMIZE
    INPUT prompt "length in y direction = ": Ly
    LET Lx = 2*Ly
    OPEN #1: screen 0.01,0.65,0.2,0.99
    CALL compute_aspect_ratio(Lx,xwin,ywin)
    SET WINDOW 0,xwin,0,ywin
    SET COLOR "blue"
    BOX AREA 0,1,0,1
    FLOOD 0.5,0.5
    BOX KEEP 0,1,0,1 in water$
    CLEAR
    SET COLOR "black"
    BOX LINES 0.5,Lx+0.5,0.5,Ly+0.5
    FOR y = 1 to Ly
        FOR x = 1 to Lx
            PLOT POINTS: x,y
        NEXT x
    NEXT y
END SUB

SUB assign(site(,),perx(),pery(),nper,Lx,Ly,water$,#1)
    FOR y = 1 to Ly
        LET site(1,y) = 1         ! occupy first column
        BOX SHOW water$ at 0.5,y-0.5
    NEXT y
    ! assign random numbers to remaining sites
    FOR y = 1 to Ly
        FOR x = 2 to Lx
            LET site(x,y) = rnd
        NEXT x
    NEXT y
    ! sites in second column are initial perimeter sites
    ! site(x,y) greater than 2 if perimeter site
    LET x = 2
    LET nper = 0
    FOR y = 1 to Ly
        LET site(x,y) = 2 + site(x,y)
        LET nper = nper + 1       ! number of perimeter sites
        ! sort perimeter sites
        ! order perimeter list
        CALL insert(site(,),perx(),pery(),nper,x,y)
    NEXT y
END SUB

SUB insert(site(,),perx(),pery(),nper,x,y)
    ! call linear or binary sort
    CALL binary_search(site(,),perx(),pery(),nper,x,y,ninsert)
    ! move sites with smaller random numbers to next higher array index
    FOR ilist = nper to ninsert + 1 step - 1
        LET perx(ilist) = perx(ilist-1)
        LET pery(ilist) = pery(ilist-1)
    NEXT ilist
    LET perx(ninsert) = x         ! new site inserted in list
    LET pery(ninsert) = y
END SUB

SUB binary_search(site(,),perx(),pery(),nper,x,y,ninsert)
    ! divide list in half and determine half in which random # belongs
    ! continue this division until position of number is determined
    LET nfirst = 1                ! beginning of list
    LET nlast = nper - 1          ! end of list
    IF nlast < 1 then LET nlast = 1
    LET nmid = int((nfirst + nlast)/2)      ! middle of list
    ! determine which half of list new number is located
    DO
       IF nlast - nfirst <= 1 then     ! exact position equal to nlast
          LET ninsert = nlast
          EXIT SUB
       END IF
       LET xmid = perx(nmid)
       LET ymid = pery(nmid)
       IF site(x,y) > site(xmid,ymid) then
          LET nlast = nmid        ! search upper half
       ELSE
          LET nfirst = nmid       ! search lower half
       END IF
       LET nmid = int((nfirst + nlast)/2)
    LOOP
END SUB

SUB linear_search(site(,),perx(),pery(),nper,x,y,ninsert)
    IF nper = 1 then
       LET ninsert = 1
    ELSE
       FOR iper = 1 to nper - 1
           LET xperim = perx(iper)
           LET yperim = pery(iper)
           IF site(x,y) > site(xperim,yperim) then
              LET ninsert = iper  ! insert new site
              EXIT SUB
           END IF
       NEXT iper
    END IF
    LET ninsert = nper
END SUB

SUB invade(site(,),perx(),pery(),nper,Lx,Ly,water$,#1)
    ! nx and ny are components of vectors pointing to nearest neighbors
    DIM nx(4),ny(4)
    DATA 1,0,-1,0,0,1,0,-1
    FOR i = 1 to 4
        READ nx(i),ny(i)
    NEXT i
    DO
       LET x = perx(nper)
       LET y = pery(nper)
       LET nper = nper - 1
       ! mark site occupied and no longer perimeter site
       LET site(x,y) = site(x,y) - 1    
       BOX SHOW water$ at x-0.5,y-0.5
       FOR i = 1 to 4             ! find new perimeter sites          
           LET xnew = x + nx(i)
           LET ynew = y + ny(i)
           IF ynew > Ly then      ! periodic boundary conditions in y 
              LET ynew = 1
           ELSE IF ynew < 1 then
              LET ynew = Ly
           END IF
           IF site(xnew,ynew) < 1 then      ! new perimeter site
              LET site(xnew,ynew) = site(xnew,ynew) + 2
              LET nper = nper + 1
              CALL insert(site(,),perx(),pery(),nper,xnew,ynew)
           END IF
       NEXT i
    LOOP until x >= Lx            ! until cluster reaches right boundary
END SUB

SUB average(site(,),Lx,Ly,#1)
    ! compute probability density P(r)
    DIM P(0 to 20),nr(0 to 20)
    LET Lmin = Lx/3
    LET Lmax = 2*Lmin
    LET n = (Lmax - Lmin + 1)*Ly  ! # sites in middle half of lattice
    LET dr = 0.05
    LET nbin = 1/dr
    FOR x = Lmin to Lmax
        FOR y = 1 to Ly
            LET ibin = nbin*(mod(site(x,y),1))
            LET nr(ibin) = nr(ibin) + 1
            IF site(x,y) >= 1 and site(x,y) < 2 then
               LET occupied = occupied + 1  ! total # of occupied sites
               LET P(ibin) = P(ibin) + 1
            END IF
        NEXT y
    NEXT x
    WINDOW #1
    PRINT "# occupied sites ="; occupied
    OPEN #2: screen 0.66,1.0,0.01,0.99
    PRINT " r","  P(r)"
    PRINT
    FOR ibin = 0 to nbin
        LET r = dr*ibin
        IF nr(ibin) > 0 then PRINT r,P(ibin)/nr(ibin)
    NEXT ibin
    PRINT
END SUB