program two_particles_in_box_animate
      real x(2),y(2),vx(2),vy(2),x_old(2),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

      do i=1,2
       x(i) = 5.0
       y(i) = 5.0
      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)
      call pgpap(4.0,0.75)
      dt = 0.00001

      do j =1, 10000000
       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(i) = x(i) + vx(i) * dt
        y(i) = y(i) + vy(i) * 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