c
c Simulation of trajectory for a curve/knuckle ball or a line drive hit
c        - Euler or Runge-Kutta (2nd order usual one)
c Fortran version written by H. Nakanishi, need to be compiled with "-lpepl".
c
	program ball
c
c       Declare the arrays we will need
c
	dimension x(5003),y(5003),z(5003)
c
c       Use subroutines to do the work
c
	call initialize(dt,vinit,om,th0,mp,bm,sm,lsym,nsym,method)
	call calculate(x,y,z,dt,vinit,om,th0,mp,bm,sm,n,method)   
	call display(x,z,n,om,th0,mp,bm,sm,dt,lsym,nsym,method)
	stop
	end
c
	subroutine initialize(dt,vinit,om,th0,mp,bm,sm,lsym,nsym,method)
c
c       Initialize variables
c
	character ans,yes
	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), Knuckle (2), or line drive hit (3) ?'
	read(5,*) mp
	if(mp.eq.1) then
		print *,'magnus/m -> '
		read(5,*) sm
		th0=0
	else
		print *,'theta0 (radians) -> '
		read(5,*) th0
		sm=0
	endif
	print *,'set line, symbol?'
	read(5,14) ans
14	format(a1)
	if(ans.eq.yes) then
		print *,'lsym, nsym ? -> '
		read(5,*) lsym,nsym
	else
		lsym=-1
		nsym=1
	endif
	return
	end
c
	subroutine calculate(x,y,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),y(1),z(1)
	x(1)=0
	y(1)=1.2
	z(1)=0
	th=th0
	vx=vinit
	vy=0
	if(mp.eq.3) then
		vx=vinit*cos(th0)
		vy=vinit*sin(th0)
	endif
	vz=0
	nmax=5000
	if(method.eq.1) then
		do 10 i = 2,nmax
		call dv(x(i-1),y(i-1),z(i-1),th,om,vx,vy,vz,mp,bm,sm
     c                 ,dx,dy,dz,dvx,dvy,dvz)
		x(i)=x(i-1)+dt*dx
		y(i)=y(i-1)+dt*dy
		z(i)=z(i-1)+dt*dz
		th=th+dt*om*6.2831853
		vx=vx+dt*dvx
		vy=vy+dt*dvy
		vz=vz+dt*dvz
		if((mp.eq.1.or.mp.eq.2).and.x(i).ge.18.29) then
			n=i
			go to 15
		elseif((mp.eq.3).and.y(i).lt.0) then
			n=i
			go to 16
		endif
10		continue
	elseif(method.eq.2) then
		do 30 i = 2,nmax
		call dv(x(i-1),y(i-1),z(i-1),th,om,vx,vy,vz,mp,bm,sm
     c                 ,dx,dy,dz,dvx,dvy,dvz)
		x1=x(i-1)+0.5*dt*dx
		y1=y(i-1)+0.5*dt*dy
		z1=z(i-1)+0.5*dt*dz
		th1=th+0.5*dt*om*6.2831853
		vx1=vx+0.5*dt*dvx
		vy1=vy+0.5*dt*dvy
		vz1=vz+0.5*dt*dvz
		call dv(x1,y1,z1,th1,om,vx1,vy1,vz1,mp,bm,sm
     c                 ,dx2,dy2,dz2,dvx2,dvy2,dvz2)
		x(i)=x(i-1)+dt*dx2
		y(i)=y(i-1)+dt*dy2
		z(i)=z(i-1)+dt*dz2
		th=th+dt*om*6.2831853
		vx=vx+dt*dvx2
		vy=vy+dt*dvy2
		vz=vz+dt*dvz2
		if((mp.eq.1.or.mp.eq.2).and.x(i).ge.18.29) then
			n=i
			go to 15
		elseif((mp.eq.3).and.y(i).lt.0) then
			n=i
			go to 16
		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))
  	y(n)=(y(n-1)*x(n)-x(n-1)*y(n)+18.29*(y(n)-y(n-1)))/(x(n)-x(n-1))
	x(n)=18.29
	return
16	a=-y(n)/y(n-1)
	x(n)=(x(n)+a*x(n-1))/(1+a)
	y(n)=0
	z(n)=(z(n)+a*z(n-1))/(1+a)
	return
	end
c
	subroutine dv(x0,y0,z0,th0,om,vx0,vy0,vz0,mp,bm,sm
     c  ,dx,dy,dz,dvx,dvy,dvz)
	g=9.8
	vr=28.9
	dx=vx0
	dy=vy0
	dz=vz0
	v0=sqrt(vx0**2+vy0**2+vz0**2)
	dvx=-bm*v0*vx0
	dvy=-g-bm*v0*vy0
	if(mp.eq.1) then
		dvz=-sm*vx0*om*6.2831853-bm*v0*vz0
	else
		dvz=((v0/vr)**2)*4.9*(sin(4*th0)-0.25*sin(8*th0)
     c		+0.08*sin(th0)-0.025*sin(16*th0))-bm*v0*vz0
	endif
	return
	end
c
	subroutine display(x,y,nmax,om,th0,mp,bm,sm,dt,lsym,nsym,method)
c
c       First set up title and label axes for graph. Plotting is a lot
c       of work in fortran.
c       This version displays output as well as writes to a file "graph.out".
c
	dimension x(1),y(1)
	call usrmon(.true.,.false.,-0.5,9.,-0.5,12.)
	call pltlun(19,.true.,.false.)
	call pltlfn('graph.out')
	call plots
	call plot(0.,3.5,-3)
	if(method.eq.1) then
		if(mp.eq.1.or.mp.eq.2) then
	call text(0.2,5.8,0.2,'Curve/Knuckle ball: Euler',0.,25,0)
		else
	call text(0.2,5.8,0.2,'Line drive hit: Euler',0.,21,0)
		endif
	elseif(method.eq.2) then
		if(mp.eq.1.or.mp.eq.2) then
	call text(0.2,5.8,0.2,'Curve/Knuckle ball: Runge-Kutta2',0.,32,0)
		else
	call text(0.2,5.8,0.2,'Line drive hit: Runge-Kutta2',0.,28,0)
		endif
	endif
	call text(0.2,5.5,0.2,'Deflection vs Distance',0.,22,0)
	call scalex(x,5.,nmax,1)
	call axctl(0.15,0.04,0.15,0.2,0.2,-1)
	call axisx(0.,0.,'Distance (m)',-12,5.,0.,
     c	x(nmax+1),x(nmax+2),x(nmax+3),4)
	call axisx(0.,5.,' ',1,5.,0.,
     c	x(nmax+1),x(nmax+2),x(nmax+3),20)
	call scalex(y,5.,nmax,1)
	call axctl(0.15,0.04,0.15,0.2,0.2,2)
	call axisx(0.,0.,'Deflection (m)',14,5.,90.,
     c	y(nmax+1),y(nmax+2),y(nmax+3),-5)
	call axisx(5.,0.,' ',-1,5.,90.,
     c	y(nmax+1),y(nmax+2),y(nmax+3),-21)
	call line(x,y,nmax,1,lsym,nsym)
	call number(0.5,3.3,0.2,om,0.,'''omega= '',f5.2')
	call number(0.5,3.6,0.2,bm,0.,'''bm= '',1pe8.2')
	call number(0.5,3.9,0.2,dt,0.,'''dt = '',f6.3')
	if(mp.eq.1) then
		call number(0.5,4.2,0.2,sm,0.,'''sm= '',1pe8.2')
	else
		call number(0.5,4.2,0.2,th0,0.,'''th0 = '',f5.2')
	endif
	call plot(0.,0.,999)
	return
	end

