c
c Simulation of trajectory for a curve or knuckle ball
c        - Euler or Runge-Kutta (2nd order usual one)
c Fortran version written by H. Nakanishi, need to be compiled with "-lpepl".
c
	program pitch
c
c       Declare the arrays we will need
c
	dimension x(5003), z(5003)
	character*12 out
c
c       Use subroutines to do the work
c
	call initialize(dt,vinit,om,th0,mp,bm,sm,out,method)
	call calculate(x,z,dt,vinit,om,th0,mp,bm,sm,n,method)   
	call save(x,z,n,om,th0,mp,bm,sm,dt,out,method)
	stop
	end
c
	subroutine initialize(dt,vinit,om,th0,mp,bm,sm,out,method)
	character ans,yes
	character*12 out
c
c       Initialize variables
c
	yes='y'
	print *,'Euler (1) or Runge-Kutta 2nd order (2) ? -> '
	read(5,*) method
	if(method.ne.1.and.method.ne.2) then
		print *,'must select 1 or 2 ..'
		stop
	endif
	print *,'initial velocity, time step, rev/s, drag/m -> '
	read(5,*) vinit,dt,om,bm
	print *,'Curve (1) or Knuckle (2) ?'
	read(5,*) mp
	if(mp.eq.1) then
		print *,'magnus/m -> '
		read(5,*) sm
		th0=0
	else
		print *,'theta0 -> '
		read(5,*) th0
		sm=0
	endif
	print *,'set output data file name ?'
	read(5,14) ans
14	format(a1)
	if(ans.eq.yes) then
		print *,'name ? -> '
		read(5,15) out
15		format(a12)
	else
		out = 'pitch.data'
	endif
	return
	end
c
	subroutine calculate(x,z,dt,vinit,om,th0,mp,bm,sm,n,method)   
c
c       Now use the Euler method or the Runge-Kutta (2nd order)
c
	dimension x(1),z(1)
	x(1)=0
	z(1)=0
	th=th0
	vx=vinit
	vz=0
	nmax=5000
	if(method.eq.1) then
		do 10 i = 2,nmax
		call dv(x(i-1),z(i-1),th,om,vx,vz,mp,bm,sm,dx,dz,dvx,dvz)
		x(i)=x(i-1)+dt*dx
		z(i)=z(i-1)+dt*dz
		th=th+dt*om*6.2831853
		vx=vx+dt*dvx
		vz=vz+dt*dvz
		if(x(i).ge.18.29) then
			n=i
			go to 15
		endif
10		continue
	elseif(method.eq.2) then
		do 30 i = 2,nmax
		call dv(x(i-1),z(i-1),th,om,vx,vz,mp,bm,sm,dx,dz,dvx,dvz)
		x1=x(i-1)+0.5*dt*dx
		z1=z(i-1)+0.5*dt*dz
		th1=th+0.5*dt*om*6.2831853
		vx1=vx+0.5*dt*dvx
		vz1=vz+0.5*dt*dvz
		call dv(x1,z1,th1,om,vx1,vz1,mp,bm,sm,dx2,dz2,dvx2,dvz2)
		x(i)=x(i-1)+dt*dx2
		z(i)=z(i-1)+dt*dz2
		th=th+dt*om*6.2831853
		vx=vx+dt*dvx2
		vz=vz+dt*dvz2
		if(x(i).ge.18.29) then
			n=i
			go to 15
		endif
30		continue
	endif
	n=nmax
15	z(n)=(z(n-1)*x(n)-x(n-1)*z(n)+18.29*(z(n)-z(n-1)))/(x(n)-x(n-1))
	x(n)=18.29
	return
	end
c
	subroutine dv(x0,z0,th0,om,vx0,vz0,mp,bm,sm,dx,dz,dvx,dvz)
	dx=vx0
	dz=vz0
	dvx=-bm*sqrt(vx0**2+vy0**2)*vx0
	if(mp.eq.1) then
		dvz=-sm*vx0*om*6.2831853
	else
		dvz=4.9*(sin(4*th0)-0.25*sin(8*th0)
     c		+0.08*sin(th0)-0.025*sin(16*th0))
	endif
	return
	end
c
	subroutine save(x,y,nmax,om,th0,mp,bm,sm,dt,out,method)
	character*12 out
	dimension x(1),y(1)
	open(1,file=out)
	rewind(1)
	write(1,10) (x(i),y(i),i=1,nmax)
10	format(1x,1p,e12.5,2x,e12.5)
	return
	end

