c  test of runge-kutta  method
      real t,x(4),dt,t0,x0(4)
      real xreal,error
      real s,s2
      real tlim
      real pi
      real mass
      integer n,nwrite, i
      common /para/alpha,beta
      common /velocity/vabs
c
      open(unit=2,file='ball.dat')

      pi=3.1415926535897932385
      write(6,*) ' v0 (km/h)?'
      read(5,*) v0
      write(6,*) ' theta ?'
      read(5,*) theta

      rho=1.2
      rball=0.0715/2.0
      mass=0.1445
      eta=1.8e-5
      v0=v0*1000.0/3600.0
      alpha=6.0*pi*rball*eta/mass
      beta=0.25*pi*rho*rball**2/mass

      t0=0.0                                                            
      x0(1)=0.0
      x0(2)=0.0
      x0(3)=v0*cos(pi*theta/180)
      x0(4)=v0*sin(pi*theta/180)
      write(6,*) ' dt ?'
      read(5,*) dt
      i=0
      t=t0                                                              
      x(1)=x0(1)                        
      x(2)=x0(2)                        
      x(3)=x0(3)                        
      x(4)=x0(4)                        
   10 continue
      i=i+1
      vabs=sqrt(x(3)**2+x(4)**2)
      call runge2(4,t,x,dt)                                               
      write(2,601) t,x(1),x(2),x(3),x(4)                                     
      write(6,601) t,x(1),x(2),x(3),x(4)                                     
      if(x(2).lt.0.0) then
        stop
      endif
      go to 10                                                          
  601 format(1h ,5(1pe12.5,2x))                                         
      end                                                               
*********************************************************************** 
      real function f(i,t,x)                                                   
      integer i
      real t
      real x(4)
      real pi
      common /para/alpha,beta
      common /velocity/vabs
      pi=3.1415926535897932385
      if(i.eq.1) then	
        f=x(3)                                          
      elseif(i.eq.2) then
        f=x(4)
      elseif(i.eq.3) then
        f=0.0 -alpha*x(3) -beta*x(3)*vabs
      elseif(i.eq.4) then
        f=-9.8 -alpha*x(4) -beta*x(4)*vabs
      endif
      return                                                            
      end                                                               
*********************************************************************** 
c  subroutine of runge-kutta method of many valuables
      subroutine runge2(n,x,y,dx)                                       
      real h1,k(10,4),k1(10),y1(10),y2(10),y3(10),y(10)
c
      h1=dx/2.0                                                         
      do 10 i1=1,n
        k(i1,1)=f(i1,x,y)*dx                                              
        y1(i1)=y(i1)+k(i1,1)/2.0
  10  continue                                                          
      do 20 i2=1,n
        k(i2,2)=f(i2,x+h1,y1)*dx                                          
        y2(i2)=y(i2)+k(i2,2)/2.0
  20  continue                                                          
      do 30 i3=1,n
        k(i3,3)=f(i3,x+h1,y2)*dx                                          
        y3(i3)=y(i3)+k(i3,3)
  30  continue                                                          
      do 40 i4=1,n
        k(i4,4)=f(i4,x+dx,y3)*dx                                          
   40 continue                                      
c
      x=x+dx                                                            
      do 50 i5=1,n
        k1(i5)=(k(i5,1)+2.0*k(i5,2)+2.0*k(i5,3)+k(i5,4))/6.0              
        y(i5)=y(i5)+k1(i5)
  50  continue                                                          
      return
      end