program two_particles_in_box_animate
      real x(2),y(2),x_old(2),y_old(2),vx(2),vy(2)
      integer pgopen
      do i=1,2
        write(*,*)'What is initial vx of particle',i,'?'
        read(*,*) vx(i)
        write(*,*)'What is initial vy of particle',i,'?'
        read(*,*) vy(i)
      end do

      x_range = 10.0
      y_range = 10.0

c     do i=1,2
       x(1) = 5.0
       y(1) = 5.0
       x(2) = 0.0
       y(2) = 0.0
c     end do

      if( pgopen('/xwin') <= 0 ) stop
c     call pgenv(0.0,x_range,0.0,y_range,0,-1)
      call pgenv(-0.15,x_range+0.15,-0.15,y_range+0.15,0,-1)
      call pgsfs(2)

      dt = 0.0001
      dt2 = dt*dt

      do j =1, 10000000
c      do while(
       do i=1,2
        if(x(i)>x_range) vx(i) = (-1)*vx(i)
        if(x(i)<0.0) vx(i) = (-1)*vx(i)
        if(y(i)>y_range) vy(i) = (-1)*vy(i)
        if(y(i)<0.0) vy(i) = (-1)*vy(i)

          x_self=x(i)
          y_self=y(i)
        if(i==1)then
          x_other=x(i+1)
          y_other=y(i+1)
         else
          x_other=x(i-1)
          y_other=y(i-1)
        endif

        call lj_force(x_self,y_self,x_other,y_other,force_x,force_y)

        x(i) = x(i) + vx(i) * dt + 0.5 * force_x * dt2
        y(i) = y(i) + vy(i) * dt + 0.5 * force_y * dt2
        vx(i) = vx(i) + force_x * dt
        vy(i) = vy(i) + force_y * dt

        call pgbbuf
        call pgsci(0)
        call pgcirc(x_old(i),y_old(i),0.1)
        call pgsci(1)
        call pgcirc(x(i),y(i),0.1)
        call pgebuf

        x_old(i) = x(i)
        y_old(i) = y(i)
       enddo
      enddo
      call pgclos
      end

      subroutine lj_force(x1,y1,x2,y2,force_x,force_y)
      sigma = 1.0
      epsilon = 1.0

      r = sqrt( (x2-x1)**2 + (y2-y1)**2 )
      s_o_r = sigma / r

      force = (-1)*24*epsilon*(2*s_o_r**12-s_o_r**6)/r

      force_x = force*(x2-x1)/r
      force_y = force*(y2-y1)/r
      end subroutine