PROGRAM perc_cluster
*     cluster generated by Hammersley, Leath, and Alexandrowicz algorithm
*    
      IMPLICIT NONE
      REAL*8 p
      INTEGER IR_, L, N, xs(10000), ys(10000), status(0:2)
      DATA xs,ys/20000*0/
      CALL GWopen(IR_, 0)
      CALL initial(p,L,status)
      CALL grow(p,L,N,xs,ys,status)
      CALL mass_dist(L,N,xs,ys)
      CALL GWquit(IR_)
      END
C
      SUBROUTINE initial(p,L,status)
      IMPLICIT NONE
      REAL*8 p
      REAL xwin, ywin, sz
      REAL dmy, rnd
      INTEGER IR_, L, x, y, status(0:2)
      dmy = rnd(-1)
      WRITE(*,'(A,$)') 'value of L (odd) = '
      READ(*,*) L
      WRITE(*,'(A,$)') 'site occupation probability = '
      READ(*,*) p
      CALL compute_aspect_ratio(REAL(L),xwin,ywin)
      CALL GWindow(IR_, 0.0, 0.0, xwin, ywin)
      sz = ywin/30
      CALL GWsettxt(IR_, sz, 0.0, 5, -1, -1, '')
      CALL GWsetpen(IR_, 0, -1, -1, -1)
      CALL GWrect(IR_, 0.5, 0.5, L+0.5, L+0.5)
      DO y = 1, L
         DO x = 1, L
            CALL GWsetpxl(IR_, REAL(x), REAL(y), 0)
         END DO
      END DO
      CALL GWsetpen(IR_, 13, -1, -1, 4)
      CALL GWcmbmrk(status(1), 2, 6, 1.0, 13, -1, 4)  ! occupied site
      CALL GWputcmb(IR_, status(1), 0.5+sz/2, L+0.5+sz, sz, sz, 0)
      CALL GWputtxt(IR_, 0.5+sz*1.5, L+0.5+sz, 'occupied site')
      CALL GWsetpen(IR_, 16, -1, -1, 4)
      CALL GWcmbmrk(status(2), 3, 6, 1.0, 16, -1, 4)  ! perimeter site
      CALL GWputcmb(IR_, status(2), 
     +              L/3.0+0.5+sz/2, L+0.5+sz, sz, sz, 0)
      CALL GWputtxt(IR_, 
     +              L/3.0+0.5+sz*1.5, L+0.5+sz, 'perimeter site')
      CALL GWsetpen(IR_, 0, -1, -1, 4)
      CALL GWcmbmrk(status(0), 4, 6, 1.0, 0, -1, 4)   ! tested, unoccupied site
      CALL GWputcmb(IR_, status(0), 
     +              L/3.0*2+0.5+sz/2, L+0.5+sz, sz, sz, 0)
      CALL GWputtxt(IR_, 
     +              L/3.0*2+0.5+sz*1.5, L+0.5+sz, 'tested site')
      END
C
      SUBROUTINE grow(p,L,N,xs,ys,status)
*     generate single percolation cluster          
*     set up boundary sites
      IMPLICIT NONE
      REAL*8 p
      REAL rnd
      INTEGER nx(4), ny(4), site(0:131,0:131), L, xs(10000),
     +       ys(10000), status(0:2)
      INTEGER IR_, i, iper, nn, nper, perx(25000), pery(25000), x,
     +       xnew, xseed, y, ynew, yseed, N
      LOGICAL Ltmp1_
      INTEGER DP_, DATA_(8)
      DATA DP_/1/,DATA_/1,0,-1,0,0,1,0,-1/
      DO i = 1, L
         site(0,i) = -1
         site(L+1,i) = -1
         site(i,L+1) = -1
         site(i,0) = -1
      END DO
*     seed at center of lattice
      xseed = L/2 + 1
      yseed = xseed
      site(xseed,yseed) = 1                 ! seed site
      xs(1) = xseed
      ys(1) = yseed
      N = 1                                 ! number of sites in the cluster
      CALL GWputcmb(IR_, status(1), REAL(xseed), REAL(yseed), 
     +              1.0, 1.0, 0)
      DO i = 1, 4
*        nx,ny direction vectors for new perimeter sites
         nx(i) = DATA_(DP_)
         ny(i) = DATA_(DP_ + 1)
         DP_ = DP_ + 2
*        perx,pery, positions of perimeter sites of seed
         perx(i) = xseed + nx(i)
         pery(i) = yseed + ny(i)
*        perimeter sites labeled by 2
         site(perx(i),pery(i)) = 2          ! site placed on perimeter list
         CALL GWputcmb(IR_, status(2), REAL(perx(i)), REAL(pery(i)), 
     +                 1.0, 1.0, 0)
      END DO
      nper = 4                              ! initial number of perimeter sites
      Ltmp1_ = .TRUE.
      DO WHILE(Ltmp1_)
*        randomly choose perimeter site
         iper = int(rnd(0)*nper) + 1
         x = perx(iper)                     ! coordinate of a perimeter site
         y = pery(iper)
*        relabel remaining perimeter sites so that
*        last perimeter site in array replaces newly chosen site
         perx(iper) = perx(nper)
         pery(iper) = pery(nper)
         nper = nper - 1
         IF(rnd(0) .LT. p) THEN
            site(x,y) = 1
            N = N + 1
            xs(N) = x                       ! save position of occupied site
            ys(N) = y
            CALL GWputcmb(IR_, status(1), REAL(x), REAL(y), 
     +                    1.0, 1.0, 0)
            DO nn = 1, 4
               xnew = x + nx(nn)
               ynew = y + ny(nn)
               IF(site(xnew,ynew) .EQ. 0) THEN
                  nper = nper + 1
                  perx(nper) = xnew
                  pery(nper) = ynew
*                 place site on perimeter list
                  site(xnew,ynew) = 2
                  CALL GWputcmb(IR_, status(2), REAL(xnew), REAL(ynew), 
     +                          1.0, 1.0, 0)
               END IF
            END DO
         ELSE
            site(x,y) = -1
            CALL GWputcmb(IR_, status(0), REAL(x), REAL(y), 
     +                    1.0, 1.0, 0)
         END IF
         Ltmp1_ = .NOT.(nper .LT. 1)
      END DO
      CALL GWsetbrs(IR_, 0, 0, -1)
      CALL GWrect(IR_, 0.5, 0.5, L+0.5, L+0.5)  ! redraw box
      END
C
      SUBROUTINE mass_dist(L,N,xs,ys)
      IMPLICIT NONE
      REAL*8 dx, dy, xcm, ycm
      INTEGER i, r, rprint, N, L, 
     +        masstotal, mass(10000), xs(10000), ys(10000)
      DATA xcm, ycm/2*0.0D0/
      DATA masstotal, mass/10001*0/
      DO i = 1, N
         xcm = xcm + xs(i)                  ! compute center of mass
         ycm = ycm + ys(i)
      END DO
      xcm = xcm/N
      ycm = ycm/N
      DO i = 1, N
         dx = xs(i) - xcm
         dy = ys(i) - ycm
         r = int(sqrt(dx*dx + dy*dy))
*        distance from center of mass
*        mass(r) = number of sites at distance r from center of mass
         IF(r .GT. 1) mass(r) = mass(r) + 1
      END DO
      rprint = 2
      WRITE(*,'(2A9,2A13)') 'r', 'm', 'ln(r)', 'ln(m)'
      DO r = 2, L/2
         masstotal = masstotal + mass(r)
         IF(r .EQ. rprint) THEN
            WRITE(*,'(2I9,2F13.6)') 
     +         r, masstotal, log(REAL(r)), log(REAL(masstotal))
            rprint = 2*rprint               ! use logarithmic scale for r
         END IF
      END DO
      END