      program walk
      implicit none
      include 'walk.cm'

      integer i,j,n,m,w,delaytimer,oldnfood,a,b,foodtimer,oldval
     &  ,length,viewtimer,temp,oldeat,newval,decayamt,addamt
      real*8 decaymult,x,random,foodmult,decaytimer,writetimer
      character*45 infile

c read in base runid
      write(*,*) 'Enter base runid (less than 30 characters): '
      read(*,*) runid
      length = index(runid,' ')-1

c open input file and read in parameters
      infile = runid(1:length)//'.in'
      call readinput(infile)

c get output filenames
      dfile = runid(1:length)//'.d'
      efile = runid(1:length)//'.e'
      outfile = runid(1:length)//'.out'
      psfile = runid(1:length)//'.ps'
      vfile = runid(1:length)//'.v'
      fdfile = runid(1:length)//'.fd'
      grfile = runid(1:length)//'.gr'

c begin the run!

      seed = initseed
      call initialize           ! initializes arrays and variables
      call getnbrs              ! calculates neighbor sites, with PBCs
      call newwalker            ! add a walker

      delaytimer = 0
      decaymult = 1. - 1./tau      
      foodmult = 1. - 1./fooddecay      
      foodtimer = 0
      decaytimer = 0
      writetimer = 0
      viewtimer = 0
      oldeat = 0
      decayamt = 0
      addamt = 0

      write(*,*) ''
      write(*,*) 'fooddecay =',fooddecay
      write(*,*) 'foodmult =',foodmult
      i = foodmult*foodpersite
      write(*,*) 'foodmult*foodpersite =',i
      write(*,*) ''
      write(*,*) 'tau =',tau
      write(*,*) 'pheromone decay:'
      write(*,*) '1.0 -->',decaymult,' after 1 step'
      write(*,*) '    -->',decaymult**10,' after 10 steps'
      write(*,*) '    -->',decaymult**20,' after 20 steps'
      write(*,*) 'foodpersite =',foodpersite,' (# of food ptcles/chunk)'

      call openapp(lun(4),outfile)

      a = foodrate/foodpersite  ! # of food chunks per step
      b = foodpersite/foodrate  ! steps per chunk
      if (a.ge.1) then ! add food every step
       write(*,*) 'Adding',a,' chunks of food to the system every step'
       write(lun(4),*) 'Adding',a,' chunks of food every step'
      else
       write(*,*) 'Adding 1 chunk of food to the system every',b
     &  ,' steps'
       write(lun(4),*) 'Adding 1 chunk of food to the system every',b
     &  ,' steps'
      endif

      close(lun(4))

      oldnfood = foodeaten

cccccccccccccccccccccccc  MAIN LOOP  ccccccccccccccccccccccccccccc

      write(*,*) ''
      write(*,*) 'Running',msteps,' steps...'
      do n=1,msteps

c pheromone decay
        do m=1,msites
          pheromone(m,1) = decaymult*pheromone(m,1)
        enddo
        if (npher.gt.1) then
          do m=1,msites
            pheromone(m,2) = decaymult*pheromone(m,2)
          enddo
        endif


cewd
c        write(*,*) 'About to call movewalkers...',n

c move walkers
        call movewalkers
        if (nofood.eq.1.and.mfood.gt.0) then 
          write(*,*) 'food exhausted at step',n
          call openapp(lun(4),outfile)
          write(lun(4),*) 'food exhausted at step',n
          close(lun(4))
          goto 99
        endif

c add a walker?
        delaytimer = delaytimer + 1
        if (delaytimer.ge.wdelay .and. nwalkers.lt.mwalkers) then
          delaytimer = 0
          call newwalker
        endif

c add food at random
        a = foodrate/foodpersite ! # of food chunks per step
        b = foodpersite/foodrate ! steps per chunk
        foodtimer = foodtimer + 1
        if (a.ge.1) then        ! add food every step
          do i=1,a
            x = nsites*random(1)+1
            food(x) = food(x) + foodpersite + random(1)
            nfood = nfood + foodpersite
            addamt = addamt + foodpersite   ! keep track for output
          enddo
        else
          if (foodtimer.ge.b) then
            foodtimer = 0
            x = nsites*random(1)+1
            food(x) = food(x) + foodpersite + random(1)
            nfood = nfood + foodpersite
            addamt = addamt + foodpersite  ! keep track for output
          endif
        endif

c food decay
        decaytimer = decaytimer + 1
        if (decaytimer.ge.decaydelay) then
          decaytimer = 0
          do m=1,nsites
            oldval = food(m)
            newval = food(m)*foodmult

cewd
c            if (newval.gt.oldval) then
c              write(*,*) 'Here is the problem!',newval,oldval
c            endif

            if (newval.ne.oldval) then 
              decayamt = decayamt + (oldval-newval)
            endif
            food(m) = food(m)*foodmult
            nfood = nfood - (oldval-food(m)*foodmult)
          enddo
        endif

c compute nfood (ignore fractional amounts of food on each site)
        temp = 0
        do m=1,nsites
          temp = temp + food(m)
        enddo
        nfood = temp        

c calculate diffusion constant

        if (mod(n,dwrite).eq.0) then
          call diffusion(n)
          call efficiency(n,oldnfood)
          call gofrcalc
          oldnfood = foodeaten
        endif


c write out food supply
        writetimer = writetimer + 1
        if (writetimer.ge.writefreq) then
          writetimer = 0
          write(*,*) 'step',n,', nfood =',nfood,'  ('
     &       ,addamt,' added,',foodeaten-oldeat,' eaten,'
     &       ,decayamt,' decayed )'
          call openapp(lun(4),outfile)
          write(lun(4),*) 'step',n,', nfood =',nfood,'  ('
     &       ,addamt,' added,',foodeaten-oldeat,' eaten,'
     &       ,decayamt,' decayed )'
          close(lun(4))
          call openapp(lun(6),fdfile)
          write(lun(6),*) n,nfood
          close(lun(6))
          oldeat = foodeaten
          decayamt = 0
          addamt = 0
        endif

c Viewer output
        if (viewerout.eq.1.and.n.gt.viewdelay) then
          viewtimer = viewtimer + 1
          if (n.eq.1) then
            call openapp(lun(5),vfile)
            do w=1,mwalkers
              write(lun(5),60) n,0,0
            enddo
            write(lun(5),*) ''
            close(lun(5))
          endif

          if (viewtimer.ge.viewfreq.and.n.gt.1) then
            viewtimer = 0 
            call openapp(lun(5),vfile)
            do w=1,nwalkers
              call conv(walker(w),i,j)
              write(lun(5),60) n,i-1,j-1
 60           format(i8,2i10)
            enddo
            write(lun(5),*) ''
            close(lun(5))
          endif
        endif

      enddo   ! nsteps

 99   continue


cccccccccccccccccccccccc  OUTPUT SECTION  cccccccccccccccccccccccccc


      call wgofr

c      write(*,*) 'nest sites: ',(nestid(i),i=1,nnest)
c      write(*,*) ''
      write(*,*) 'number of food particles left: ',nfood
      call openapp(lun(4),outfile)
      write(lun(4),*) 'number of food particles left: ',nfood
      close(lun(4))

c color in special sites
c   set array site input to pswrite (e.g. nhits) to: 
c      -10 for red
c      -20 for yellow
c      -30 for green
c      -40 for cyan
c      -50 for blue
c      -60 for magenta
      do i=1,nnest
        nhits(nestid(i)) = -40
      enddo

      if (showwalkers.eq.1) then
        do i=1,mwalkers
          nhits(walker(i)) = -10
        enddo
      endif

c graphical output section

      if (write_ps.eq.1) then
        call pswrite(nhits,nsites)
        write(*,*) 'Writing output to ',psfile
      endif
      
      end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccc     SUBROUTINES    ccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      subroutine readinput(infile)
      implicit none
      include 'walk.cm'

      integer i,len,ncount
      real*8 value
      character*45 infile,name

c set default values
      mfood = 1000
      msteps = 1000
      dwrite = msteps/20
      writefreq = 50
      tau = 20
      wdelay = 1
      lattype = 1
      wint = 1.0
      wint2 = 1.0
      foodpersite = 50
      outres = 800
      initseed = 2002
      write_ps = 0
      showwalkers = 0
      foodrate = 3.0
      fooddecay = 100
      decaydelay = 25
      showfood = 0
      viewerout = 0
      viewfreq = 1000
      npher = 1
      nscouts = mwalkers/2
      binwidth = 1.
      viewdelay = 0
      inertia = 0.

c read input file
      len = index(infile,' ')-1
      open(20,file=infile(1:len),status='old',form='formatted')
      ncount = 0
      do i=1,100000
        read(20,*,end=91) value,name
        ncount = ncount+1
        len=index(name,' ')-1

        if(name(1:len).eq.'mfood') then
          mfood = value
        elseif (name(1:len).eq.'msteps') then
          msteps = value
        elseif(name(1:len).eq.'tau') then
          tau = value
        elseif(name(1:len).eq.'wdelay') then
          wdelay = value
        elseif(name(1:len).eq.'lattype') then
          lattype = value
        elseif(name(1:len).eq.'wint') then
          wint = value
        elseif(name(1:len).eq.'wint2') then
          wint2 = value
        elseif(name(1:len).eq.'foodpersite') then
          foodpersite = value
        elseif(name(1:len).eq.'outres') then
          outres = value
        elseif(name(1:len).eq.'dwrite') then
          dwrite = value
        elseif(name(1:len).eq.'initseed') then
          initseed = value
        elseif(name(1:len).eq.'write_ps') then
          write_ps = value
        elseif(name(1:len).eq.'showwalkers') then
          showwalkers = value
        elseif(name(1:len).eq.'foodrate') then
          foodrate = value
        elseif(name(1:len).eq.'fooddecay') then
          fooddecay = value
        elseif(name(1:len).eq.'decaydelay') then
          decaydelay = value
        elseif(name(1:len).eq.'writefreq') then
          writefreq = value
        elseif(name(1:len).eq.'showfood') then
          showfood = value
        elseif(name(1:len).eq.'viewerout') then
          viewerout = value
        elseif(name(1:len).eq.'viewfreq') then
          viewfreq = value
        elseif(name(1:len).eq.'npher') then
          npher = value
        elseif(name(1:len).eq.'nscouts') then
          nscouts = value
        elseif(name(1:len).eq.'binwidth') then
          binwidth = value
        elseif(name(1:len).eq.'viewdelay') then
          viewdelay = value
        elseif(name(1:len).eq.'greynorm') then
          greynorm = value
        elseif(name(1:len).eq.'inertia') then
          inertia = value
        elseif(name(1:len).eq.'end') then 
          goto 91
        else
          write(*,*) 'WARNING!!  unrecognized parameter: ',name(1:len)
        endif
      enddo
 91   close(20)

      write(*,*) ncount,' lines read from ',infile

      return
      end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      subroutine initialize
      implicit none
      include 'walk.cm'

      integer i,j,x,nestsite,cnt,center,counter,npixels
      real*8 random,a

c initialize random number generator
      a = random(0)

c set up luns
      do i=1,10
        lun(i) = 30+i
      enddo

c open viewer file
      if (viewerout.eq.1) then
        open(lun(5),file=vfile)
        write(lun(5),'(a8,i8,a8)') '#PERIOD ',viewfreq,'   steps'
        write(lun(5),'(a)') '#   step     walker_x     walker_y'
        write(*,*) 'AntViewer walker information will be stored in:'
     &    ,vfile
        close(lun(5))
      endif

c open food file (nfood vs. step)
      open(lun(6),file=fdfile)
      write(lun(6),'(a8,i8,a8)') '#PERIOD ',writefreq,'   steps'
      write(lun(6),'(a)') '#NAMES  time  nfood'
      write(lun(6),'(a)') '#UNITS  steps   null'
      write(lun(6),*) 0,mfood
      write(*,*) 'Food information will be stored in:',fdfile
      close(lun(6))

c open diffusion constant file
      cumd_w = 0.
      open(lun(2),file=dfile)
      write(lun(2),'(a8,i8,a8)') '#PERIOD ',dwrite,'   steps'
      write(lun(2),'(a)') '#NAMES  time  diffusion   cum_diffusion'
      write(lun(2),'(a)') '#UNITS  steps   sites^2/step   sites^2/step'
      write(*,*) 'Diffusion rate information will be stored in:',dfile
      close(lun(2))

c open food collection efficiency file
      cumeff = 0.
      open(lun(3),file=efile)
      write(lun(3),'(a8,i8,a8)') '#PERIOD ',dwrite,'   steps'
      write(lun(3),'(a)') '#NAMES  time  efficiency   cum_effic'
      write(lun(3),'(a)') '#UNITS  steps   food/step    food/step'
      write(*,*) 'Food collection efficiency info will be stored in:'
     & ,efile
      close(lun(3))

c open output file
      open(lun(4),file=outfile)
      write(lun(4),*) 'seed =',seed
      write(lun(4),*) 'msites =',msites
      write(lun(4),*) 'mwalkers =',mwalkers
      write(lun(4),*) 'mfood =',mfood
      write(lun(4),*) 'msteps =',msteps
      write(lun(4),*) 'tau =',tau
      write(lun(4),*) 'wdelay =',wdelay
      write(lun(4),*) 'lattype =',lattype
      write(lun(4),*) 'wint =',wint
      write(lun(4),*) 'wint2 =',wint2
      write(lun(4),*) 'foodpersite =',foodpersite
      write(lun(4),*) 'outres =',outres
      write(lun(4),*) 'dwrite =',dwrite
      write(lun(4),*) 'foodrate =',foodrate
      write(lun(4),*) 'fooddecay =',fooddecay
      write(lun(4),*) 'decaydelay =',decaydelay
      write(lun(4),*) 'writefreq =',writefreq
      write(lun(4),*) 'nscouts =',nscouts
      write(lun(4),*) 'npher =',npher
      write(lun(4),*) 'binwidth =',binwidth
      close(lun(4))

c initialize various variables
      grnorm = 0
      nwalkers = 0
      do i=1,mwalkers
        lastsite(i) = 0
        timeout(i) = 0
        hasfood(i) = 0
        scout(i) = 0
      enddo
      nofood = 0
      foodeaten = 0

c make scouts
      if (nscouts.gt.mwalkers) then
        write(*,*) 'nscouts > mwalkers!!'
        call openapp(lun(4),outfile)
        write(lun(4),*) 'nscouts > mwalkers!!'
        close(lun(4))
        stop
      endif

      if (npher.gt.1) then
        do i=1,nscouts
          scout(i) = 1
        enddo
      endif

c calculate the number of sites
      n1d = sqrt(1.0*msites)  ! integer truncation
      nsites = n1d*n1d    ! actual number of sites (may not = msites)

c initialize lattice

      do i=1,msites
        nhits(i) = 0
        food(i) = 0.
        initfood(i) = 0
        pheromone(i,1) = 0.
        pheromone(i,2) = 0.
        nest(i) = 0
      enddo

      nfood = mfood    ! total number of food
      write(*,*) 'Total number of food particles (t=0): ',nfood

c make a nest ! we'll hardwire it as 7x7 for now
      nnest = 49
      center = (n1d+1)/2 + nsites/2
      cnt = 0
      do i=-3,3,1
        do j=-3,3,1
          nestsite = center + i + j*n1d
          nest(nestsite) = 1
          cnt = cnt+1
          nestid(cnt) = nestsite
        enddo
      enddo

c place food
      if (mfood.eq.0) goto 2
      counter = nfood
      do i=1,nsites*nsites         ! essentially an infinite loop
        x = nsites*random(1)+1      ! pick a site at random
        if (food(x).lt.foodpersite.and.nest(x).ne.1) then
          food(x) = foodpersite + random(1)    ! place food on x
          initfood(x) = 1                      ! (random just staggers
          counter = counter-foodpersite        !  init decay cycle)
        endif
        if (counter.lt.foodpersite) goto 2  ! quit when food exhausted
      enddo
 2    continue

c set up bins for g(r)

      call binsites

c set up graphics output

      if (write_ps.eq.1) then
        npixels = outres/n1d    ! number of pixels per lattice site
        call psinit(npixels,lun(1),psfile)
      endif

      return
      end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine newwalker
      implicit none
      include 'walk.cm'

      integer i,j,w,freesite
      real*8 random

c make sure nest sites aren't all full
      do j=1,25
        freesite = 1
        i=nestid(j)
        do w=1,nwalkers
          if (walker(w).eq.i) then
            freesite = 0
          endif
        enddo
        if (freesite.eq.1) then   ! there is a free site
          goto 5
        endif
      enddo
      
      goto 6  ! there isn't a free site -- try again later

c pick a random nest site
 5    continue
      j = nnest*random(1)+1
      i = nestid(j)
      do w=1,nwalkers                ! make sure no one is on site i
        if (walker(w).eq.i) goto 5
      enddo

      if (nwalkers.lt.mwalkers) then
        nwalkers = nwalkers + 1
        walker(nwalkers) = i
        lastsite(nwalkers) = i
        timeout(nwalkers) = 0
      endif

 6    continue

      return
      end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine movewalkers
      implicit none
      include 'walk.cm'

      integer i,j,n,w,y,z,site,newsite,foodflag,siteid(6),inertflag
      integer a,b,c,d
      real*8 prob,cump(0:6),cumf(0:6),sample
      real*8 random

      if (npher.eq.1) then   ! SINGLE PHEROMONE

        do w=1,nwalkers

c calculate the relative probability of all possible moves
          site = walker(w)      ! current location of walker
          cump(0) = 0.
          cumf(0) = 0.
          foodflag = 0
          inertflag = 1
          do j=1,nnbr
            z = nbr(j,site)     ! neighboring site
            prob = 1. + wint*pheromone(z,1)
            
            do n=1,w-1         ! don't move on top of another walker
              if (walker(n).eq.z) prob = 0.
            enddo
            if (prob.lt.0.) prob = 0.         ! no negative probability
            if (z.eq.lastsite(w)) prob = 0.   ! can't backtrack
            if (food(z).ge.1.) then            ! food always wins
              foodflag = foodflag + 1
              siteid(foodflag) = z
              cumf(foodflag) = cumf(foodflag-1) + food(z)
            endif

c inertia?
            if (inertia.gt.0.) then
              if (z.ne.lastsite(w)) then
                if (pheromone(z,1).gt.0.) inertflag = 0
                if (foodflag.gt.0) inertflag = 0
                if (walker(n).eq.z) inertflag = 0
                if (nest(z).eq.1) inertflag = 0
              endif
            endif

            cump(j) = cump(j-1) + prob        ! add to cumulative prob
          enddo

c for inertial moves, re-compute the cumulative probability
          if (inertia.gt.0.) then
            if (inertia.gt.1.) then
              write(*,*) 'inertia > 1!',inertia
              stop
            endif
            if (inertflag.gt.0) then
              call conv(lastsite,a,b)
              do j=1,nnbr
                cump(j) = 0.
                z = nbr(j,site)
                if (z.ne.lastsite(w)) then
                  call conv(z,c,d) 
                  if (a.eq.c.or.b.eq.d) then  ! this is the forward site
                    prob = 2.*inertia/(1.-inertia)
                  else
                    prob = 1.
                  endif
                else
                  prob = 0.
                endif
                cump(j) = cump(j-1) + prob
              enddo
            endif
          endif

c make sure cump > 0
          if (cump(nnbr).le.0.) then
            do j=1,nnbr
              cump(j) = cump(j-1) + 1.0
            enddo
          endif
              
c sample the probability

c just to be safe, pick a default site
          newsite = nbr(1,site)

c now sample

c if food was found, then the probability is weighted by how much food
c is on each neighboring site
          if (foodflag.gt.0) then
            sample = random(1)*cumf(foodflag)
            do j=1,foodflag
              if (cumf(j).gt.sample) then
                newsite = siteid(j)
                sample = cumf(foodflag) + 1.     ! end the loop
              endif
            enddo

c if no food was found, then the probability is weighted by how much
c pheromone was found on each site
          else
            sample = random(1)*cump(nnbr)
            do j=1,nnbr
              if (cump(j).gt.sample) then
                newsite = nbr(j,site)
                sample = cump(nnbr) + 1.        ! end the loop
              endif
            enddo
          endif

c move the walker

          if (foodflag.gt.0) then

cewd
            if (newsite.gt.nsites) then
              write(*,*) 'newsite > nsites!',newsite,nsites
            endif


c            if (food(newsite).lt.1.) then
c              write(*,*) 'Trying to drink tea from an empty cup!'
c     &           ,newsite,food(newsite)
c            endif
cewd
            nhits(newsite) = nhits(newsite) + 1
            pheromone(site,1) = pheromone(site,1) + 1.
            food(newsite) = food(newsite) - 1.       ! YUMMY
            nfood = nfood - 1 
            foodeaten = foodeaten + 1
            if (nfood.eq.0) nofood = 1

c teleport walker back to the nest
 55         continue
            j = nnest*random(1)+1
            i = nestid(j)
            do y=1,nwalkers     ! make sure no one is on site i
              if (walker(y).eq.i) goto 55
            enddo
            walker(w) = i
            lastsite(w) = i
            timeout(w) = 0

          else
            walker(w) = newsite
            pheromone(site,1) = pheromone(site,1) + 1.
            nhits(newsite) = nhits(newsite) + 1
            lastsite(w) = site
            timeout(w) = timeout(w) + 1
          endif

        enddo

      else   ! MULTIPLE PHEROMONES

        do w=1,nwalkers

c calculate the relative probability of all possible moves
          site = walker(w)      ! current location of walker
          cump(0) = 0.
          cumf(0) = 0.
          inertflag = 1
          foodflag = 0
          do j=1,nnbr
            z = nbr(j,site)     ! neighboring site

c            if (scout(w).eq.1.and.hasfood(w).eq.0) then
c              prob = 1.      ! random walk
c            elseif (scout(w).eq.1.and.hasfood(w).eq.1) then
c              prob = 1. + wint*pheromone(z,1)
c            else
c              prob = 0.01 + wint2*pheromone(z,2)   ! worker
c            endif

            if (scout(w).eq.1.and.hasfood(w).eq.0) then
              prob = 0.1 - wint*pheromone(z,1)  ! scout, no food
            elseif (scout(w).eq.1.and.hasfood(w).eq.1) then
              prob = 0.1 + wint*pheromone(z,1)  ! scout, w. food
            elseif(scout(w).ne.1.and.hasfood(w).eq.0) then
              prob = 0.1 + wint2*pheromone(z,2) ! worker, no food
            else
              prob = 0.1 + wint2*pheromone(z,1)  ! worker, w. food
            endif


            do n=1,w-1          ! don't move on top of another walker
              if (walker(n).eq.z) prob = 0.
            enddo
            if (prob.lt.0.) prob = 0.         ! no negative probability
c            if (z.eq.lastsite(w).and.scout(w).eq.1) then
c              prob = 0.                       ! scouts can't backtrack
c            endif
            if (z.eq.lastsite(w)) then
              prob = 0.                    ! no backtracking
            endif
            if (hasfood(w).eq.1.and.nest(z).eq.1) then
              prob = 10000000.                ! put food in nest
            elseif (nest(z).eq.1.and.scout(w).ne.1) then
              prob = 0.2                      ! workers stay near nest
            endif

            if (food(z).ge.1..and.hasfood(w).eq.0) then   ! food wins
              foodflag = foodflag + 1
              siteid(foodflag) = z
              cumf(foodflag) = cumf(foodflag-1) + food(z)
            endif


c inertia?  (i.e. walker moves forward preferentially)
            if (inertia.gt.0.) then
              if (scout(w).eq.1) then
                if (z.ne.lastsite(w)) then
                  if (pheromone(z,1).gt.0.) inertflag = 0
                  if (foodflag.gt.0) inertflag = 0
                  if (walker(n).eq.z) inertflag = 0
                  if (nest(z).eq.1) inertflag = 0
                endif
              else
                if (pheromone(z,1).gt.0.) inertflag = 0
                if (foodflag.gt.0) inertflag = 0
                if (walker(n).eq.z) inertflag = 0
                if (nest(z).eq.1) inertflag = 0
              endif
            endif


            cump(j) = cump(j-1) + prob        ! add to cumulative prob
          enddo


c for inertial moves, re-compute the cumulative probability
          if (inertia.gt.0.) then
            if (inertia.gt.1.) then
              write(*,*) 'inertia > 1!',inertia
              stop
            endif
            if (inertflag.gt.0) then
              call conv(lastsite,a,b)
              do j=1,nnbr
                cump(j) = 0.
                z = nbr(j,site)
                if (z.ne.lastsite(w)) then
                  call conv(z,c,d) 
                  if (a.eq.c.or.b.eq.d) then  ! this is the forward site
                    prob = 2.*inertia/(1.-inertia)
                  else
                    prob = 1.
                  endif
                else
                  prob = 0.
                endif
                cump(j) = cump(j-1) + prob
              enddo
            endif
          endif

c make sure cump > 0
          if (cump(nnbr).le.0.) then
            do j=1,nnbr
              cump(j) = cump(j-1) + 1.0
            enddo
          endif
              
c sample the probability

c just to be safe, pick a default site
          newsite = nbr(1,site)

c if food was found, then the probability is weighted by how much food
c is on each neighboring site
          if (foodflag.gt.0) then
            sample = random(1)*cumf(foodflag)
            do j=1,foodflag
              if (cumf(j).gt.sample) then
                newsite = siteid(j)
                sample = cumf(foodflag) + 1.     ! end the loop
              endif
            enddo

c if no food was found, then the probability is weighted by how much
c pheromone was found on each site
          else
            sample = random(1)*cump(nnbr)
            do j=1,nnbr
              if (cump(j).gt.sample) then
                newsite = nbr(j,site)
                sample = cump(nnbr) + 1.        ! end the loop
              endif
            enddo
          endif

c move the walker

          if (foodflag.gt.0) then
            food(newsite) = food(newsite) - 1.
            hasfood(w) = 1
            lastsite(w) = site
            timeout(w) = timeout(w) + 1

          elseif (nest(newsite).eq.1.and.hasfood(w).eq.1) then
            walker(w) = newsite
            lastsite(w) = site
            nhits(newsite) = nhits(newsite) + 1
            nfood = nfood - 1 
            foodeaten = foodeaten + 1
            hasfood(w) = 0
            if (nfood.eq.0) nofood = 1
            timeout(w) = 0

          else
            walker(w) = newsite
            nhits(newsite) = nhits(newsite) + 1
            lastsite(w) = site
            timeout(w) = timeout(w) + 1
          endif

c lay down pheromone
c          if (scout(w).eq.1.and.hasfood(w).eq.0) then
c            pheromone(site,1) = pheromone(site,1) + 1.
c          elseif (scout(w).eq.1.and.hasfood(w).eq.1) then
c            pheromone(site,2) = pheromone(site,2) + 1.
c          endif
          if (hasfood(w).eq.0.and.scout(w).eq.1) then
            pheromone(site,1) = pheromone(site,1) + 1.
          elseif (hasfood(w).eq.1) then
            pheromone(site,2) = pheromone(site,2) + 1.
          endif

        enddo
      endif

      return
      end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine diffusion(n)
      implicit none
      include 'walk.cm'

      integer n,w,x,y,t,a,b,walkcnt
      real*8 mindist,dist,d_w

      d_w = 0.
      walkcnt = 0
      do w=1,nwalkers

        call conv(walker(w),x,y)   ! converts cellnumber into coords
        mindist = nsites
        do t=1,nnest
          call conv(nestid(t),a,b)
          ! find distance to nest site (can ignore PBCs, since nest is
          !                             in the center)
          dist = sqrt( 1.*(x-a)**2 + 1.*(y-b)**2)
          mindist = min(mindist,dist)
        enddo
        if (timeout(w).gt.0) then
          d_w = d_w + mindist**2/timeout(w)  ! d_w = <R^2>/nsteps  
          walkcnt = walkcnt + 1
        endif
      enddo
      if (walkcnt.gt.0) then
        d_w = d_w/walkcnt    ! get the average value

c add to cumulative average
        cumd_w = cumd_w*(1.*n-1.*dwrite)/dwrite   ! get old total
        cumd_w = cumd_w + d_w                     ! add new point
        cumd_w = cumd_w/(1.*n/dwrite)             ! get new average 

        call openapp(lun(2),dfile)
        write(lun(2),75) n,d_w,cumd_w
 75     format(i8,6X,f10.7,9X,f10.7)
        close(lun(2))
      endif

      return
      end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine conv(cell,x,y)
      implicit none
      include 'walk.cm'

      integer cell,x,y

      x = mod(cell,n1d)
      if (x.eq.0) x = n1d
      y = cell/n1d + 1

      return
      end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine efficiency(n,oldnfood)
      implicit none
      include 'walk.cm'

      integer n,oldnfood
      real*8 eff

      eff = (1.0*foodeaten-1.0*oldnfood)/dwrite  ! % food collected in 
      eff = eff/nwalkers                         ! last dwrite steps
      eff = eff/nfood

      cumeff = 1.0*foodeaten/n            ! total fraction of food 
      cumeff = cumeff/nwalkers            ! collected per step
      cumeff = cumeff/nfood

      call openapp(lun(3),efile)
      write(lun(3),76) n,eff,cumeff
 76   format(i8,6X,e13.6,5X,e13.6)
      close(lun(3))

      return
      end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine gofrcalc
      implicit none
      include 'walk.cm'

      integer b,j,site,sum

      do b=1,maxbins
        sum = 0.
        do j=1,nbins(b)
          site = binsite(j)
          sum = sum + food(site)
        enddo

        gofr(b) = gofr(b) + (1.0*sum/nbins(b))/(1.0*nfood/nsites)

      enddo
      grnorm = grnorm + 1

      return
      end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine wgofr
      implicit none
      include 'walk.cm'

      integer b
      real*8 val

      write(*,*) 'g(r) information will be stored in:',grfile
      open(lun(7),file=grfile)

      write(lun(7),'(a)') '#NAMES r g(r)'
      write(lun(7),'(a)') '#UNITS sites null'
      write(lun(7),'(a)') '#   r     g(r)'

      do b=1,maxbins
        if (nbins(b).gt.0) then
          val = gofr(b)/grnorm
          write(lun(7),39) bindist(b),val
 39       format(f12.8,5X,e13.6)
        endif
      enddo

      close(lun(7))      

      return
      end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine binsites
      implicit none
      include 'walk.cm'

      integer i,t,a,b,x,y,bin,j
      real*8 mindist,dist

      do j=1,maxbins
        nbins(j) = 0
      enddo

      do j=1,maxbins
        bindist(j) = (j-0.5)*binwidth
      enddo


      do i=1,nsites

        call conv(i,x,y)        ! converts cellnumber into coords
        mindist = nsites
        do t=1,nnest
          call conv(nestid(t),a,b)
          ! find distance to nest site (can ignore PBCs, since nest is
          !                             in the center)
          dist = sqrt( 1.*(x-a)**2 + 1.*(y-b)**2)
          mindist = min(mindist,dist)
        enddo

        bin = mindist/binwidth + 1

        if (bin.gt.maxbins) then
          write(*,*) 'error in binsites:  maxbins too small!'
          stop
        endif

        nbins(bin) = nbins(bin) + 1
        binsite(nbins(bin)) = i
      enddo


      return
      end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine getnbrs
      implicit none
      include 'walk.cm'

c calculates (once) which sites neighbor each site, using periodic 
c boundary conditions (PBCs)

      integer i,a,b
      integer unconv

      if (lattype.eq.1) then                   ! square lattice
        nnbr = 4
        do i=1,nsites
          a = mod(i,n1d)   ! Cartesian coordinates
          b = i/n1d + 1
          nbr(1,i) = unconv(a+1,b,n1d)
          nbr(2,i) = unconv(a-1,b,n1d)
          nbr(3,i) = unconv(a,b+1,n1d)
          nbr(4,i) = unconv(a,b-1,n1d)
        enddo

      elseif (lattype.eq.2) then               ! triangular lattice
        nnbr = 6
        do i=1,nsites
          a = mod(i,n1d)   ! Cartesian coordinates
          b = i/n1d + 1
          nbr(1,i) = unconv(a+1,b,n1d)
          nbr(2,i) = unconv(a-1,b,n1d)
          nbr(3,i) = unconv(a-1,b+1,n1d)
          nbr(4,i) = unconv(a,b+1,n1d)
          nbr(5,i) = unconv(a,b-1,n1d)
          nbr(6,i) = unconv(a+1,b-1,n1d)
        enddo
      else
        write(*,*) 'Invalid lattice type!'
        stop
      endif

      return
      end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      integer function unconv(x,y,n1d)
      implicit none

c takes the Cartesian coordinates of a site, applies periodic boundary
c conditions to them, and converts back to site number.

      integer x,y,n1d

      if (x.eq.0) x = n1d
      if (y.eq.0) y = n1d
      if (x.gt.n1d) x = x-n1d
      if (y.gt.n1d) y = y-n1d
        
      unconv = x + n1d*(y-1)

      return
      end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      real*8 FUNCTION random(init)
c   Generates random numbers  0 < r < 1 
c   multipl=524287 => Period=max=2**31-2 , seed=any int.

      include 'walk.cm'

      integer init
      real*8 scope, multipl, scopinv
      parameter(scope=2.0d0**31-1.0d0,scopinv=1.0/scope
     &    ,multipl=524287.0d0) 
      real*8 r_int
      SAVE r_int

      if (init.eq.0) r_int = seed

      r_int=dmod(r_int*multipl, scope)
      random=r_int*scopinv
      return
      end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      ! this subroutine removes all spaces from a string by shifting
      !   characters left
      ! Author: W. Magro
      ! Date:  July 11, 1995
      subroutine compactstr(string)
      character*(*) string

      ! locals
      integer i, j

      j = 1
      do i=1,len(string)
         if(string(i:i).ne." ") then
             if(j.ne.i) string(j:j) = string(i:i)
             j = j + 1
         endif
      enddo

      do i=j,len(string)
          string(i:i) = " "
      enddo

      return
      end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c openapp opens up a file at the end of the file, so that lines
c can be appended to it.  You need to input the lun and the filename.

      subroutine openapp(l,filename)
      implicit none
      
      integer i,l,length
      character*50 filename

      length = index(filename,' ')-1
      open (l,file=filename(1:length))
      do i=1,400000000
        read(l,*,END=99)
      enddo
 99   continue

      return
      end

