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 Uses Golden Section algorithm. c program maxgold 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) call calc(dt,vinit,th0,Am,gm,g,ng,y0,n,method,range0) gold=0.5d0*(dsqrt(5.d0)-1.d0) dth=10.d0 th=th0+dth call calc(dt,vinit,th,Am,gm,g,ng,y0,n,method,range1) if(range1.le.range0) then temp=range0 range0=range1 range1=temp dth=-dth th=th0 endif isuc=0 do 30 i=1,50 th=th+dth call calc(dt,vinit,th,Am,gm,g,ng,y0,n,method,range2) if(range2.le.range1) then dth2=-gold*dth th2=th+dth2 call calc(dt,vinit,th2,Am,gm,g,ng,y0,n,method,range3) if(range3.gt.range1) then range0=range1 range1=range2 range=range3 else th2=th2-dth call calc(dt,vinit,th2,Am,gm,g,ng,y0,n,method,range4) if(range4.gt.range1) then th=th-dth range=range4 else temp=th2 th2=th-dth th=temp dth=-dth range0=range3 range=range1 range1=range4 endif endif call golden(dt,vinit,Am,gm,g,ng,y0,n,method,iter c ,range0,range1,range,th,th2,dth,gold,err,isuc) go to 35 else range0=range1 range1=range2 endif 30 continue 35 if (isuc.gt.0) then print *,'Maximum range = ',range,' m',' at ',th2,' degrees' print *,iter,' iterations, dth=',dth else print *,'failed to obtain maximum range...' print *,iter,' iterations, dth=',dth endif 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 c subroutine golden(dt,vinit,Am,gm,g,ng,y0,n,method,iter c ,range0,range1,range,th,th2,dth,gold,err,isuc) implicit real*8 (a-h,o-z) iter=0 10 iter=iter+1 if(dabs(dth).lt.err) then isuc=1 return elseif(iter.ge.100) then isuc=0 return else th3=th-dth*gold**2 call calc(dt,vinit,th3,Am,gm,g,ng,y0,n,method,range3) if(range3.gt.range) then th2=th3 dth=dth*gold range0=range range=range3 else th=th-dth dth=-dth*gold range1=range0 range0=range3 endif go to 10 endif return end