c===========================================================
c     dla: 2d diffusion-limited-aggregation with option 
c     for "central bias" to accelerate cluster growth.
c===========================================================
c     usage: dla2d <size> <npart> [<r0> <bias>]
c
c     size:   Number of lattice sites on a side of the 
c             arena.
c     npart:  Number of particles to evolve. Each particle
c             is evolved until one of its eight NN is 
c             fixed.  It then becomes fixed and a new 
c             particle is launched.
c     r0:     Relative launch diameter for new particles 
c             (fraction of size).
c     bias:   0 <= bias <= 1.  Amount of bias towards 
c             center of arena.  Current default is no 
c             bias.
c
c     Program reads initial fixed particle positions 
c     (x_i,y_i) from standard input (two numbers per line) 
c     and writes final fixed particle positions to standard 
c     output in same format.
c
c     The evolution (update) per se is done in-place in
c     the main-program, but separate routines for reading and 
c     writing state, generating an initial particle position,
c     and generating a random move have been coded.
c
c     Refer to class notes for further details.
c===========================================================
      program         dla

      implicit        none

      integer         iargc,        i4arg,
     &                randstep
      real*8          r8arg,        rand
      real*8          r8_never
      parameter     ( r8_never = -1.0d-60 )

      integer         maxsize
      parameter     ( maxsize = 2000 )
      integer         arena(maxsize,maxsize)

      integer         size,         npart
      real*8          r0,           bias
      real*8          default_r0,   default_bias
      parameter     ( default_r0 = 0.75d0,
     &                              default_bias = 0.0d0 )

      real*8          threshold
      integer         nfixed
      integer         xfree,        yfree,       ipart,
     &                nstep,        i,           j,
     &                dxfree,       dyfree
      real*8          theta
      logical         free 

      logical         ltrace
      parameter     ( ltrace = .true. )

c-----------------------------------------------------------
c     Argument parsing.
c-----------------------------------------------------------
      if( iargc() .lt. 2 ) go to 900
      size = i4arg(1,-1)
      if( size .lt. 50  .or.  size .gt. maxsize ) then   
         write(0,*) 'dla: Specify size between 50 and ',
     &              maxsize
         stop
      end if
      npart = i4arg(2,-1)
      if( npart .lt. 1 ) go to 900
      r0    = r8arg(3,default_r0)
      bias  = r8arg(4,default_bias)
      threshold = 1.0d0 - bias
      if( ltrace ) then
         write(0,*) 'dla: size   ', size
         write(0,*) 'dla: npart  ', npart
         write(0,*) 'dla: r0     ', r0 
         write(0,*) 'dla: bias   ', bias
      end if

c-----------------------------------------------------------
c     Initialize arena and read fixed particle positions
c     from standard input.
c-----------------------------------------------------------
      call getfixed(arena,maxsize,size,nfixed)
      if( nfixed .gt. 0 ) then
         write(0,*) 'dla: Read ', nfixed, 
     &              ' particle positions.'
      else
         write(0,*) 'dla: No valid fixed particle ',
     &              'positions read.  Exiting.'
         stop
      end if

c-----------------------------------------------------------
c     For number of requested particles ...
c-----------------------------------------------------------
      do ipart = 1 , npart
         nstep = 0 
c-----------------------------------------------------------
c        Generate random initial position.
c-----------------------------------------------------------
         call initposrand(size,r0,xfree,yfree)
c-----------------------------------------------------------
c        Until the particle is fixed ... 
c-----------------------------------------------------------
         free = .true.
         do while ( free )
c-----------------------------------------------------------
c           Take a random step (+1,0,-1) in both directions.
c-----------------------------------------------------------
            xfree = xfree + randstep()
            yfree = yfree + randstep()
c-----------------------------------------------------------
c           If bias is non zero, take a step towards 
c           the origin with probability 'bias'.
c-----------------------------------------------------------
            if( bias .ne. 0.0d0 ) then
               dxfree = xfree - (size - 1) / 2
               dyfree = yfree - (size - 1) / 2
               theta  = atan2(1.0d0 * dxfree,1.0d0 * dyfree)
               if( rand() .le. abs(sin(theta)) ) then
                  if( rand() .gt. threshold ) then
                     if( dxfree .gt. 0 ) then
                        xfree = xfree - 1
                     else 
                        xfree = xfree + 1
                     end if
                  end if
               end if
               if( rand() .le. abs(cos(theta)) ) then
                  if( rand() .gt. threshold ) then
                     if( dyfree .gt. 0 ) then
                        yfree = yfree - 1
                     else 
                        yfree = yfree + 1
                     end if
                  end if
               end if
            end if
c-----------------------------------------------------------
c           Check if particle is outside arena.
c-----------------------------------------------------------
            if( xfree .lt. 1  .or.  xfree .gt. size  .or.
     &          yfree .lt. 1  .or.  yfree .gt. size ) then
c-----------------------------------------------------------
c              If it is, reinitialize
c-----------------------------------------------------------
               call initposrand(size,r0,xfree,yfree)
            end if
c-----------------------------------------------------------
c           Check if particle should be fixed.
c-----------------------------------------------------------
            do i = max(xfree-1,1) , min(xfree+1,size)
               do j = max(yfree-1,1) , min(yfree+1,size)
c-----------------------------------------------------------
c                 If it is, update corresponding arena 
c                 site and set flag.
c-----------------------------------------------------------
                  if( arena(i,j) .ne. 0 ) then      
                     arena(xfree,yfree) = 1
                     free = .false.
                  end if
               end do
            end do
            nstep = nstep + 1
         end do
         write(0,*) 'dla: Particle ', ipart, ' fixed after ',
     &              nstep, ' steps'
      end do

c-----------------------------------------------------------
c     Write fixed particle positions to standard output.
c-----------------------------------------------------------
      call putfixed(arena,maxsize,size,nfixed)
      write(0,*) 'dla: Wrote ', nfixed, ' particle positions.'

      stop

900   continue
         write(0,*) 'usage: dla <size> <npart> [<r0> <bias>]'
         write(0,9000) default_r0, default_bias
9000     format(/  
     &      '        Current default <r0>: ',f13.2/
     &      '        Current default <bias>: ',1p,e11.4,0p//
     &      '        Program reads initial fixed-particle coordinates '/
     &      '        (integers x,y; 1 <= x,y <= size)  from standard'/
     &      '        input, writes final fixed positions to standard'/
     &      '        output.')
      stop 

      end

c===========================================================
c     Returns -1, 0 or 1 chosen randomly
c===========================================================
      integer function randstep()
         implicit     none

         real*8       rand

         randstep = min(2,int(3.0d0 * rand())) - 1

         return
      end 

c===========================================================
c     Initialize arena then read fixed particle positions 
c     from standard input.  Ignore particles lying outside 
c     of current arena.  Returns number of fixed particles 
c     inside arena.
c===========================================================
      subroutine getfixed(arena,maxsize,size,nfixed)
         implicit     none
         
         integer      maxsize,      size,       nfixed
         integer      arena(maxsize,maxsize)

         integer      x,            y,          rc,
     &                i,            j

         do j = 1 , size
            do i = 1 , size
               arena(i,j) = 0
            end do
         end do

         nfixed = 0
 100     continue
            read(*,*,end=200,iostat=rc)  x,  y
            if( rc .eq. 0  ) then
               if( 1 .le. x  .and.  x .le. size   .and.  
     &             1 .le. y  .and.  y .le. size ) then
                  arena(x,y) = 1
                  nfixed = nfixed + 1
               end if
            end if
         go to 100

 200     continue
         return

      end

c===========================================================
c     Writes fixed particle positions to standard output.
c     Returns number of fixed particles.
c===========================================================
      subroutine putfixed(arena,maxsize,size,nfixed)
         implicit     none

         integer      maxsize,      size,       nfixed
         integer      arena(maxsize,maxsize)

         integer      i,            j

         nfixed = 0
         do j = 1 , size
            do i = 1 , size
               if( arena(i,j) .ne. 0 ) then
                  nfixed = nfixed + 1
                  write(*,*) i, j
               end if
            end do
         end do

         return
      end

c-----------------------------------------------------------
c     Generates initial particle position 0.5 * r0 * size
c     from arena center, randomly positioned in angle.
c-----------------------------------------------------------
      subroutine initposrand(size,r0,xfree,yfree)
         implicit     none

         real*8       rand
         integer      size,    xfree,     yfree
         real*8       r0
         real*8       r,       theta
c-----------------------------------------------------------
c        Generate a random angle from 0 to 2 Pi.
c-----------------------------------------------------------
         theta = rand() * 8.0d0 * atan(1.0d0)
         r     = 0.5d0 * r0 * (size - 1)
         xfree = 0.5d0 * (size - 1)  + r * cos(theta)
         yfree = 0.5d0 * (size - 1)  + r * sin(theta)
         return
      end