      PROGRAM continuum
*     generate continuum percolation configurations
      IMPLICIT NONE
      INTEGER IR_, L
      CALL GWopen(IR_, 0)
      CALL initial(L)
      CALL disks(L)                         ! place disks at random inside box
      CALL GWquit(IR_)
      END
C
      SUBROUTINE initial(L)
      IMPLICIT NONE
      CHARACTER*80 Stmp1_
      REAL xwin, ywin
      REAL dmy, rnd
      INTEGER IR_, L
      dmy = rnd(-1)
      CALL GWinput(IR_, 'linear dimension of lattice = ', Stmp1_)
      READ(Stmp1_,*) L
      CALL compute_aspect_ratio(REAL(L),xwin,ywin)
      CALL GWindow(IR_, -0.1*xwin, -0.05*ywin, xwin, 1.05*ywin)
      CALL GWcolor(IR_, 0, 2)
      CALL GWcolor(IR_, 19, 7)
      CALL GWsetpen(IR_, 13, -1, -1, 4)     ! erase disk
      CALL GWsetbrs(IR_, 13, 0, -1)
      CALL GWrect(IR_, 0.0, 0.0, REAL(L), REAL(L))
      CALL GWanchor(IR_, 1)
      END
C
      SUBROUTINE disks(L)
      IMPLICIT NONE
      CHARACTER*80 Stmp1_
      REAL*8 area, density, phi, pi, r, x(10000), y(10000)
      PARAMETER(pi = 3.141592653589793D0)
      REAL xs, ys, rnd
      INTEGER IR_, L, n, number, i
      LOGICAL Ltmp1_
      DATA n, number/2*0/
      r = 1                                 ! radius of disks
      area = L*L
      CALL GWsetpen(IR_, 15, -1, -1, 4)
      CALL GWsetmrk(IR_, 6, 1.0, 15, -1, 4)
      DO WHILE(n .GE. 0)
         CALL GWsetogn(IR_, 0,-10)
         CALL GWinput(IR_, 'number of additional disks = ', Stmp1_)
         IF(IR_ .EQ. 0 ) RETURN
         READ(Stmp1_, *) n
         density = (n + number)/area
         phi = 1 - exp(-density*pi*r*r)
         WRITE(Stmp1_, *) 'number of disks =', n + number
         CALL GWputtxt(IR_, 0.0, REAL(L*25.0/24), Stmp1_)
         WRITE(Stmp1_, *) 'density of disks =', density
         CALL GWputtxt(IR_, 0.0, REAL(L*27.0/24), Stmp1_)
         WRITE(Stmp1_, *) 'phi =', phi
         CALL GWputtxt(IR_, 0.0, REAL(L*26.0/24), Stmp1_)
         DO i = 1, number
            CALL GWputmrk(IR_, REAL(x(i)), REAL(y(i)))
         END DO
         DO i = number+1, n+number
            x(i) = L*rnd(0)                 ! center of disk inside box
            y(i) = L*rnd(0)
            CALL GWputmrk(IR_, REAL(x(i)), REAL(y(i)))
         END DO
         number = number + n
         CALL GWflush(IR_, -10)
         CALL GWsetogn(IR_, 0, 2)
         CALL GWsetbrs(IR_, 16, 1, -1)
         Ltmp1_ = .TRUE.
         DO WHILE(Ltmp1_)
            CALL GWcappnt(IR_, xs, ys, '')
*           click on occupied site to see cluster
*           click mouse to right of lattice to increase p by delta
            IF(xs .LT. L .AND. ys .LE. L) THEN
               CALL GWflood(IR_, REAL(xs), REAL(ys))
            END IF
            Ltmp1_ = xs .LE. L
         END DO
         CALL GWerase(IR_, 2, -1)
      END DO
      END
