PROGRAM continuum
! generate continuum percolation configurations
LIBRARY "csgraphics"
CALL initial(L,occup$)
CALL disks(L,occup$)              ! place disks at random inside box
END

SUB initial(L,occup$)
    RANDOMIZE
    INPUT prompt "linear dimension of lattice = ": L
    CALL compute_aspect_ratio(L,xwin,ywin)
    CLEAR
    SET WINDOW -0.1*xwin,xwin,-0.05*ywin,1.05*ywin
    SET BACKGROUND COLOR "black"
    SET COLOR "red"
    BOX AREA 0,L,0,L
    SET COLOR "yellow"
    BOX CIRCLE 2,3,2,3
    FLOOD 2.5,2.5
    BOX KEEP 2,3,2,3 in occup$
    SET COLOR "red"               ! erase disk
    BOX CIRCLE 2,3,2,3
    FLOOD 2.5,2.5
END SUB

SUB disks(L,occup$)
    DIM x(10000),y(10000)
    LET r = 1                     ! radius of disks
    LET area = L*L
    DO while n >= 0
       SET CURSOR 1,1
       SET COLOR "white"
       LET density = number/area
       LET phi = 1 - exp(-density*pi*r*r)
       PRINT using "density of disks = #.###": density
       PRINT using "phi = #.###": phi
       SET CURSOR 3,29
       PRINT "     ";
       SET CURSOR 3,1
       INPUT prompt "number of additional disks = ": n
       SET COLOR "red"
       BOX AREA 0,L,0,L
       FOR i = 1 to number
           BOX SHOW occup$ at x(i)-0.5,y(i)-0.5 using "or"
       NEXT i
       FOR i = number+1 to n+number
           LET x(i) = L*rnd       ! center of disk inside box
           LET y(i) = L*rnd
           BOX SHOW occup$ at x(i)-0.5,y(i)-0.5 using "or"
       NEXT i
       LET number = number + n
       DO
          DO
             GET MOUSE xs,ys,state
          LOOP until state = 2
          ! click on occupied site to see cluster
          ! click mouse to right of lattice to increase p by delta
          IF xs < L and ys <= L then
             SET COLOR "blue"
             FLOOD xs,ys
          END IF
       LOOP until xs > L
    LOOP
END SUB