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) character*12 out c c Use subroutines to do the work c call initialize(dt,vinit,om,th0,mp,bm,sm,out,method) call calculate(x,z,dt,vinit,om,th0,mp,bm,sm,n,method) call save(x,z,n,om,th0,mp,bm,sm,dt,out,method) stop end c subroutine initialize(dt,vinit,om,th0,mp,bm,sm,out,method) character ans,yes character*12 out c c Initialize variables c 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 output data file name ?' read(5,14) ans 14 format(a1) if(ans.eq.yes) then print *,'name ? -> ' read(5,15) out 15 format(a12) else out = 'pitch.data' 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 save(x,y,nmax,om,th0,mp,bm,sm,dt,out,method) character*12 out dimension x(1),y(1) open(1,file=out) rewind(1) write(1,10) (x(i),y(i),i=1,nmax) 10 format(1x,1p,e12.5,2x,e12.5) return end