Program PPlane * Program to compute the path of a paper plane being influenced by * gravity, lift and drag. implicit none * Parameters need for program (Physical constants) and the "constants" * in the program (Mass, XArea, Cd, Targ_dist, Terr) include 'PPlane.h' integer*4 max_iter ! Maximum number of iterations allowed parameter ( max_iter = 1000 ) * MAIN PROGRAM VARIABLES complex*16 Vel ! Complex X.Z components of velocities real*8 Int_dist ! Distance travel computed by numerial ! integration (m) real*8 Trav_time ! Travel time (sec) real*8 Dist1, Dist2 ! Distance traveled by plane (m) real*8 step ! Integration step size (secs). We initially ! set this value but check that it is OK and ! reset as needed. **** Let the user know what this program does: write(*,120) 120 format(/,'12.010 Program to solve Paper Plane Trajectory ', . 'problem',/, . 'Given initial height and velocity, program computes ', . 'the trajectory of the plane and distance traveled ',/, . 'to reach the ground') **** Find out from the user the necessary inputs. We use a variety * of input schemes here. The characteristics of the object can * be put into a file thus allowing easy runs with different * distances and launch angles. call read_input ***** Report to the user the setup parameters call report_setup **** Starting with the initial value of Launch_vel, iterate the * solution until we hit the targ_dist. In case there is a problem * keep track of number of iterations so that the program will * eventually end. Int_dist = 0.d0 step = 0.1d0 ! Initially set step size to 1-sec (tested in int_soln) outT = .false. **** Write out the solution write(*,220) Launch_vel, step 220 format('Initial Velocity needed ',F12.3,' m/s',/, . 'Step size ',F6.3,' sec' ) outT = .false. Trav_time = 0.d0 do while ( abs(dist1-dist2).gt.terr ) Vel = dcmplx(Launch_vel*cos(Ltheta_deg*deg_to_rad), . Launch_vel*sin(Ltheta_deg*deg_to_rad)) call runge_kutta(Vel, step, Dist1, Trav_time ) step = step/2 Vel = dcmplx(Launch_vel*cos(Ltheta_deg*deg_to_rad), . Launch_vel*sin(Ltheta_deg*deg_to_rad)) call runge_kutta(Vel, step, Dist2, Trav_time ) print *,'Step ',step, dist1,dist2, abs(dist1-dist2)*1000.0 end do outT = .true. Trav_time = 0.d0 Vel = dcmplx(Launch_vel*cos(Ltheta_deg*deg_to_rad), . Launch_vel*sin(Ltheta_deg*deg_to_rad)) call runge_kutta(Vel, step, Dist1, Trav_time ) **** Thats all end CTITLE REPORT_SETUP subroutine report_setup implicit none * Routine to report the parameters being used in the run * Parameters need for program (Physical constants) include 'PPlane.h' * **** Write out all the major parameters: write(*,120) Launch_hgt, Terr*1000, Ltheta_deg 120 format(/,'Launch Height: ',F5.3,' m, within ',F6.1,' mm, ', . 'Launch Angle ',F6.2,' deg') write(*,140) Mass, Area, Cl, Cd 140 format('MASS ',f7.4,' kg, Area ',F6.3, ' m^2, Cl ',F7.4, . ' Cd ',F7.4) write(*,160) g_0, rho_air 160 format('GRAVITY ',F8.5,' m/s**2 ', . 'RHO AIR ',F8.5,' kg/m**3') ***** Thats all return end CTITLE READ_INPUT subroutine read_input implicit none * Routine to read the input parameters of the program include 'PPlane.h' integer*4 ierr ! IOSTAT error reading user inputs * LOCAL VARIABLES real*8 terr_mm ! Integration accuracy mm real*8 Mass_gm ! Mass in grams **** Set default values. User can then get defaults by typing / Terr = 0.001d0 ! 1 mm error Ltheta_deg = -10.d0 ! Initial launch angle Launch_vel = 3.5d0 ! m/s Launch_hgt = 2.d0 ! m Mass = 0.003d0 ! kg Area = 0.017 ! m**2 Cl = 0.22d0 ! Lift coefficient m**2 Cd = 0.04d0 ! Drag coefficient m**2 **** Get input from user: write(*,120) 120 format(/,'Program Parameter Input. [Defaults are printed ', . 'use / to accept default') Terr_mm = terr*1000.d0 write(*,140) Launch_hgt, Terr_mm 140 format('Launch Height (m) and error (mm) [',F6.3,' km', . F5.2,' mm] ',$) read(*,*,iostat=ierr) Launch_hgt, Terr_mm call report_error(ierr,'Reading launch height and error') terr = terr_mm/1000.d0 write(*,160) Launch_vel, Ltheta_deg 160 format('Launch Velocity [',F6.3,' m/s] and Launch Angle (deg) [', . F6.2,' deg] ',$) read(*,*, iostat=ierr) Launch_vel, Ltheta_deg call report_error(ierr,'Reading launch angle') Mass_gm = Mass*1000.d0 write(*,180) Mass_gm, Area, Cl, Cd 180 format('Plane Mass (g), Area (m**2), Lift and Drag Coefficient ,', . '[',F8.1,' g, ',F6.3,' m^2, ', F8.5,' and ', . F8.5,' mm**2] ',$) read(*,*,iostat=ierr) Mass_gm, Area, Cl, Cd call report_error(ierr,'Reading Mass, Xarea and Cd') Mass = Mass_gm/1000.d0 return end CTITLE RUNGE_KUTTA subroutine runge_kutta(Vel, step, Dist, int_time ) implicit none * Parameters need for program (Physical constants) include 'PPlane.h' * PASSED VARIABLES complex*16 Vel ! X (horizontal) and Z (vertical) velocity of ! object (m/s) real*8 step ! Integration step size (secs). real*8 Dist ! Computed distances in time for Z-coordinate ! to return to zero. real*8 int_time ! Integration time. If this is zero initially ! integrate on position, otherwise integrate on ! on time. * LOCAL VARIABLES complex*16 Pos ! Complex representation of 2-components complex*16 P, V ! Intermediate position and velocities needed complex*16 k1, k2, k3, k4 ! Coefficients needed for Runge Kutta integration complex*16 accel ! Complex function that return x acceleration in real part ! and z-acceleration imaginary part real*8 time,T ! Time into integration real*8 tpast ! Time past the zero point for the last step. real*8 ustep ! Used step for integration. As the object approaches the ! ground, we take smaller steps in the last interval. This is ! is judged logical done ! Set true when the integration is complete based either ! on position or flight time integer*4 n ***** Initialize position Pos = dcmplx(0.0d0,Launch_hgt) time = 0.d0 done = .false. ustep = step if( step.eq.0 ) then print *,' Zero step ' step = 0.05d0 ustep = step end if n = 0 if( outT ) then write(*,160) time, Pos, Vel end if do while ( .not.done ) n = n + 1 if( n.gt.5000 ) done = .true. **** Compute the accelerations we need k1 = ustep*accel(time, Pos, Vel ) P = Pos + ustep*Vel/2.d0 +ustep*k1/8.d0 V = Vel + k1/2.d0 k2 = ustep*accel(time+ustep/2.d0, P, V) P = Pos + ustep*V/2.d0 + ustep*k1/8.d0 V = Vel + k2/2.d0 k3 = ustep*accel(time+ustep/2.d0, P, V) P = Pos + ustep*Vel + ustep*k3/2.d0 V = Vel + k3 k4 = ustep*accel(time+ustep, P, V) **** Update the time, position and velocity T = time + ustep P = Pos + ustep*(Vel + (k1+k2+k3)/6.d0) V = Vel + (k1+2*k2+2*k3+k4)/6.d0 **** See if we have passed below zero. If we have then * take smaller steps. if( Imag(P).lt.0 ) then tpast = Imag(Pos)/Imag(Vel) ustep = abs(tpast/2.d0) if( abs(Imag(P)).lt. terr/10 ) then time = T Pos = P Vel = V done = .true. if( outT ) then write(*,160) time, Pos, Vel end if endif else * Continue with integration time = T Pos = P Vel = V if( outT ) then write(*,160) time, Pos, Vel 160 format('T ',F10.3, 2F15.4, 2F15.4) end if end if **** See if we have completed. Based either on position * of time if( int_time.ne.0 ) then * Test on time should be = but to avoid rounding * error we test when is it within step/10. if( abs(time-int_time).lt.step/10 ) done =.true. end if end do **** Return the Distance traveled. Correct for the overshoot i.e., * the last Z-position is below the surface, so compute how many * seconds to get back to surface at last velocity and correct the * horizontal distance for this effect. tpast = Imag(Pos)/Imag(Vel) c Dist = realpart(Pos) - realpart(Vel)*tpast Dist = dble(Pos) - dble(Vel)*tpast int_time = time ***** Thats all return end CTITLE ACCEL complex*16 function accel( time, P, V ) implicit none **** Routine to compute the acceleration of the body at position * P and with Velocity V. Time is passed into the routine but * ir is not needed for the forces considered here. (In a more * general problem with mass wastage we may need time to compute * mass). * Parameters need for program (Physical constants) include 'PPlane.h' * PASSED VARIABLES real*8 time ! Integration time (not used in this version) complex*16 P, V ! X (horizontal) and Z (vertical) position and ! velocity of object (m, m/s) in complex form ! so that vector calculations can be done. * LOCAL VARIABLES complex*16 gacc, dacc, lacc ! Gravity, Drag and Lift accelerations complex*16 ldir ! Direction perpendicular to motion direction real*8 Vmag ! Magnitude of velocity real*8 Rho_h ! Density of air at our height ***** Compute the accelations: * GRAVITY: Depends of Height gacc = dcmplx(0.d0, -g_0) * DRAG: Acts opposite to the motion vector (since wind is fixed * relative to motion direction) Vmag = abs(V) Rho_h = rho_air dacc = -(Cd*Rho_h*Vmag*V/2)/Mass*Area ! Note: This is Vmag^2*V_unit_vector * LIFT: Acts perpendicular to the motion vector ldir = dcmplx(-Imag(V),dble(V))/Vmag lacc = (Cl*Rho_h*Vmag**2*ldir)/2/Mass*Area **** Add the acceleration together accel = gacc + dacc + lacc c accel = dcmplx(0.d0, -g_0 ) ***** Thats all return end CTITLE REPORT_ERROR subroutine report_error(ierr,mess) * Routine to report IOSTAT errors. Initially developed * for 12.010 HW2 Problem 3. * PASSED VARIABLES integer*4 ierr character*(*) mess * LOCAL VARIABLES * None ***** See if the IOSTAT error was non-zero if( ierr.ne.0 ) then write(*,120) ierr, mess 120 format('IOSTAT Error ',i4,' occurred ',a) stop 'HW2_3: IO Error in program' end if return end