c
c Simulation to obtain the maximum range of a large cannon
c        - Euler or Runge-Kutta (2nd or 4th order usual ones)
c Fortran version written by H. Nakanishi, output to files
c Changed the density correction to be appropriate for troposphere.
c
	program maxrange
c
c       Declare the arrays we will need
c
	implicit real*8 (a-h,o-z)
	character ans,yes
	yes='y'
c
c       Use subroutines to do the work
c
20	call init(dt,vinit,th0,err,Am,gm,g,ng,y0,method)   
	dth=10.d0
	th=th0-dth
	range0=0.d0
	do 30 i=1,50
	th=th+dth
	call calc(dt,vinit,th,Am,gm,g,ng,y0,n,method,range)
	if(range.le.range0) then
		dth=-0.5d0*dth
		if(4.*abs(dth).lt.err) go to 40
	endif
	range0=range
30	continue
	print *,'failed to obtain maximum range...'
	go to 50
40	print *,'Maximum range = ',range,' m',' at ',th,' degrees'
	print *,i,' iterations, dth=',dth
50      print *,'another try?'
        read(5,10) ans
10      format(a1)
        if(ans.eq.yes) go to 20
	stop
	end
c
	subroutine init(dt,vinit,th0,err,Am,gm,g,ng,y0,method)
c
c       Initialize variables
c
	implicit real*8 (a-h,o-z)
	print *,'Euler (1) or Runge-Kutta 2nd order (2), 4th (3)? -> '
	read(5,*) method
	if(method.ne.1.and.method.ne.2.and.method.ne.3) then
		print *,'must select 1, 2, or 3 ..'
		stop
	endif
	print *,'initial v, time step, initial angle, error (deg) -> '
	read(5,*) vinit,dt,th0,err
	print *,'g and drag/m at sea level, -dT/dh, gamma -> '
	read(5,*) g,Am,y0,gm
	print *,'vary g with altitude [1=yes,0=no] ?'
	read(5,*) ng
	return
	end
c
	subroutine calc(dt,vinit,theta,Am,gm,g,ng,y0,n,method,range)
c
c       Now use the Euler method or the Runge-Kutta (2nd order)
c
	implicit real*8 (a-h,o-z)
	xx1=0
	yy1=0
	vx=vinit*dcos(3.141592*theta/180.)
	vy=vinit*dsin(3.141592*theta/180.)
	if(method.eq.1) then
10		continue
		xx0=xx1
		yy0=yy1
		call der(xx1,yy1,vx,vy,0.d0,Am,gm,g,ng,y0,dx,dy,dvx,dvy)
		xx1=xx1+dt*dx
		yy1=yy1+dt*dy
		vx=vx+dt*dvx
		vy=vy+dt*dvy
		if(yy1.le.0) then
			n=i
			go to 15
		endif
		go to 10
	elseif(method.eq.2) then
30		continue
		xx0=xx1
		yy0=yy1
		call der(xx1,yy1,vx,vy,0.d0,Am,gm,g,ng,y0,dx,dy,dvx,dvy)
		x1=xx1+0.5d0*dt*dx
		y1=yy1+0.5d0*dt*dy
		vx1=vx+0.5d0*dt*dvx
		vy1=vy+0.5d0*dt*dvy
		call der(x1,y1,vx1,vy1,0.d0,Am,gm,g,ng,y0,dx2,dy2,dvx2,dvy2)
		xx1=xx1+dt*dx2
		yy1=yy1+dt*dy2
		vx=vx+dt*dvx2
		vy=vy+dt*dvy2
		if(yy1.le.0) then
			n=i
			go to 15
		endif
		go to 30
	else
40		continue
		xx0=xx1
		yy0=yy1
		call der(xx1,yy1,vx,vy,0.d0,Am,gm,g,ng,y0,dx,dy,dvx,dvy)
		x1=xx1+0.5d0*dt*dx
		y1=yy1+0.5d0*dt*dy
		vx1=vx+0.5d0*dt*dvx
		vy1=vy+0.5d0*dt*dvy
		call der(x1,y1,vx1,vy1,0.d0,Am,gm,g,ng,y0,dx2,dy2,dvx2,dvy2)
		x2=xx1+0.5d0*dt*dx2
		y2=yy1+0.5d0*dt*dy2
		vx2=vx+0.5d0*dt*dvx2
		vy2=vy+0.5d0*dt*dvy2
		call der(x2,y2,vx2,vy2,0.d0,Am,gm,g,ng,y0,dx3,dy3,dvx3,dvy3)
		x3=xx1+dt*dx3
		y3=yy1+dt*dy3
		vx3=vx+dt*dvx3
		vy3=vy+dt*dvy3
		call der(x3,y3,vx3,vy3,0.d0,Am,gm,g,ng,y0,dx4,dy4,dvx4,dvy4)
		xx1=xx1+0.16666667*dt*(dx+2*dx2+2*dx3+dx4)
		yy1=yy1+0.16666667*dt*(dy+2*dy2+2*dy3+dy4)
		vx=vx+0.16666667*dt*(dvx+2*dvx2+2*dvx3+dvx4)
		vy=vy+0.16666667*dt*(dvy+2*dvy2+2*dvy3+dvy4)
		if(yy1.le.0) then
			n=i
			go to 15
		endif
		go to 40
	endif
15	a=-yy1/yy0
	range=(xx1+a*xx0)/(1.d0+a)
	return
	end
c
	subroutine der(x0,y0,vx0,vy0,t0,Am,gm,g,ng,yy,dx,dy,dvx,dvy)
	implicit real*8 (a-h,o-z)
	dx=vx0
	dy=vy0
	if(yy.eq.0) then
		f=Am*dsqrt(vx0**2+vy0**2)
	else
		f=Am*dsqrt(vx0**2+vy0**2)*(1-yy*y0/300)**(1/(gm-1))
	endif
	dvx=-f*vx0
	if(ng.ne.1) then
		dvy=-f*vy0-g
	else
		r0=6.38d6
		dvy=-f*vy0-g*(r0/(y0+r0))**2
	endif
	return
	end
