c
c Simulation of velocity vs. time for a large cannon
c        - Euler or 2nd order Runge-Kutta
c Program to accompany "Computational Physics" by N. Giordano and H. Nakanishi
c Fortran version written by H. Nakanishi, need to be compiled with "-lpepl".
c
	program cannon2 
c
c       Declare the arrays we will need
c
	dimension x(5003), y(5003)
c
c       Use subroutines to do the work
c
	call initialize(dt,vinit,theta,Am,lsym,nsym,method)   
	call calculate(x,y,dt,vinit,theta,Am,n,lsym,nsym,method)   
	call display(x,y,n,theta,Am,dt,lsym,nsym,method)
	stop
	end
c
	subroutine initialize(dt,vinit,theta,Am,
     c		lsym,nsym,method)
c
c       Initialize variables
c
	character ans,yes
	yes='y'
	print *,'Euler (1) or Runge-Kutta (2)? -> '
	read(5,*) method
	if(method.ne.1.and.method.ne.2) then
		print *,'must select 1 or 2 ..'
		stop
	endif
	print *,'initial velocity -> '
	read(5,*) vinit
	print *,'time step -> '
	read(5,*) dt 
	print *,'drag/m -> '
	read(5,*) Am
	print *,'firing angle -> '
	read(5,*) theta
	print *,'set line, symbol?'
	read(5,14) ans
14	format(a1)
	if(ans.eq.yes) then
		print *,'line and symbol numbers -> '
		read(5,*) lsym,nsym
	else
		lsym=-1
		nsym=1
	endif
	return
	end
c
	subroutine calculate(x,y,dt,vinit,theta,Am,n,
     c		lsym,nsym,method)
c
c       Now use the Euler method or the Runge-Kutta (2nd order)
c
	dimension x(1),y(1)
	x(1)=0
	y(1)=0
	vx=vinit*cos(3.141592*theta/180)
	vy=vinit*sin(3.141592*theta/180)
	nmax=5000
	if(method.eq.1) then
		do 10 i = 2,nmax
		x(i)=x(i-1)+vx*dt
		y(i)=y(i-1)+vy*dt
		f=Am*sqrt(vx**2+vy**2)
		vx=vx-f*vx*dt
		vy=vy-9.8*dt-f*vy*dt
		if(y(i).le.0) then
			n=i
			go to 15
		endif
10		continue
	else
		do 30 i = 2,nmax
		call mid(x1,x(i-1),y1,y(i-1),vx1,vx,
     c			vy1,vy,dt,Am)
		x(i)=x(i-1)+vx1*dt
		y(i)=y(i-1)+vy1*dt
		f=Am*sqrt(vx1**2+vy1**2)
		vx=vx-f*vx1*dt
		vy=vy-9.8*dt-f*vy1*dt
		if(y(i).le.0) then
			n=i
			go to 15
		endif
30		continue
	endif
	n=nmax
15	a=-y(n)/y(n-1)
	x(n)=(x(n)+a*x(n-1))/(1+a)
	y(n)=0
	return
	end
c
	subroutine mid(x1,x0,y1,y0,vx1,vx0,vy1,vy0,dt,Am)
	x1=x0+vx0*dt*0.5
	y1=y0+vy0*dt*0.5
	f=Am*sqrt(vx0**2+vy0**2)
	vx1=vx0-f*vx0*0.5
	vy1=vy0-4.9*dt-f*vy0*0.5
	return
	end
c
	subroutine display(x,y,nmax,theta,Am,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
	call text(0.2,5.8,0.2,'Cannonball in air: Euler',0.,24,0)
	else
	call text(0.2,5.8,0.2,'Cannonball in air: Runge-Kutta',0.,30,0)
	endif
	call text(0.2,5.5,0.2,'Height vs Distance',0.,20,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,-1)
	call axisx(0.,0.,'Height (m)',11,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(1.5,2.0,0.2,theta,0.,'''theta = '',f5.2')
	call number(1.5,2.3,0.2,Am,0.,'''A_m= '',1pe8.2')
	call number(1.5,2.6,0.2,dt,0.,'''dt = '',f5.2')
	call plot(0.,0.,999)
	return
	end

