program pinball
      implicit none
      integer pgopen,i_again,n_nail,ii,num_of_nails_left
      parameter (n_nail = 150)
      real x,vx,ax, y,vy,ay, vx_ini,vy_ini, x_range,y_range,delta_t
      real mass,k,radius_2
      real nail_x(n_nail),nail_y(n_nail),old_pos_x,old_pos_y,radius
      parameter (radius = 0.5)

      radius_2 = radius - 0.02
      call input_n_init(vx_ini,vy_ini,delta_t,mass,k)

      if (pgopen('/xwin').le.0) stop
      call pgpap(4.0, 2.0)
      call pgsfs(2)

 100  write(*,*) 'How far (x-dir) and how high (y-dir) is the plot ?'
      read (*,*) x_range, y_range
      call pgeras
      call pgenv(0.0,x_range,0.0,y_range,1,0)

      call nail_line_up (x_range,y_range,n_nail,nail_x,nail_y,
     &     num_of_nails_left,radius)

      x = 0.0+radius
      y = 0.9*y_range
      vx = vx_ini
      vy = vy_ini
      old_pos_x = x
      old_pos_y = y

      do while (y.gt.0.0)
        call fly_drift (x,vx,ax,y,vy,ay,delta_t,mass,k)
        call pgbbuf
        call pgsci (0)
        call pgcirc (old_pos_x,old_pos_y,radius_2)
        call pgsci (1)
        call pgcirc (x,y,radius_2)
        call pgpt (n_nail-num_of_nails_left,nail_x,nail_y,-1)
        call pgebuf
        old_pos_x = x
        old_pos_y = y
        call test_collide(x,vx,y,vy,nail_x,nail_y,n_nail,radius,ii,
     &       x_range,y_range)
        if (ii.ne.0) call bounch (x,vx,y,vy,nail_x(ii),nail_y(ii))
      end do

 102  write(*,*) 'Again ? (1=Yes;0=No)'
      read(*,*) i_again
      if (i_again.eq.1) goto 100
      if (i_again.eq.0) stop
      goto 102

      call pgclos

      end



      subroutine input_n_init(vx,vy,delta_t,mass,k)
      implicit none
      real vx,vy,v,angle,ke,delta_t,mass,k
      integer i_gun_powder
      write(*,*) 'How maqny packs of gun powder (integer) ?'
      read (*,*) i_gun_powder
      write (*,*) 'What is the angle of shooting ?'
      read (*,*) angle
      write(*,*) 'What is the small time-interval delta_t ?'
      read (*,*) delta_t
      write(*,*) 'What is the mass of the particle ?'
      read (*,*) mass
      write(*,*) 'What is the coefficient of air resistence ?'
      read (*,*) k
      angle = angle * 3.14159268/180
      ke = 1.0*i_gun_powder
      v = sqrt(2.0*ke/mass)
      vx = v*cos(angle)
      vy = v*sin(angle)
      end



      subroutine fly_drift(x,vx,ax,y,vy,ay,delta_t,mass,k)
      implicit none
      real x,vx,ax,y,vy,ay,delta_t,theda,g
      real v_old,v_old_x,v_new_x,v_old_y,v_new_y,mass,v_mid_x
      real a_mid_x,a_mid_y,a_drift_x,a_drift_y,a_drift,k,v_mid_y

      g = -9.8

      v_old_x = vx
      v_old_y = vy

      v_old = sqrt(v_old_x**2 + v_old_y**2)
      a_drift = -1*k*v_old**2/mass
c     theda = atan(v_old_y/v_old_x)
c     theda = asin(v_old_y/v_old)
c     a_drift_x = a_drift*cos(theda)
c     a_drift_y = a_drift*sin(theda)
      a_drift_x = a_drift*(v_old_x/v_old)
      a_drift_y = a_drift*(v_old_y/v_old)
      a_mid_x = a_drift_x
      a_mid_y = a_drift_y + g
      v_mid_x = v_old_x + a_mid_x*0.5*delta_t
      v_mid_y = v_old_y + a_mid_y*0.5*delta_t

      v_new_x = v_old_x + a_mid_x*delta_t
      x = x + v_mid_x*delta_t

      v_new_y = v_old_y + a_mid_y*delta_t
      y = y + v_mid_y*delta_t

      vx = v_new_x
      vy = v_new_y

      end


      subroutine test_collide(x,vx,y,vy,nail_x,nail_y,n_nail,radius,ii
     &,x_range,y_range)
      implicit none
      integer n_nail,i,ii
      real nail_x(n_nail),nail_y(n_nail),x,vx,y,vy,radius
      real d(n_nail),x_range,y_range
c assuming hitting to one of the nails only
      ii = 0
      do i=1,n_nail
        d(i) = sqrt( (x-nail_x(i))**2 + (y-nail_y(i))**2 )
        if (d(i).le.radius) ii = i
      end do
      if ((x.le.radius).or.(x.ge.x_range-radius)) vx = -1*vx
      if ((y.le.radius).or.(y.ge.y_range-radius)) vy = -1*vy
      end

      subroutine bounch (x,vx,y,vy,nail_x,nail_y)
      real x,vx,y,vy,nail_x,nail_y
      real normal_x,normal_y,v_para_x,v_para_y,v_perp_x,v_perp_y,temp
      real v_dot_n
c N_vec = nail_vec - ball_vec, then normalise N_vec to become n_vec
      normal_x = nail_x - x
      normal_y = nail_y - y
      temp = sqrt(normal_x**2 + normal_y**2)
      normal_x = normal_x / temp
      normal_y = normal_y / temp
c inner product of V_ball and n_vec
      v_dot_n = vx*normal_x + vy*normal_y
c V_perpendicular = Vdotn n
      v_perp_x = v_dot_n*normal_x
      v_perp_y = v_dot_n*normal_y
c V_parallele = V_ball - V_peprp
      v_para_x = vx - v_perp_x
      v_para_y = vy - v_perp_y
c bouch back on n diection
      v_perp_x = -1*v_perp_x
      v_perp_y = -1*v_perp_y
c New vx anc vy after bouch back
      vx = v_para_x + v_perp_x
      vy = v_para_y + v_perp_y
      end

      subroutine nail_line_up (x_range,y_range,n_nail,nail_x,nail_y,
     &           num_of_nails_left,radius)
      implicit none
      integer n_nail
      real x_range, y_range, nail_x(n_nail), nail_y(n_nail),shift
      integer n_col, n_row, row_size,col_size,num_of_nails_left,i
      integer ix,iy
      real dx,dy,h_space,v_space,radius
 200  write(*,*) "What is the horizontal and vertical nail spacing with 
     &respect to ball diameter ?"
      read (*,*) h_space, v_space
c     if (h_space.le.1.0)) stop ('Nails too close !')
c     if (h_space.le.1.0) then
c       write(*,*) 'Nails too close !'
c       goto 200
c     endif
      h_space = h_space * 2 * radius
      v_space = v_space * 2 * radius
      row_size = int(x_range/h_space) - 1 
      col_size = int(0.75*y_range/v_space) - 1
c     col_size = int(y_range/v_space) - 1
      if ((row_size.eq.0).or.(col_size.eq.0)) stop ('table too small')
      dx = x_range / row_size
      dy = 0.75*y_range / col_size
c     dy = y_range / col_size
      num_of_nails_left=n_nail
c     shift = 1.0
      shift = 2.0*radius
      do iy = 1, col_size
        do ix = 1, row_size
          i = (iy-1)*row_size + ix
          nail_x(i) = shift + (ix-1)*dx
          nail_y(i) = iy* dy
c         call pgpt (1,nail_x,nail_y,-1)
          num_of_nails_left = num_of_nails_left - 1
          if (num_of_nails_left.eq.0) exit 
        end do
        if (num_of_nails_left.eq.0) exit 
c       shift = 1.0 - shift
        shift = 2.0*radius - shift
      end do
      call pgpt (n_nail-num_of_nails_left,nail_x,nail_y,-1)
      end