      PROGRAM TO2
***********************************************************************
*     This program calculates takeoff for commercial jet aircraft.    *
*     Both All Engines Operating (AEO) (normal) and One Engine Out    *
*     (OEI) takeoff are calculated.  OEI calculations are for         *
*     Balanced Field Length (BFL).                                    *
*                                                                     *
*     Methodology:  A.R. Krenkel and A. Salzman                       *
*                   "Journal of Aircraft", Vol. 5, No. 5,             *
*                   Sept-Oct 1968, pp. 429-436                        *
*                                                                     *
*     Algorithm:    Sean R. Lynn                                      *
*                   Original Scheme March 19, 1994                    *
*                   Revision of BFL methodology May 11, 1994          *
*                                                                     *
*     The original methodology solved balance of forces equations of  *
*     motions for the aircraft parametrically.  This program uses     *
*     time-step integration of a 4th order Runge-Kutta scheme to solve*
*     the equations.  The Runge-Kutta scheme is semi-adaptive,        *
*     reducing the user-input step size until the error tolerance is  *
*     met.                                                            *
*                                                                     *
*     The program flows as follows:                                   *
*     1) Initialize variables and open i/o devices.                   *
*     2) Read in user aircraft data from TAKEOFF2.IN using subroutine *
*        GETDAT.  Put input data into array DATARY(24)                *
*     3) Set i/o devices, calculate quadratic thrust coefficients     *
*        using QINTRP, calculate the rotation velocity and set the    *
*        per step error tolerance.                                    *
*     4) Echo input data to user-specified i/o device                 *
*     5) Solve for normal takeoff.                                    *
*        a) Incremental output.                                       *
*           i)   V=0 to V=Vr (up to rotation)                         *
*           ii)  V=Vr to V=Vlo (rotation for user-specified time)     *
*           iii) Vlo to Vobs (climb to user-specified obstacle height)*
*        b) Summary output of normal takeoff.                         *
*     6) Solve for OEI BFL.                                           *
*        a) Two initial guesses are made of the critical engine       *
*           failure velocity.                                         *
*        b) Calculate both OEI takeoff and aborted takeoff for both   *
*           guess velocities.  For BFL these two distances should be =*
*        c) Use bisection method to converge on a Vcrit such that     *
*           the BFL definition is solved                              *
*        d) Output summary data of BFL calculation                    *
*                                                                     *
*     Important subroutines:                                          *
*     GETDAT - gets user input data from TAKEOFF2.IN                  *
*     QINTRP - finds the quadratic coefficients for thrust using      *
*              polynomial interpolation                               *
*     LINTRP - performs linear interpolation between 2 points.  Used  *
*              to interpolate between incremental velocities.         *
*     GROUND - the equations of motion for the ground roll segment    *
*     AIR    - the equations of motion for the climb segment          *
*     ADPTRK - reduces user-input time step size to make sure error   *
*              tolerance is met. (This is the semi-adaptive portion)  *
*     RK34   - third and fourth order Runge-Kutta scheme used in      *
*              conjunction with ADPTRK to solve aircraft EOMs			*  
*																	*
*	Mods:															*
*			1. Added a routine to prompt the user for input and     *
*			   output filenames, Leifur Thor Leifsson 03-02-04.     *
*			2. A header with information on the program is written  *
*			   to screen/file/printer,Leifur Thor Leifsson 03-02-04.*
***********************************************************************

      REAL DATARY(24),T0,T1,T2,VLO,RPARM(11),ERR,Y4(4),Y3(4),T,TOUT,
     .     TOLD,Y1OLD,Y2OLD,TEMP1,TEMP2,LINTRP,DYDX(4),VR,VOBS,XR,XLO,
     .     XOBS,TR,TLO,TOBS,XDIF1,XDIF2,TGES,TGES1,TGES2,Y3OLD,XDIF,
     .     XCRIT,VCRIT,TCRIT,TDESC,VDESC,XDESC,BFL,TFL
      INTEGER INFO(1),IOUNIT,IPARM(1),ONETWO,I
      CHARACTER*40 TITLE

      EXTERNAL GROUND, AIR

      OPEN(6,FILE='PRN')
C      OPEN(7,FILE='TAKEOFF2.OUT')
      OPEN(8,FILE='CON')

C     Get aircraft data from the input file
      CALL GETDAT(DATARY,INFO,TITLE)

C     Set the output device
      IOUNIT = INFO(1)

C     Find quadratic thrust coefficients T0,T1,T2
      CALL QINTRP(DATARY(20),DATARY(17),T0,T1,T2)

C     Calculate the rotation velocity VR
      VR = DATARY(13) * SQRT((2.*DATARY(3))/(DATARY(2)*DATARY(4)*
     .      DATARY(5)))

C     Set step-wise tolerance for Runge-Kutta scheme
      ERR = 1.E-6

C     Write out case title, headings and echo input data
      WRITE (IOUNIT,'(20X,A)') TITLE
      WRITE (IOUNIT,'(//)')
      WRITE (IOUNIT,'(A)') '<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< CASE INPUT D
     .ATA >>>>>>>>>>>>>>>>>>>>>>>>>>>>>'
      WRITE (IOUNIT,*)
      WRITE (IOUNIT,100)'Gravitational accel. (g)=',DATARY(1),'(ft/s^2)'
      WRITE(IOUNIT,100)'Air density at takeoff  =',DATARY(2),'(sl/ft^3)'
      WRITE (IOUNIT,100)'Aircraft weight      (W)=',DATARY(3),'(lbs)'
      WRITE (IOUNIT,100)'Aircraft wing area   (S)=',DATARY(4),'(ft/s^2)'
      WRITE (IOUNIT,100)'Max lift coeff   (CLmax)=',DATARY(5),'(nondim)'
      WRITE (IOUNIT,100)'Ground lift coeff(CLgrd)=',DATARY(6),'(nondim)'
      WRITE (IOUNIT,100)'Air lift coeff   (CLair)=',DATARY(7),'(nondim)'
      WRITE (IOUNIT,100)'Ground drag coeff(CDgrd)=',DATARY(8),'(nondim)'
      WRITE (IOUNIT,100)'Air drag coeff   (CDair)=',DATARY(9),'(nondim)'
      WRITE(IOUNIT,100)'Roll frict coeff (MUgrd)=',DATARY(10),'(nondim)'
      WRITE(IOUNIT,100)'Brake frict coeff(MUbrk)=',DATARY(11),'(nondim)'
      WRITE(IOUNIT,100)'Angle of thrust (LAMBDA)=',DATARY(12),'(rad)'
      WRITE(IOUNIT,100)'Stall margin         (k)=',DATARY(13),'(nondim)'
      WRITE(IOUNIT,100)'Descision time    (TIME)=',DATARY(14),'(sec)'
      WRITE(IOUNIT,100)'Obstacle height  (OBSHT)=',DATARY(15),'(ft)'
      WRITE(IOUNIT,100)'EI pwr remaining (PLOSS)=',DATARY(16),'(fract)'
      WRITE(IOUNIT,*)
      WRITE(IOUNIT,110)'THRUST = ',T0,' + ',T1,'*V + ',T2,'*V^2'
      WRITE (IOUNIT,'(//)')

C     Solve for normal takeoff
c     Set variables passed to subroutines: rho,W,S,CL,CD,mu,T0,T1,T2,
c                                          LAMDA,g - used for grd roll
      RPARM(1) = DATARY(2)
      RPARM(2) = DATARY(3)
      RPARM(3) = DATARY(4)
      RPARM(4) = DATARY(6)
      RPARM(5) = DATARY(8)
      RPARM(6) = DATARY(10)
      RPARM(7) = T0
      RPARM(8) = T1
      RPARM(9) = T2
      RPARM(10) = DATARY(12)
      RPARM(11) = DATARY(1)

c     Ground roll
c     Set number of equations for call to ADPTRK and init dist,spd,time
      NEQ = 2
      Y4(1) = 0.
      Y4(2) = 0.
      Y4(3) = 0.
      Y4(4) = 0.
      T = 0.

c     Find solution between V=0 and V=Vr
      WRITE (IOUNIT,'(A)') '<<<<<<<<<<<<<<<<<<<<<<<<<<<<< NORMAL TAKEOFF
     . DATA >>>>>>>>>>>>>>>>>>>>>>>>>'
      WRITE (IOUNIT,'(A)')'Time (s)      X-Dist (ft)      X-Vel (ft/s)
     .   Y-Dist (ft)      Y-Vel (ft/s)'
      WRITE (IOUNIT,120) T,Y4(1),Y4(2),Y4(3),Y4(4)
      DO 10 I = 1,1000
        TOUT = T + DATARY(23)
        TOLD = T
        Y1OLD = Y4(1)
        Y2OLD = Y4(2)
        CALL ADPTRK(GROUND,T,TOUT,Y4,Y3,DYDX,NEQ,IPARM,RPARM,ERR)
        IF (Y4(2) .LT. VR) THEN
          WRITE (IOUNIT,120) TOUT,Y4(1),Y4(2),Y4(3),Y4(4)
        ELSE
          TEMP1 = LINTRP(VR,Y2OLD,Y4(2),TOLD,TOUT)
          TEMP2 = LINTRP(VR,Y2OLD,Y4(2),Y1OLD,Y4(1))
          WRITE (IOUNIT,120) TEMP1,TEMP2,VR,Y4(3),Y4(4)
          TR = TEMP1
          XR = TEMP2
          T = TOUT
          GOTO 20
        END IF
        T = TOUT
10    CONTINUE

c     Find solution for rotation phase (ground run for T=TIME seconds)
20    Y4(1) = TEMP2
      Y4(2) = VR
      TEMP2 = TEMP1 + DATARY(24)
      CALL ADPTRK(GROUND,TEMP1,TEMP2,Y4,Y3,DYDX,NEQ,IPARM,RPARM,ERR)
      WRITE (IOUNIT,120) TEMP1,Y4(1),Y4(2),Y4(3),Y4(4)
      TLO = TEMP1
      XLO = Y4(1)
      VLO = Y4(2)

c     Climb phase
c     Change vars. passed to subroutines: CL,CD - used for climb
      RPARM(4) = DATARY(7)
      RPARM(5) = DATARY(9)

c     Set number of equations for call to ADPTRK and init dist,spd,time
      NEQ = 4

c     Solve from end of rotation to next even increment
      T = T + DATARY(24)
      CALL ADPTRK(AIR,TEMP2,T,Y4,Y3,DYDX,NEQ,IPARM,RPARM,ERR)
      WRITE (IOUNIT,120) T,Y4(1),Y4(2),Y4(3),Y4(4)

c     Solve from even increment to OBSHT
      DO 30 I = 1,1000
        TOUT = T + DATARY(23)
        TOLD = T
        Y1OLD = Y4(1)
        Y2OLD = Y4(2)
        Y3OLD = Y3(3)
        Y4OLD = Y4(4)
        CALL ADPTRK(AIR,T,TOUT,Y4,Y3,DYDX,NEQ,IPARM,RPARM,ERR)
        IF (Y4(3) .LT. DATARY(15)) THEN
          WRITE (IOUNIT,120) TOUT,Y4(1),Y4(2),Y4(3),Y4(4)
        ELSE
          T     = LINTRP(DATARY(15),Y3OLD,Y4(3),TOLD,TOUT)
          Y4(1) = LINTRP(DATARY(15),Y3OLD,Y4(3),Y1OLD,Y4(1))
          Y4(2) = LINTRP(DATARY(15),Y3OLD,Y4(3),Y2OLD,Y4(2))
          Y4(4) = LINTRP(DATARY(15),Y3OLD,Y4(3),Y4OLD,Y4(4))
          WRITE (IOUNIT,120) T,Y4(1),Y4(2),DATARY(15),Y4(4)
          TOBS = T
          XOBS = Y4(1)
          VOBS = SQRT(Y4(2)*Y4(2)+Y4(4)*Y4(4))
          GOTO 40
        END IF
        T = TOUT
30    CONTINUE

c     Write normal takeoff data to output device
40    WRITE (IOUNIT,*)
      WRITE (IOUNIT,'(A)') '<<<<<<<<<<<<<<<<<<<<<<<<<< NORMAL TAKEOFF SU
     .MMARY >>>>>>>>>>>>>>>>>>>>>>'
      WRITE (IOUNIT,*)
      WRITE (IOUNIT,100)'Rotation Velocity   (Vr)=',VR,'(ft/s)'
      WRITE (IOUNIT,100)'Liftoff Velocity   (Vlo)=',VLO,'(ft/s)'
      WRITE (IOUNIT,100)'Velocity over obs.(Vobs)=',VOBS,'(ft/s)'
      WRITE (IOUNIT,100)'Rotation Distance   (Xr)=',XR,'(ft)'
      WRITE (IOUNIT,100)'Liftoff Distance   (Xlo)=',XLO,'(ft)'
      WRITE (IOUNIT,100)'Distance to obst. (Xobs)=',XOBS,'(ft)'
      WRITE (IOUNIT,100)'Rotation Time       (Tr)=',TR,'(s)'
      WRITE (IOUNIT,100)'Liftoff Time       (Tlo)=',TLO,'(s)'
      WRITE (IOUNIT,100)'Time to obst.     (Tobs)=',TOBS,'(s)'
      WRITE (IOUNIT,*)
      WRITE (IOUNIT,100)'TOTAL TAKEOFF DIST (Xto)=',XOBS,'(ft)'
      WRITE (IOUNIT,100)'TOTAL TAKEOFF TIME (Tto)=',TOBS,'(s)'
      WRITE (IOUNIT,*)



C     FIND THE BALANCED FIELD LENGTH
c     Find the critical engine failure velocity, Vcrit
      WRITE (8,'(A)') 'PLEASE WAIT, CONVERGING ON BFL!:  USU @ 2MIN'

c     Set the inital Vcrit guess velocities
      TGES1 = .25*TR
      TGES2 = .95*TR

c     Solve for engine out takeoff at both guess times and use the
c     bisection method to iterate and converge on Xbrak = XOEI

c     Iterate to converge on Vcrit
      DO 93 K = 1, 1000
c     Iterate through two guesses for bisection routine
      DO 92 ONETWO = 1,2

      IF (ONETWO .EQ. 1) THEN
        TGES = TGES1
      ELSE
        TGES = TGES2
      END IF

      RPARM(1) = DATARY(2)
      RPARM(2) = DATARY(3)
      RPARM(3) = DATARY(4)
      RPARM(4) = DATARY(6)
      RPARM(5) = DATARY(8)
      RPARM(6) = DATARY(10)
      RPARM(7) = T0
      RPARM(8) = T1
      RPARM(9) = T2
      RPARM(10) = DATARY(12)
      RPARM(11) = DATARY(1)
      NEQ = 2
      Y4(1) = 0.
      Y4(2) = 0.
      Y4(3) = 0.
      Y4(4) = 0.
      T = 0.

C     Solution for OEI ground roll
c     Solution between V=0 and V=Vcrit (Engine Failure)
      CALL ADPTRK(GROUND,T,TGES,Y4,Y3,DYDX,NEQ,IPARM,RPARM,ERR)

      XCRIT = Y4(1)
      VCRIT = Y4(2)
      TCRIT = TGES
      RPARM(7) = T0*DATARY(16)
      RPARM(8) = T1*DATARY(16)
      RPARM(9) = T2*DATARY(16)
      TOLD = TGES

c     Solution from V=Vcrit to V=Vr
      DO 50 I = 1, 1000
      T = TOLD + 0.5
      Y1OLD = Y4(1)
      Y2OLD = Y4(2)
      CALL ADPTRK(GROUND,TOLD,T,Y4,Y3,DYDX,NEQ,IPARM,RPARM,ERR)
      IF (Y4(2) .GE. VR) THEN
        TEMP1 = LINTRP(VR,Y2OLD,Y4(2),Y1OLD,Y4(1))
        TEMP2 = VR
        TEMP3 = LINTRP(VR,Y2OLD,Y4(2),TOLD,T)
        GOTO 60
      END IF
      TOLD = T
50    CONTINUE

c     Solution from V=Vr to V=Vlo (Engine Out Rotation)
60    Y4(1) = TEMP1
      Y4(2) = TEMP2
      T = TEMP3
      TOUT = T + DATARY(24)
      CALL ADPTRK(GROUND,T,TOUT,Y4,Y3,DYDX,NEQ,IPARM,RPARM,ERR)

c     Solution for OEI Climb phase
      NEQ = 4
      RPARM(4) = DATARY(7)
      RPARM(5) = DATARY(9)

      DO 70 I = 1, 1000
      TOLD = T
      T = T + .5
      Y1OLD = Y4(1)
      Y3OLD = Y4(3)
      CALL ADPTRK(AIR,TOLD,T,Y4,Y3,DYDX,NEQ,IPARM,RPARM,ERR)
      IF (Y4(3) .GE. DATARY(15)) THEN
        XDIF = LINTRP(DATARY(15),Y3OLD,Y4(3),Y1OLD,Y4(1))
        GOTO 80
      END IF
70    CONTINUE

C     Solution for OEI braking
c     Ground run V=Vcrit to V=V1 (Descision time)
80    NEQ = 2
      RPARM(4) = DATARY(6)
      RPARM(5) = DATARY(8)
      Y4(1) = XCRIT
      Y4(2) = VCRIT
      Y4(3) = 0.
      Y4(4) = 0.
      T = TGES
      TDESC = T + DATARY(14)
      CALL ADPTRK(GROUND,T,T+DATARY(14),Y4,Y3,DYDX,NEQ,IPARM,RPARM,ERR)
      XDESC = Y4(1)
      VDESC = Y4(2)

c     Braking Solution V=V1 to V=0.
c     Set thrust to 0. and friction coefficient to braking value
      RPARM(7) = 0.
      RPARM(8) = 0.
      RPARM(9) = 0.
      RPARM(6) = DATARY(11)
      T = TGES + DATARY(14)

      DO 90 I = 1, 1000
      TOLD = T
      T = T + .5
      Y1OLD = Y4(1)
      CALL ADPTRK(GROUND,TOLD,T,Y4,Y3,DYDX,NEQ,IPARM,RPARM,ERR)
      IF (Y4(2) .LE. 0.) THEN
        Y4(1) = LINTRP(0.,Y2OLD,Y4(2),Y1OLD,Y4(1))
        TFL = LINTRP(0.,Y2OLD,Y4(2),TOLD,T)
        GOTO 91
      END IF
90    CONTINUE

	

91    IF (ONETWO .EQ. 1) THEN
        XDIF1 = XDIF - Y4(1)
      ELSE
        XDIF2 = XDIF - Y4(1)
      END IF

92    CONTINUE

      IF (ABS(XDIF2) .LE. 1E-2) THEN		
        BFL = Y4(1)						
        GOTO 94
      ELSE
        TEMP1 = .8*LINTRP(0.,XDIF1,XDIF2,TGES1,TGES2)
        TEMP2 = LINTRP(0.,XDIF1,XDIF2,TGES1,TGES2)
      END IF
        TGES1 = TEMP1
        TGES2 = TEMP2
	
	IF (K == 1000) THEN
	  WRITE(IOUNIT,*) 'BFL CALCULATION DID NOT CONVERGE!'
	END IF
93    CONTINUE

c     Write OEI takeoff data to output device
94    WRITE (IOUNIT,*)
      WRITE (IOUNIT,'(A)') '<<<<<<<<<<<<<<<<<<<<<<<<<< OEI TAKEOFF SUMMA
     .RY >>>>>>>>>>>>>>>>>>>>>>'
      WRITE (IOUNIT,*)
      WRITE (IOUNIT,100)'Critical Velocity   (Vcrit) =',VCRIT,'(ft/s)'
      WRITE (IOUNIT,100)'Decision Velocity      (V1) =',VDESC,'(ft/s)'
      WRITE (IOUNIT,100)'Velocity over obs.   (Vobs) =',VOBS,'(ft/s)'
      WRITE (IOUNIT,100)'Critical Distance   (Xcrit) =',XCRIT,'(ft)'
      WRITE (IOUNIT,100)'Decision Distance      (X1) =',XDESC,'(ft)'
      WRITE (IOUNIT,100)'Balanced Field Length (BFL) =',BFL,'(ft)'
      WRITE (IOUNIT,100)'Critical Time       (Tcrit) =',TCRIT,'(s)'
      WRITE (IOUNIT,100)'Decision Time          (T1) =',TDESC,'(s)'
      WRITE (IOUNIT,100)'OEI Takeoff Time     (TBFL) =',TFL,'(s)'
      WRITE (IOUNIT,*)

	PAUSE 'Press ENTER to close program'

100   FORMAT(A30,F15.5,T44,A8)
110   FORMAT(A,3(F10.3,A))
120   FORMAT(F7.3,6X,4(F10.3,7X))
      END


      SUBROUTINE GETDAT(DATARY,IOINFO,TITLE)
C     This subroutine is used to load data from TAKEOFF2.IN
      REAL DATARY(24)
      INTEGER IOINFO(*),I
      CHARACTER*40 TITLE
	CHARACTER*20 INFILE
	CHARACTER*20 OUTFILE

	WRITE(8,3000)
	WRITE(8,3001)
	WRITE(8,3002)
	WRITE(8,3003)

	WRITE(8,1000) 
      READ(5,'(a)') INFILE
      OPEN(UNIT = 10, FILE = INFILE, STATUS = 'OLD')

C      OPEN (10, FILE='TAKEOFF2.IN')

c     Set acceleration due to gravity
      DATARY(1) = 32.174

c     Read in run title information
      READ (10,'(A)') TITLE

c     Read in the rest of the data (see TAKEOFF2.IN for order)
      DO 10 I = 2,16
        READ(10,*) DATARY(I)
10    CONTINUE
      READ (10,*) DATARY(17),DATARY(18),DATARY(19)
      READ (10,*) DATARY(20),DATARY(21),DATARY(22)
      READ (10,*) DATARY(23)
      READ (10,*) DATARY(24)
      READ (10,*) IOINFO(1)

	IF (IOINFO(1) == 7) THEN
		WRITE(8,1010) 
		READ(5,'(a)') OUTFILE
		OPEN(UNIT = 7, FILE = OUTFILE, STATUS = 'REPLACE')
		WRITE(7,3000)
		WRITE(7,3001)
		WRITE(7,3002)
		WRITE(7,3003)
		WRITE(7,3004)
	END IF

1000	FORMAT(/2x,'Enter input filename:'/)
1010	FORMAT(/2x,'Enter output filename:'/)
3000  FORMAT('  **************************************************'/
     1        '    VT Aerospace Aircraft Design Software Series'/
     2		'  **************************************************'/)
3001  FORMAT('                   Program Takeoff'/
     1		/'  This program calculates takeoff for commercial jet'/
     2		'  aircraft. Both All Engines Operating (AEO) (normal)'/
     3		'  and One Engine Out (OEI) takeoff are calculated. OEI'/
     4		'  calculations are for Balanced Field Length (BFL).'/ 
     5		/'  Methodology:  A.R. Krenkel and A. Salzman'/
     6		 '                Journal of Aircraft, Vol. 5, No. 5,'/            
     7         '                Sept-Oct 1968, pp. 429-436.'/)
3002  FORMAT('  Code written by Sean Lynn, March, 1994.'/
     1			'  Modified by Leifur Thor Leifsson, March, 2004.'/)
3003  FORMAT('  The Department of Aerospace and Ocean Engineering,'/
     1		'  Virginia Tech, Blacksburg, VA 24061.'/
     2		'  http://www.aoe.vt.edu, email: whmason@vt.edu')
3004  FORMAT(/'  **************************************************'/)

      CLOSE(10)
      RETURN
      END


      SUBROUTINE QINTRP(V,T,T0,T1,T2)
C     This subroutine uses polynomial interpolation to find the
C     quadratic contants for the engine thrust/velocity curve, given
C     the three dimensional arrays for thrust (T) and velocity (V).
      REAL V(3),T(3),T0,T1,T2,A,B,C

      A = T(1)
      B = (T(2) - T(1))/(V(2)-V(1))
      C = ((T(3)-T(2))/(V(3)-V(2))-(T(2)-T(1))/(V(2)-V(1)))/(V(3)-V(1))

      T0 = T(1) - B*V(1) + C*V(1)*V(2)
      T1 = B - C*(V(1)+V(2))
      T2 = C

      RETURN
      END


      FUNCTION LINTRP(XIN,X1,X2,Y1,Y2)
C     This function finds the linear interpolation between the points
C     [X0,Y0] and [X1,Y1].
      REAL X1,X2,Y1,Y2,LINTRP,XIN

      LINTRP = Y1 + ((Y2-Y1)/(X2-X1))*(XIN-X1)

      RETURN
      END


      SUBROUTINE GROUND(X,Y,IPARM,RPARM,DYDX)
C     These are the equations for ground roll
      REAL X,Y(*),RPARM(*),DYDX(*),CONST,K0,K1,K2,DUMMY
      INTEGER IPARM(*)

      IPARM(1) = 0
      DUMMY = X

C     Equation is of form du/dt = g/W(K0 + K1*u + K2*u*u), where
C     g/W is 1/mass, and K's are constants and u is the velocity
      CONST = COS(RPARM(10)) + RPARM(6)*SIN(RPARM(10))
      K0 = RPARM(7)*CONST - RPARM(6)*RPARM(2)
      K1 = RPARM(8)*CONST
      K2 = RPARM(9)*CONST + .5*RPARM(1)*RPARM(3)*(RPARM(6)*RPARM(4) -
     .     RPARM(5))

      DYDX(1) = Y(2)
      DYDX(2) = RPARM(11)/RPARM(2)*(K0+K1*Y(2)+K2*Y(2)*Y(2))

      RETURN
      END


      SUBROUTINE AIR(X,Y,IPARM,RPARM,DYDX)
C     These are the equations for CLIMB segment
      REAL X,Y(*),RPARM(*),DYDX(*),K11,K12,K13,K21,K22,K23,DUMMY
      INTEGER IPARM(*)

      IPARM(1) = 0
      DUMMY = X

C     Equations are of the form:  du/dt = g/W*[k11+k12*V+k13*V*V]
C                                 dv/dt = g/W*[k21+k22*V+k23*V*V]
C     where, V = sqrt (u*u + v*v) and tan gamma = v/u
      GAMMA = ATAN(Y(4)/Y(2))
      VSQR = Y(4)*Y(4) + Y(2)*Y(2)

      K11 = RPARM(7)*COS(RPARM(10) + GAMMA)
      K12 = RPARM(8)*COS(RPARM(10) + GAMMA)
      K13 = RPARM(9)*COS(RPARM(10) + GAMMA) - .5*RPARM(1)*RPARM(3)*
     .      (RPARM(5)*COS(GAMMA)+RPARM(4)*SIN(GAMMA))
      K21 = RPARM(7)*SIN(RPARM(10) + GAMMA) - RPARM(2)
      K22 = RPARM(8)*SIN(RPARM(10) + GAMMA)
      K23 = RPARM(9)*SIN(RPARM(10) + GAMMA) + .5*RPARM(1)*RPARM(3)*
     .      (RPARM(4)*COS(GAMMA)-RPARM(5)*SIN(GAMMA))

      DYDX(1) = Y(2)
      DYDX(2) = (RPARM(11)/RPARM(2))*(K11 + K12*SQRT(VSQR) + K13*VSQR)
      DYDX(3) = Y(4)
      DYDX(4) = (RPARM(11)/RPARM(2))*(K21 + K22*SQRT(VSQR) + K23*VSQR)

      RETURN
      END



      SUBROUTINE ADPTRK(F,X,XOUT,Y4,Y3,DYDX,NEQ,IPARM,RPARM,ERR)
************************************************************************
*     This subroutine makes subroutine RK34 a semiadaptive differential*
*     equation solver.  Semi adaptive because it checks only to see if *
*     it can reduce the step size, not make it larger.  This is useful *
*     if the calling program needs to know the solution on specific    *
*     intervals.                                                       *
*                                                                      *
*     Since the a posteriori error of RK34 is for the RK3 solutions,   *
*     the use of RK4 values by the user insures that the true error    *
*     is probably less than prescribed by ERR.                         *
*     THIS IS A VERY CONSERVATIVE ROUTINE                              *
*                                                                      *
*     INPUT:                                                           *
*           X     - (REAL) Initial independent variable position       *
*           XOUT  - (REAL) Independent variable position for which     *
*                   solution is desired                                *
*           Y4    - (REAL ARRAY) Initial conditions for dependent vars *
*           NEQ   - (INTEGER) Number of equations in the system        *
*           IPARM - (INTEGER ARRAY) Used to pass constants to the      *
*                   subroutine calculating the system derivatives      *
*           RPARM - (REAL ARRAY) Used to pass constants to the         *
*                   subroutine calculating the sytem derivatives       *
*           ERR   - (REAL) Prescribed error tolerance.  This will be   *
*                   used to decide whether to initiate adaptive proces *        *
*     OUTPUT:                                                          *
*           Y4    - RK4 solutions                                      *
*           Y3    - RK3 solutions for whatever they're worth (ie. RK4  *
*                   values will be more accurate                       *
************************************************************************
      INTEGER NEQ,IPARM(*),I,N
      REAL X,XOUT,Y4(NEQ),Y3(NEQ),DYDX(NEQ),RPARM(*),ERR,APERR,H,XEND,
     .     TEMPY(10)
      EXTERNAL F

C     Ensure that the user input error tolerance is positive
      ERR = ABS(ERR)

C     Save the initial conditions otherwise if APERR is .GT. ERR the
C     initial conditions end up being run through twice
      DO 10 I = 1, NEQ
         TEMPY(I) = Y4(I)
10    CONTINUE

C     Get an initial a posteriori error estimate (APERR)
      CALL RK34 (F,X,XOUT,TEMPY,Y3,DYDX,NEQ,IPARM,RPARM,APERR)

C     Decide whether error is too large and find new step size if it is
      IF (APERR .GT. ERR) THEN
         H = (XOUT - X)/FLOAT(INT(SQRT(SQRT(APERR/ERR)))+1)
         N = INT(SQRT(SQRT(APERR/ERR)))+1

         DO 20 I = 1,N
            XEND = X+ H
            CALL RK34(F,X,XEND,Y4,Y3,DYDX,NEQ,IPARM,RPARM,APERR)
            X = XEND
20       CONTINUE
      ELSE
         DO 30 I = 1,NEQ
            Y4(I) = TEMPY(I)
30       CONTINUE
      END IF

      RETURN
      END


      SUBROUTINE RK34( F, X, XOUT, Y4, Y3, DYDX, NEQ, IPARM, RPARM, ERR)
************************************************************************
*     RK4   - PERFORMS THIRD AND FOURTH ORDER RUNGE-KUTTA INTEGRATION  *
*             OF A SYSTEM OF NEQ FIRST ORDER EQUATIONS OF THE FORM:    *
*                            DY/DX = F(X,Y)                            *
*             Returns an a posteriori error estimate (RK4 - RK3) for   *
*             RK3 for use in adaptive techniques.                      *
*     INPUT VARIABLES                                                  *
*             Form: SUBROUTINE F( X, Y, IPARM, RPARM, DYDX )           *
*             The subroutine F evaluates the NEQ components of the     *
*             system of diferential equuations DY/DX = F(X,Y).         *
*             IPARM and RPARM are INTEGER and REAL arrays used to      *
*             communicate any information between calling program and F*
*             NOTE: THE SUBROUTINE NAME OF THE SUBROUTINE F MUST       *
*             BE DECLARED "EXTERNAL" IN THE CALLING PROGRAM.           *
*     X     - (REAL) The starting value of the independent variable    *
*     XOUT  - (REAL) The point at which the solution is desired        *
*     Y4    - (ARRAY of length NEQ, REAL) The array of initial condit  *
*     NEQ   - (INTEGER) Number of equations in the system              *
*     K     - (ARRAY of Dimension (NEQ,4), REAL) Used to store the     *
*             values of K1,K2,K3,K4,K5 for each diff eq in the system  *
*             K1,K2,K3 are used for RK3 & K1,K2,K4,K5 are used for RK4 *
*     IPARM - (ARRAY of unknown length, INTEGER) Used to input         *
*             parameters to subroutine F                               *
*     RPARM - (ARRAY of unknown length, REAL) Used to input parameters *
*             to subroutine F                                          *
************************************************************************
*     OUTPUT VARIABLES                                                 *
*                                                                      *
*     XOUT  - Same as input                                            *
*     Y     - Solution array at XOUT                                   *
*     DYDX  - Derivative array AT XOUT                                 *
************************************************************************
      INTEGER J, NEQ, IPARM(*)
      REAL X,XOUT,Y4(NEQ),DYDX(NEQ),H,RPARM(*),H2,TEMPX,TEMPY(10),
     .     K(10,5),Y3(NEQ),ERR

C     Calculate the step size
      H = XOUT - X

C     Calculate H/2
      H2 = .5*H

C     Calculate the K values for the functions K(Jth eqn,#of K)
c     Find K1
      CALL F(X,Y4,IPARM,RPARM,DYDX)
      DO 10 J = 1, NEQ
         K(J,1) = DYDX(J)
10    CONTINUE

c     Find K2
      DO 15 J = 1, NEQ
         TEMPY(J) = Y4(J) + H2*K(J,1)
15    CONTINUE
      TEMPX = X + H2
      CALL F(TEMPX,TEMPY,IPARM,RPARM,DYDX)
      DO 20 J = 1, NEQ
         K(J,2) = DYDX(J)
20    CONTINUE

c     Find K3
      DO 25 J = 1, NEQ
         TEMPY(J) = Y4(J) + .75*H*K(J,2)
25    CONTINUE
      TEMPX = X + .75*H
      CALL F(TEMPX,TEMPY,IPARM,RPARM,DYDX)
      DO 30 J = 1, NEQ
         K(J,3) = DYDX(J)
30    CONTINUE

c     Find K4
      DO 35 J = 1, NEQ
         TEMPY(J) = Y4(J) + H2*K(J,2)
35    CONTINUE
      TEMPX = X + H2
      CALL F(TEMPX,TEMPY,IPARM,RPARM,DYDX)
      DO 40 J = 1, NEQ
         K(J,4) = DYDX(J)
40    CONTINUE

c     Find K5
      DO 45 J = 1, NEQ
         TEMPY(J) = Y4(J) + H*K(J,4)
45    CONTINUE
      TEMPX = X + H
      CALL F(TEMPX,TEMPY,IPARM,RPARM,DYDX)
      DO 50 J = 1, NEQ
         K(J,5) = DYDX(J)
50    CONTINUE

C     Calculate the function values at XOUT by RK3
      DO 60 J = 1, NEQ
         Y3(J) = Y4(J) + H*(2.*K(J,1) + 3.*K(J,2) + 4.*K(J,3))/9.
60    CONTINUE

C     Calculate the function values at XOUT by RK4
      DO 70 J = 1, NEQ
         Y4(J) = Y4(J) + H*(K(J,1) + 2.*K(J,2) + 2.*K(J,4)+ K(J,5))/6.
70    CONTINUE

C     Calculate the MAX a posteriori error estimate for RK3
C     Since Y values for RK4 are output, this is the worst case error
      ERR = 0.
      DO 80 J = 1, NEQ
         TEMPX = ABS(Y4(J) - Y3(J))
         IF (TEMPX .GT. ERR) THEN
            ERR = TEMPX
         END IF
80    CONTINUE

      RETURN
      END
