program two_particles_in_box_animate
      real x(2),y(2),x_old(2),y_old(2),vx(2),vy(2)
      real force_x(2),force_y(2),force_x_old(2),force_y_old(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,1,0)
      call pgsfs(2)

      dt = 0.0001
      dt2 = dt*dt

      do i_step =1, 10000000

c find forces from pairs
        n=2
        do i=1,n
          force_x(i)=0.0
          force_y(i)=0.0
        end do
        do i=1,n-1
          do j=i+1,n
            call lj_force(x(i),y(i),x(j),y(j),force_x_tmp,force_y_tmp)
            force_x(i)=force_x(i)+force_x_tmp
            force_y(i)=force_y(i)+force_y_tmp
            force_x(j)=force_x(j)-force_x_tmp
            force_y(j)=force_y(j)-force_y_tmp
          end do
        end do
c 


       do i=1,2

c Now move into subrutine wall_reflect START
c       if(x(i)>x_range) vx(i) = (-1)*vx(i)
c       if(x(i)<0.0) vx(i) = (-1)*vx(i)
c       if(y(i)>y_range) vy(i) = (-1)*vy(i)
c       if(y(i)<0.0) vy(i) = (-1)*vy(i)
c Now move into subrutine wall_reflect END

        call wall_reflect(x(i),y(i),vx(i),vy(i),x_range,y_range)

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

c       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(i) * dt2
        y(i) = y(i) + vy(i) * dt + 0.5 * force_y(i) * dt2
        vx(i) = vx(i) + (force_x(i)+force_x_old(i))/2 * dt
        vy(i) = vy(i) + (force_y(i)+force_y_old(i))/2 * 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)
        force_x_old(i) = force_x(i)
        force_y_old(i) = force_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


      subroutine wall_reflect(x,y,vx,vy,x_range,y_range)
        if(x>x_range) vx = (-1)*vx
        if(x<0.0) vx = (-1)*vx
        if(y>y_range) vy = (-1)*vy
        if(y<0.0) vy = (-1)*vy
      end subroutine