program shoot_target_water
      real x_draw(2), y_draw(2)
      real x_water(2), y_water(2)
      real x_wall(2), y_wall(2)
      real k
      logical no_wall

c plot related setting
      i_dummy=pgopen('/xwin')
      call pgenv(0.0,10.0,0.0,10.0,1,-1)
      call pgsfs(2)

c plot target
      y_target=2.5
      x_target_min=7
      x_target_max=8
      y_draw(1)=y_target
      y_draw(2)=y_target
      x_draw(1)=x_target_min
      x_draw(2)=x_target_max
      call pgsci(2)      
      call pgline(2,x_draw,y_draw)
      x_water(1)=5
      y_water(1)=5
      x_water(2)=10
      y_water(2)=y_water(1)
      call pgsci(4)
      call pgline(2,x_water,y_water)
      x_wall(1)=x_water(1)
      y_wall(1)=y_water(1)+0.5
      x_wall(2)=x_water(1)
      y_wall(2)=0
      call pgsci(1)
      call pgline(2,x_wall,y_wall)

c Input parameters
1003  write(*,*)'Vx=?'
      read(*,*)vx
      write(*,*)'Vy=?'
      read(*,*)vy

c cannon fly

      k=0.0
      m=1.0
      no_wall=.true.

      ax=0
      ay=-9.8 
      y=1
      x=0+0.12
      dt=0.00001
      xold=x
      yold=y

      do while (y>0)

c check wall hitting
        if((x>x_wall(1)-0.13).and.(x<x_wall(1)).and.(y<y_wall(1)))
     *  then
         no_wall=.false.
         vx=(-1)*vx
        endif
        if((vx<0).and.(x<0.13))vx=(-1)*vx

        v=sqrt(vx**2+vy**2)
        if ((x>x_water(1)).and.(y<y_water(1)))then
          k=2.0
         else
          k=0
        endif
        f_drag=(-1)*k*v**2
        f_dragx=f_drag*vx/v
        f_dragy=f_drag*vy/v
        ax=f_dragx/m
        ay=-9.8+f_dragy/m
        vx=vx+ax*dt
        vy=vy+ay*dt
        x=x+vx*dt
        y=y+vy*dt

        call pgbbuf
        call pgsci(0)
        call pgcirc(xold,yold,0.1) 
        call pgsci(1)
        call pgcirc(x,y,0.1)
        call pgsci(4)
        call pgline(2,x_water,y_water)
        call pgebuf

c check_hit
        if((y<y_target).and.(y>(y_target-0.1)))then
          if((x>x_target_min).and.(x<x_target_max)) goto 1000
        endif

        xold=x
        yold=y
      end do

      goto 1001
1000  call wining_animation(y_target,x_target_min,x_target_max)
1001  continue
      write(*,*)'Play again?(1=Yes,0=No)'
      read(*,*)i_again
      if(i_again==1) goto 1003
      end

      subroutine wining_animation(y_target,x_target_min,x_target_max)
      x_middle=(x_target_min + x_target_max)/2
      d_theta=3.14/1800
      theta=0
      r0=1
      rold=r
      do i=1,20000
        theta=theta+d_theta
        r=r0+0.5*sin(theta)
        call pgbbuf
        call pgsci(0)
        call pgcirc(x_middle,y_target,rold)
        call pgsci(3)
        call pgcirc(x_middle,y_target,r)
        call pgebuf
        rold=r
      end do
      end