program awalk

c========1=2a==================================================
c
c
c        This program simulates the random walk on a disordered
c        lattice. The disordered lattice is generated by the
c        Leath algorithm. Hence, the resulting structure is a
c        percolation cluster. There is only one cluster in the
c        system and its properties can be controlled by specifying
c        the percolation probability p
c
c        there will be different versions of this program:
c
c        1. the random walker is completely unbiased
c
c        2. there is a waiting time distribution for the
c           walker
c           a)   dichotomic
c           b)   some function, which will be realized by
c                a markov chain
c
c        3. inclusion of local anisotroy, however, isotropic
c           overall
c
c
c        This version of the program does not contain graphic!
c        For the graphic demo see
c                  view_demo.f
c
c        The program does restart with previous results for the
c        random walker
c
c        the guilty parties for this program are
c
c        Dieter  W. Heermann
c        Joachim Koehler
c        Michael Kotelyanskii
c        Andrei  Gusev
c        Simone  Arizzi
c 
c        Zuerich  july 1990
c
c
c        this is (2)
c
c        random number generator is R250
c
c================================================================
c  
c
         parameter (gmax=50,lmax=50000,tmax=500)
         parameter (numwal=10000,maxtry=2000)
c
c        the parameters are as follows:

c        gmax       : this is the size of the array holding the
c                     random lattice
c
c        lmax       : size of the stack holding the surface
c                     sites
c
c        numwal     : number of random walkers which will be
c                     started on the random lattice
c
c        maxtry     : maximum number of tries to place a walker
c                     on the random lattice
c
c        tmax       : maximum nuber of steps the random walker
c                     can perform
c
c===============================================================
c
c--------declarations
c
         integer g(-gmax:gmax,-gmax:gmax)
         integer list(0:1,1:lmax,1:2)
c
         integer mz(1:250)
c
         integer nshell,shmax
         integer iseed,seed
         integer nsurf
         integer id,ptr
         integer trial,try,again
c
         integer x0,y0,x,y,k,const
         integer time
c
         real    r2(tmax),count(tmax)
         real    xt(1:2,1:tmax)
         real    rana(1:tmax)
c
         real    p
         real    mass,radius,rgyr
c
         double precision rand,ran

c--------read in or set the initial values

         write(*,*) 'what is the percolation probability        : '
         read(*,*)  p
         write(*,*) 'give the number of shell to be grown       : '
         read(*,*)  shmax
         write(*,*) 'seed for the random number generator       : '
         read(*,*)  iseed
         write(*,*) 'how times should we try                    : '
         read(*,*)  trial
         write(*,*) 'shall we read in previous results (0,1)    : '
         read(*,*) again


         if ( again .eq. 0 ) then
c
c-----------here we initialize 'stuff'

            do 5009 i=1,tmax
               r2(i)=0
               count(i)=0
 5009        continue
c
           else
             write(*,*) 'restarting ... '
             read(11,*) iseed
             do   6010  i=1,tmax
                read(11,*) r2(i),count(i)
 6010        continue
             write(*,*) 'finished reading r2 and count '
             iseed = iseed + 3
          end if
c
c
c--------initialize the random number generator
c
         write(*,*) 'initializing the random number generator ... '
         seed = iseed
         call inr250(mz,iseed)
         nr = tmax
         write(*,*) 'done '

         write(*,*) 'starting to set the percolation cluster ... '

         do  111  ntrial=1,trial
           
           call srand((seed + ntrial))
           write(*,*) 'this trial number : ',ntrial,' ... '

c--------initialize the entire lattice with zeros
c
         do 100  i=-gmax,gmax
           do  100  j=-gmax,gmax
             g(i,j) = 0
 100     continue

         rgyr   = 0.0
         radius = 0.0
         mass   = 1.0

c===============================================================
c
c        now we set up the percolation lattice using the
c        leath idea to create just the percolation cluster
c
c===============================================================

c        start out with the seed in the center of the lattice
c
         g(0,0) = 2
c
c        create the initial list of surface sites which will
c        potentialy belong to the percolation cluster
c
         list(0,1,1) = 0
         list(0,1,2) = 1
         list(0,2,1) = -1
         list(0,2,2) = 0
         list(0,3,1) = 0
         list(0,3,2) = -1
         list(0,4,1) = 1
         list(0,4,2) = 0
c
         id    = 0
         ptr   = 0
         nsurf = 4
c

         do  400  nshell=1,shmax


c        the meaning of the lattice settings are as follows:
c        virgin site     0
c        unoccupied      1
c        occupied        2

c        now with grow the cluster one layer
c
         do  500  i=1,nsurf
c
c          here we go
           modid = mod(id+1,2)
           ix = list(id,i,1)
           iy = list(id,i,2)
c
c          check if the surface site is already to close for comfort
           if ( iabs(ix) .eq. (gmax-1) )  goto 554
           if ( iabs(iy) .eq. (gmax-1) )  goto 554
c
c          ok, we can press on to occupy or not the site
           if ( g(ix,iy) .eq. 0 ) then
               ran = rand()
               if ( ran .lt. p ) then
                  g(ix,iy) = 2
                  mass     = mass + 1
                  rad      = ix * ix + iy * iy
                  rad      = sqrt( rad )
                  if ( rad .gt. radius ) radius = rad

c                 search for the new surface sites
                  if ( g(ix,iy+1) .eq. 0 ) then
                         ptr = ptr + 1
                         if ( ptr .gt. lmax ) goto 555
                         list(modid,ptr,1) = ix
                         list(modid,ptr,2) = iy+1
                      end if
                  if ( g(ix-1,iy) .eq. 0 ) then
                         ptr = ptr + 1
                         if ( ptr .gt. lmax ) goto 555
                         list(modid,ptr,1) = ix-1
                         list(modid,ptr,2) = iy
                      end if
                  if ( g(ix,iy-1) .eq. 0 ) then
                         ptr = ptr + 1
                         if ( ptr .gt. lmax ) goto 555
                         list(modid,ptr,1) = ix
                         list(modid,ptr,2) = iy-1
                      end if
                  if ( g(ix+1,iy) .eq. 0 ) then
                         ptr = ptr + 1
                         if ( ptr .gt. lmax ) goto 555
                         list(modid,ptr,1) = ix+1
                         list(modid,ptr,2) = iy
                      end if
                else
                  g(ix,iy) = 1
               end if           
             end if
 500      continue          

c         we now check if the cluster growth has already
c         terminated. If this is the case then we go back to
c         square one and start all over again

          if ( ptr .ne. 0 ) then

              nsurf = ptr
              ptr   = 0
c
c             toggle the list id
              id = mod(id+1,2)
            else
c             well, nice try but did not produce the desired result
              goto 553

          end if

 400     continue
c

c        well, nothing did happen, so we are in business
         goto 600
         
c        here are the stop or go conditions

 553     continue
         write(*,*) 'we are sorry! '
         write(*,*) 'we will start all over again. growth terminated '

 111     continue

 554     continue
         write(*,*) 'sorry, we reached our limits '
         write(*,*) 'will now go on '
         goto 600

 555     continue
         write(*,*) 'sorry, the stack overflow on L I S T'
         write(*,*) 'will go on with the result reached so far'
c

 600     continue
c
c        here will come all the rest of the stuff. the cluster
c        is now ready for action
c
         write(*,*) ' =============================================== '
         write(*,*) ' '
         write(*,*) ' results for the cluster is as follows :   '
         write(*,*) ' '
         write(*,*) ' radius ( max distance )  = ', radius
         write(*,*) ' mass                     = ', mass
         write(*,*) ' '
         write(*,*) ' =============================================== '

c        now we view the damm thing


         write(*,*) 'starting with the random walks ...'
c
c======================================================================
c
c        here we come: these boots are made for walking
c
c======================================================================
c
c

       do 5010  iwalk=1,numwal

c
c        generation of new walker
c
        try     = 0
        x0      = gmax
        y0      = gmax
        isearch = max0(ix,iy)
        do  5011  while ( (g(x0,y0) .ne. 2) .and. (try .lt. maxtry))

          x0  = (2.0 * rand() - 1.0) * isearch 
          y0  = (2.0 * rand() - 1.0) * isearch 
          try = try + 1

5011      continue

        if ( try .eq. maxtry ) then
           write(*,*) ' we are absolutly sorry, but we could not '
           write(*,*) ' place a walker on the random lattice     '
           stop(50)
        end if

        x=x0
        y=y0
        i=x0
        j=y0

c
c       now that we have placed a walker we let him loose 
c

      call R250(mz,rana,nr)

      do  5020  time=1,tmax
c
c      dichotomic waiting time distribution
c       
        itry = rand()*2
        if (itry.eq. 0)  then

c       select a direction and try to go to next site
        k = rana(time) * 4
        if(k.eq.0) i=i+1
        if(k.eq.1) i=i-1
        if(k.eq.2) j=j+1
        if(k.eq.3) j=j-1
c
c       check if the site is occupied, if not then we will try
c       again
        if(g(i,j).eq.2) then
           x=i
           y=j
         else
           i=x
           j=y
        end if
c
c
c       now we have to check if the walker did fall off the edge
       if( iabs(x) .eq. gmax-1) goto 5010
       if( iabs(y) .eq. gmax-1) goto 5010


       end if


c      here we do the bookeeping
       r2(time) = r2(time) + (x-x0)*(x-x0) + (y-y0)*(y-y0)
       count(time)=count(time)+1

 5020  continue

 5010  continue
c
c      we have done it!!!
c
c------now we will save the results to read them in perhaps
c      later on
       write(*,*) 'saving the rwalk results for later use ...'
       write(10,*) seed
       do  6000  i=1,tmax
          write(10,*) r2(i),count(i)
 6000  continue
       write(*,*) 'done '

c====================================================================
c               F I N I S H E D
c===================================================================
c
c      write(*,*) 'bye bye'
       write(*,*) 'timing is : ',mclock()
c
         end
      SUBROUTINE SUN(ISEED,RANF,N)
C
C     Stores N real random numbers in RANF, ISEED serves as
C     startup value.
C
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
C
      INTEGER ISEED
      REAL RANF(N)
C
C                            --------------------
C
      DATA FACTOR /41475557.0D0/, TWO28 /268435456.0D0/
C
C                            --------------------
C
      IF (ISEED .GE. 0) THEN
      R=DBLE(ISEED)/TWO28
      R=DMOD(R*FACTOR,1.0D0)
      ISEED=-1
      ENDIF
C
      DO 100 I = 1, N
      R=DMOD(R*FACTOR,1.0D0)
      RANF(I) = SNGL(R)
 100  CONTINUE
C
      RETURN
      END
C
C     -----------------------------------------------------------------
C
      SUBROUTINE INR250(MZ,ISEED)
C
C     Initialization routine for R250 random number generator
C              ----- For 32 bit computers only -----
C     Uses SUN for startup
C
C     VP100 option:
C *** CALL RANU2(ISEED,RRANF,250,ICON)
C     instead of
C     CALL SUN(ISEED,RRANF,250)
C
      INTEGER IUF(32),IGF(32),MZ(250)
      INTEGER IADDR1(1000),IADDR2(1000)
      REAL RRANF(2000)
C
C                            --------------------
C
      DO 100 I = 1, 30
      IGF(I) = 2**(I-1)
      IUF(I) = 2**I - 1
 100  CONTINUE
C
      IGF(31) = 2**30
      IGF(32) = - 2**31
      IUF(31) = 2**30 + (2**30 - 1)
      IUF(32) = -1
C
      CALL SUN(ISEED,RRANF,250)
C
C     Map the reals to integers -1, 0, 1
C
      DO 200 I = 1, 250
      MZ(I) = IDNINT(DBLE(2.*RRANF(I)-1.)*DBLE(FLOAT(IUF(31))))
 200  CONTINUE
C
C     Special treatment of the first 32 numbers
C
      DO 300 I = 1, 32
      MZ(I) = IOR(MZ(I),IGF(I))
      MZ(I) = IAND(MZ(I),IUF(I))
 300  CONTINUE
C
C     Make all numbers positive or zero
C
      DO 400 I = 1, 250
      MZ(I) = IAND(MZ(I),IUF(31))
 400  CONTINUE
C
C     Mixing
C
      CALL SUN(ISEED,RRANF,2000)
C
      DO 500 I = 1, 1000
      IADDR1(I) = 1 + 250.0*RRANF(I)
      IADDR2(I) = 1 + 250.0*RRANF(I+1000)
 500  CONTINUE
C
      DO 600 I = 1, 1000
      NTEMP = MZ(IADDR1(I))
      MZ(IADDR1(I)) = MZ(IADDR2(I))
      MZ(IADDR2(I)) = NTEMP
 600  CONTINUE
C
C     Warming up
C
      CALL R250(MZ,RRANF,2000)
C
      RETURN
      END
C
C     -----------------------------------------------------------------
C
      SUBROUTINE R250(MZ,RAN,N)
C
C     Stores N real random numbers in RAN using the integers in MZ
C
      INTEGER MZ(147,2)
      REAL RAN(N)
C
      DATA KK/1/
      static kk
C
C                            --------------------
C
      RMOD = 2147483647.0
C
C     Another possible choice is
C     RMOD = 524287.0
C
      K  = KK
      I  = 0
      L  = N+250
C
 10   L  = L-(251-K)
C
      DO 100 K=K, MIN0(L,147)
      MZ(K,1) = IEOR(MZ(K,1),MZ(K+103,1))
      I = I+1
      RAN(I) = MZ(K,1)/RMOD
 100  CONTINUE
C
      DO 200 K=K, MIN0(L,250)
      MZ(K-147,2) = IEOR(MZ(K-147,2),MZ(K-147,1))
      I = I+1
      RAN(I) = MZ(K-147,2)/RMOD
 200  CONTINUE
C
      IF  (K .EQ. 251) THEN
      K = 1
      GOTO 10
      ENDIF
C
      KK = K
      RETURN
      END