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) c c Use subroutines to do the work c call initialize(dt,vinit,om,th0,mp,bm,sm,lsym,nsym,method) call calculate(x,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) 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 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,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 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 call text(0.2,5.8,0.2,'Curve/Knuckle ball: Euler',0.,25,0) elseif(method.eq.2) then call text(0.2,5.8,0.2,'Curve/Knuckle ball: Runge-Kutta2',0.,32,0) 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