program stright_motion

      write(*,*)'Initial x ?'
      read(*,*)x
      write(*,*)'Initial y ?'
      read(*,*)y
      write(*,*)'Initial Vx ?'
      read(*,*)vx
      write(*,*)'Initial Vy ?'
      read(*,*)vy

      x_min = 0.0
      y_min = 0.0
      x_max = 10.0
      y_max = 10.0

      if(pgopen('/xwin')<=0)stop
      call pgenv(x_min,x_max,y_min,y_max,1,0)

      call pgsfs(0)
      call draw_wall()

      dt = 0.0001      
      do i=1,10000000
        x = x + vx*dt
        y = y + vy*dt
        call box(x,y,vx,vy,x_min,y_min,x_max,y_max)
c       call pbc(x,y,x_min,y_min,x_max,y_max)
        call wall(x,y,xold,yold,vx,vy)
        call pgbbuf
        call pgsci(0)
        call pgcirc(xold,yold,0.1)
        call pgsci(1)
        call pgcirc(x,y,0.1)
        call pgebuf
        xold=x
        yold=y
      end do  

      end

      subroutine box(x,y,vx,vy,x_min,y_min,x_max,y_max)
      if(x<x_min)vx=(-1)*vx
      if(y<y_min)vy=(-1)*vy
      if(x>x_max)vx=(-1)*vx
      if(y>y_max)vy=(-1)*vy
      end

      subroutine pbc(x,y,x_min,y_min,x_max,y_max)
      if(x<x_min)x=x+(x_max-x_min)
      if(y<y_min)y=y+(y_max-y_min)
      if(x>x_max)x=x-(x_max-x_min)
      if(y>y_max)y=y-(y_max-y_min)
      end

      subroutine draw_wall()
      real pt_x(4),pt_y(4)
      pt_x(1)=5.0
      pt_y(1)=5.0
      pt_x(2)=5.5
      pt_y(2)=5.0
      pt_x(3)=5.5
      pt_y(3)=0.0
      pt_x(4)=5.0
      pt_y(4)=0.0
      call pgpoly(4,pt_x,pt_y) 
      end

      subroutine wall(x,y,xold,yold,vx,vy)
c     real pt_x(4),pt_y(4)
      w_left=5.0
      w_right=5.5
      w_bottom=0.0
      w_top=5.0
c     pt_x(1)=5.0
c     pt_y(1)=5.0
c     pt_x(2)=5.5
c     pt_y(2)=5.0
c     pt_x(3)=5.5
c     pt_y(3)=0.0
c     pt_x(4)=5.0
c     pt_y(4)=0.0
c     call pgpoly(4,pt_x,pt_y) 
      if((y<w_top).and.(y>w_bottom))then
        if((xold<w_left).and.(x>w_left))vx = (-1)*vx
        if((xold>w_right).and.(x<w_right))vx = (-1)*vx
      endif
      if((x<w_right).and.(x>w_left))then 
        if((yold>w_top).and.(y<w_top))vy = (-1)*vy
        if((yold<w_bottom).and.(y>w_bottom))vy = (-1)*vy
      endif
      end