program rwalk04
c-----------------------------------------------------------------
c Random walk in 2 dimensions (self-avoiding 2D random walk)
c 2D Polymer simulation
c-----------------------------------------------------------------
c input from a keyboard:
c    number of trials
c output screen:
c    initial data, average polymer size
c output file:
c    polymer size for each trial (may be used later for histograms)
c Attention!!!
c the code call subroutine 'gettim' to get current time
c that name is not unique, and may be different for different compilers
c You may replace the seed number from current time by some seed number
c-----------------------------------------------------------------
c written by: Alex Godunov
c last revision: November 2006
c-----------------------------------------------------------------
      implicit none
      integer*4 iu, it, is, itests, iway, x, y
      integer*4 size, n, out, ntot, xold, yold, i, j
      integer*4 ihour, imin, isec, imsec
      parameter (size=1000)
      integer*4 a(1000,1000)
      real*4 rand, d, dav, nav, c

      write (*,*) ' 2D Polymer simulation'
      write (*,*) ' enter numbers of tests'
      read  (*,*) itests

      call gettim(ihour,imin,isec,imsec)
      iu = imsec*isec*imin+ihour
      call srand(iu)
      dav=0.0
      nav=0.0
      ntot=0
      open (unit=7,file='polymer.dat')

      do it=1,itests

c initialization before new trial
        x=size/2
        y=size/2
ccc        write(7,103) x-size/2,y-size/2
        out=0
        n=0

        do i=1,size
          do j=1,size
            a(i,j) = 0
          end do
        end do
c---
        do while (out.eq.0)
          xold=x
          yold=y
          c = rand()
            if(c.le.0.25) x = x+1
            if(c.gt.0.25.and.c.le.0.50) x=x-1
            if(c.gt.0.50.and.c.le.0.75) y=y+1
            if(c.gt.0.75) y=y-1
c is the (x,y) site vacant them accept the last step
          if(a(x,y).eq.0) then
            a(x,y)=1
            n = n+1
c output to plot a polymer
c            write(7,103) x-size/2,y-size/2
            else
              x=xold
              y=yold
          end if


c*** exit conditions
c condition 1: polymer is too large for a(x,y) array
c if x or y on array boundaries for a(x,y) then stop simultaion
          if(x.eq.1.or.x.eq.size.or.y.eq.1.or.y.eq.size) then
            write(*,102) x,y
            exit
          end if
c condition 2: no more available steps around new x,y where to go
          if(a(x+1,y)*a(x,y+1)*a(x-1,y)*a(x,y-1).eq.1) out=1

        end do

        d = sqrt((float(x-size/2))**2+(float(y-size/2))**2)
        dav = dav + d
        ntot = ntot + n
	  if((it/100)*100.eq.it) write(*,103) it, n
	  write(7,103) it, n
      end do

      dav = dav/float(itests)
      nav = float(ntot)/float(itests)
      write(*,100) itests, nav, dav

100   format(/,' itests: ',i7,/,' size: ',1pe14.4,/,' Ree:  ',e14.4)
101   format(2i8)
102   format('x or y exceeds a(x,y) size',2i10)
103   format(2i8)

      stop
      end