program main
      implicit none

      integer MAX_ELK
      parameter (MAX_ELK = 2000)

      real*8  x(2,MAX_ELK),v(2,MAX_ELK),f(2,MAX_ELK)
      integer n,nsteps,seed,noutput,nalive,step,i
      real*8  radius,vsum,distsq,radsq

c User inputs.

 10   write (6, *) 'Enter number of elk to simulate: '
      read (5, *) n
      if (n .gt. MAX_ELK) goto 10

      write (6, *) 'Enter radius of wilderness: '
      read (5, *) radius 

      write (6, *) 'Enter number of time steps: '
      read (5, *) nsteps

      write (6, *) 'Enter random number seed: '
      read (5, *) seed

 20   write (6, *) 'Enter which elk to print out: '
      read (5, *) noutput
      if (noutput .gt. n) goto 20

c Initialize the elk positions and velocities.

      call initialize_elk(n,x,v,radius,seed)

c Dump initial snapshot to file.

      call dump(n,x,0,radius)

c Timestep loop:
c  Compute force on each elk.
c  Move elk under influence of forces.
c  Dump snapshot of coords.

      do step = 1,nsteps

        call force_elk(n,x,f,radius)
        call move_elk(n,x,v,f,radius)
        call dump(n,x,step,radius)

      enddo

c Close snapshot file.

      call dump(n,x,-1,radius)

c Compute output.

      radsq = radius * radius
      nalive = 0
      vsum = 0.0

      do i = 1,n
        distsq = x(1,i)*x(1,i) + x(2,i)*x(2,i)
        if (distsq .le. radsq) then
          nalive = nalive + 1
          vsum = vsum + sqrt(v(1,i)*v(1,i) + v(2,i)*v(2,i))
        endif
      enddo

      if (nalive .gt. 0) vsum = vsum/nalive

c Print output.

      write (6,*) '# of Elk alive =',nalive
      write (6,*) 'Average speed of live elk =',vsum
      write (6,*) 'Coords of elk',noutput,' =',
     $     x(1,noutput),x(2,noutput)

      end


c Initialize the position and velocity of each elk.

      subroutine initialize_elk(n,x,v,radius,seed)
      implicit none

      integer   n
      real*8    x(2,n),v(2,n)
      real*8    radius
      integer   seed

      integer   i
      real*8    radsq,distsq,sigma
      real*8    rand1

      radsq = radius * radius
      sigma = sqrt(3.1415926/n) * radius

c generate random coordinates until elk is inside Gila
c generate random velocity

      do i = 1,n

 10     x(1,i) = (2.0 * rand1(seed) - 1.0) * radius
        x(2,i) = (2.0 * rand1(seed) - 1.0) * radius
        distsq = x(1,i)*x(1,i) + x(2,i)*x(2,i)
        if (distsq .gt. radsq) goto 10

        v(1,i) = (rand1(seed)-0.5) * sigma
        v(2,i) = (rand1(seed)-0.5) * sigma

      enddo

      return
      end


c Compute the forces acting on each elk.

      subroutine force_elk(n,x,f,radius)
      implicit none

      integer   n
      real*8    x(2,n),f(2,n)
      real*8    radius

      integer   i,j
      real*8    epsilon,sigma,sigmasq,radsq,distsq
      real*8    delx,dely,sig2r2,force

      epsilon = 1.0
      sigma = sqrt(3.1415926/n) * radius
      sigmasq = sigma * sigma
      radsq = radius * radius

c Zero out forces.

      do i = 1,n
        f(1,i) = 0.0
        f(2,i) = 0.0
      enddo

c N^2 double loop to accumulate forces.
c Elk outside of Gila do not contribute.

      do i = 1,n
        distsq = x(1,i)*x(1,i) + x(2,i)*x(2,i)
        if (distsq .le. radsq) then
          do j = 1,n
            distsq = x(1,j)*x(1,j) + x(2,j)*x(2,j)
            if (distsq .le. radsq .and. i .ne. j) then
              delx = x(1,i) - x(1,j)
              dely = x(2,i) - x(2,j)
              distsq = delx*delx + dely*dely
              if (distsq .le. sigmasq) then
                force = 2.0*epsilon/sigmasq
              else
                sig2r2 = sigmasq / distsq
                force = 4.0*epsilon/distsq * sig2r2 * (sig2r2-0.5)
              endif
              f(1,i) = f(1,i) + delx*force
              f(2,i) = f(2,i) + dely*force
            endif
          enddo
        endif
      enddo

      return
      end


c Move the elk using Newtonian mechanics.

      subroutine move_elk(n,x,v,f,radius)
      implicit none

      integer   n
      real*8    x(2,n),v(2,n),f(2,n)
      real*8    radius

      integer   i
      real*8    radsq,distsq,dt

c Simple 1st-order Euler integrator.

      dt = 0.5

      do i = 1,n
        x(1,i) = x(1,i) + dt*v(1,i) + 0.5*dt*dt*f(1,i)
        x(2,i) = x(2,i) + dt*v(2,i) + 0.5*dt*dt*f(2,i)
        v(1,i) = v(1,i) + dt*f(1,i)
        v(2,i) = v(2,i) + dt*f(2,i)
      enddo

c Set velocity to 0.0 for elk that wander onto rancher's land.

      radsq = radius * radius

      do i = 1,n
        distsq = x(1,i)*x(1,i) + x(2,i)*x(2,i)
        if (distsq .gt. radsq) then
          v(1,i) = 0.0
          v(2,i) = 0.0
        endif
      enddo
      
      return
      end


c Dump snapshot of coords to file.

      subroutine dump(n,x,step,radius)
      implicit none

      integer   n,step
      real*8    x(2,n)
      real*8    radius

      integer   i

      if (step .eq. -1) then
        close (1)
        return
      endif

      if (step .eq. 0) then
        open (1,file='dump.out')
        close (1,status='delete')
        open (1,file='dump.out')
        write (1,*) 'ITEM: NUMBER OF ATOMS'
        write (1,*) n
        write (1,*) 'ITEM: BOX BOUNDS'
        write (1,*) -2.0*radius,2.0*radius
        write (1,*) -2.0*radius,2.0*radius
      endif

      write (1,*) 'ITEM: TIMESTEP'
      write (1,*) step
      write (1,*) 'ITEM: ATOMS'
      do i = 1,n
        write (1,*) i,1,x(1,i),x(2,i)
      enddo
      
      return
      end


C Park/Miller RNG
 
      double precision function rand1(iseed)
      double precision aa,mm,sseed
      parameter (aa=16807.0D0,mm=2147483647.0D0)
      
      sseed = iseed
      sseed = mod(aa*sseed,mm)
      iseed = sseed
      rand1 = sseed/mm

      return
      end