C****************************
C* A I R C R A F T  T R I M * 
C*            by            *
C*      Mark E. Peters      *
C*        April 1994        *
C****************************
C
C    V A R I A B L E   D E S C R I P T I O N  
C  *******************************************

C
C       R E Q U I R E D    V A L U E S 
C    ***********************************
C  
C   REQUIRED WING INPUT 
C  *********************
C  SW ............WING AREA (ft^2) 
C  ARW............ASPECT RATIO OF THE WING
C  EW.............OSWALD/SPAN EFFICIENCY
C  AOW............2-D LIFT CURVE SLOPE (PER DEG)     
C  ALO............ALPHA AT WHICH LIFT IS ZERO FOR WING
C  CDO............DRAG COEFFICENT 
C  CMW............WING MOMENT COEFFICIENT
C  IW ............WING INCIDENCE IN DEGREES
C  CW.............MEAN AERODYNAMIC CHORD FOR THE WING (ft)
C  LBA............LAMBDA OF SMETANA. (TIP CHORD/ROOT CHORD)
C  WSPAN..........WING SPAN (ft)  (B IN SMETANAS CODE)
C   
C   REQUIRED TAIL INPUT
C  *********************      
C  ST ............TAIL AREA (ft^2)
C  ART............ASPECT RATIO FOR THE TAIL
C  ET.............OSWALD EFFICIENCY FOR THE TAIL
C  AOT............2-D LIFT CURVE SLOPE
C  ALOT...........ALPHA AT WHICH LIFT IS ZERO FOR THE WING
C  CDOT...........DRAG COEFFICIENT FOR TAIL (PARASITE)
C  CMOT...........MOMENT COEFFICIENT FOR TAIL
C  IT ............INCIDENCE OF TAIL (DEGREES)
C  CT.............MEAN AERODYNAMIC CHORD DATA FOR TAIL (ft)
C  CFC............FLAP TO CHORD RATIO
C  LBAT...........TAIL TAPER RATIO
C  TSPAN..........TAIL SPAN (SMETANA USES BT) 
C  SELEV..........AREA OF ELEVATOR (ft^2) 
  


C  DISTANCES AND FUSELAGE INPUT 
C *******************************
C                 ALL MEASUREMENTS ARE IN (ft)  
C                 THESE ARE DEFINED IN SMETANAS TEXT 
C                 
C  XA.............. DEFINED ON FIG 9 (PG 52 of Smetana's text)
C  ZA.............. DEFINED ON FIG 9 (PG 52 of Smetana's text)
C  LT.............. DEFINED ON FIG 9 (PG 52 of Smetana's text)
C  HT.............. DEFINED ON FIG 9 (PG 52 of Smetana's text)
C  ZJ.............. THE PERPINDICULAR DISTANCE FROM THE CG TO THE 
C                   THRUST LINE.
C  LB...............FUSELAGE LENGTH
C  WFUS.............FUSELAGE MAX WIDTH
C  XFUS......... ...DISTANCE OF NOSE TO WING C/4  
C  ZT. .............DIST. OF THRUST LINE BELOW C.G.
 
 
C  CDF............FUSELAGE DRAG COEFICIENT.  NONDIMENSIONALIZED 
C                 WITH WING AREA.
C  IYY............. Iyy Moment of inertia 



C   C A L C U L A T E D   V A L U E S 
C *************************************

C    WING
C **********
C  AW..............3-D LIFT CURVE SLOPE FOR WING
C  CLO.............LIFT COEF.  AT ZERO ANGLE OF ATTACK
C  CDIW............INDUCED DRAG COEF FOR WING
C  WL..............WING LIFT (lbs) 
C  DW..............WING DRAG (lbs)  
   

C   TAIL
C **********
C  AT..............3-D LIFT CURVE SLOPE FOR TAIL
C  CLOT............LIFT COEF. AT ZERO ANGLE OF ATTACK FOR TAIL
C  CDIT............INDUCED DRAG COEF. FROM TAIL
C  E...............DOWNWASH ANGLE FROM WING (deg) 
C  CMDE............MOMENT COEFFICIENT FROM ELEVATOR DEFLECTION
C  CLDE............LIFT COEFICIENT FROM ELEVATOR DEFLECTION
C  DE..............ELEVATOR DEFLECTION (deg) 
C  TL..............TAIL LIFT (lbs) 
C  DT..............TAIL DRAG (lbs) 
C  MT..............TAIL MOMENT 
 
C   O T H E R    V A R I A B L E S
C**********************************



C S M E T A N A S   V A R I A B L E S
C************************************
C  (SEE SMETANA'S TEXT)

C  THRUST.......... THRUST FOR LOPROG DATA FILE
C                   THRUST LINE.
C  GCG..............GCOS(GAMMA)
C  GSG..............GSIN(GAMMA
C  CXZ..............COS(XZ)
C  SXZ..............SIN(XZ)
C  CL2DW............CL2DW 
C  CDPIE............PARASITE DRAG COEFFICIENT


C   OTHER VARIABLES
C ******************************************************

C  V............... AIRCRAFT VELOCITY (ft/sec)
C  P............... AIR DENSITY (slug/ft^3)
C  ALOW.........   THIS IS THE LOW BOUND FOR ALPHA DURING ITERATION 
C  AHIGH.........  THIS IS THE HIGH BOUND FOR ALPHA DURING ITERATION 
C  D............   DRAG OF AIRPLANE
C  CD ............ OVERALL DRAG COEFFICIENT
C  CL ............ OVERALL LIFT COEFFICIENT



C 
C   D E C L A R E   V A R I A B L E S
C**************************************

       REAL SW,ARW,EW,AOW,ALO,CDO,CMW,IW,CW,ST,ART,ET
       REAL AOT,ALOT,CDOT,CMOT,IT,CT,CFC,XCG,ZCG,XWING,ZWING
       REAL XTAIL,ZTAIL,CDF,AW,CLO,CDIW,AT,CLOT,CDIT,E,CMDE
       REAL CLDE,DE,XA,ZA,LT,HT,V,WEIGHT,Q,WL,MW,TAIL,AA,BB,CC
       REAL DD,EE,FF,GG,QA,QB,QC,DE,TL,PI,P,L,ALOW,AHIGH 
       REAL DW,DT,FD,D,CD,CL,TMP,LBA  
       REAL MT,MCG,VLOW,VHIGH 
       REAL THRUST,ZJ,GCG,GSG,CXZ,SXZ,WSPAN,TSPAN,EC,LBAT,CL2DW
       REAL CDPIE,LB,WFUS,XFUS,ZT,SELEV,IYY,VLOW,VHIGH 
       CHARACTER * 10 FILENAME 
       CHARACTER * 10 DATAFILE
       CHARACTER * 10 VECTOR
       CHARACTER * 10 MAT
       CHARACTER * 1 CHOICE1
       CHARACTER * 1 CHOICE2
       CHARACTER * 1 CHOICE3
     
C ****************************************************************
C  I N T E R A C T I V E    P A R T   O F   P R O G R A M
C ****************************************************************
      
        print *, '                     ********************************'
        print *, '                     *    AIRCRAFT TRIM ANALYSIS    *'      
        print *, '                     *         CALCULATIONS         *'
        print *, '                     ********************************'



        PRINT *, ' '
        PRINT *, ' '
        PRINT *, '  ENTER THE NAME OF THE DATA FILE'
        READ  *, FILENAME


C        CHOOSE A VELOCITY AND WEIGHT
C *************************************


           PRINT *,'ENTER THE TRIM WEIGHT (lbs) '
           READ*, WEIGHT





        print *, '                     ********************************'
        print *, '                     *    AIRCRAFT TRIM ANALYSIS    *'
        print *, '                     *            OPTIONS           *'
        print *, '                     ********************************'

        
        print *, 'DO YOU WANT TO CREATE A SMETANA LOPROG DATA FILE?'
        print *, '(y/n)'
        read  *, CHOICE1
        print *, ' DO YOU WANT TO CREATE A MCG vs ALPHA VECTOR?'
        print *, '(y/n)'
        read  *, CHOICE2    
        print *, ' DO YOU WANT A SINGLE POINT CALCULATION OR DO '
        print *, ' YOU WANT TO ITERATE OVER A RANGE OF VELOCITIES?'
        print *, '(s/r)'
        read  *, CHOICE3    
        
        if (CHOICE3 .EQ. 'r') then 
            print *, 'ENTER RANGE (VLOW VHIGH)'
            read  *, VLOW,VHIGH
             V = VLOW
          else

           PRINT *,'ENTER THE TRIM VELOCITY (FT/SEC)'
           READ *, V

            
        endif  
       


C    END THE INTERACTIVE ROUTINE
C****************************************




C        I N P U T   D A T A 
C**************************************
C  
       OPEN(1, FILE=FILENAME)

C          REQUIRED WING DATA
C          ******************
           READ(1,*) SW
           READ(1,*) ARW
           READ(1,*) EW
           READ(1,*) AOW
           READ(1,*) ALO 
           READ(1,*) CDO 
           READ(1,*) CMW
           READ(1,*) IW
           READ(1,*) CW
           READ(1,*) LBA
           READ(1,*) WSPAN
          
C          REQUIRED TAIL DATA
C          ****************** 
           READ(1,*) ST
           READ(1,*) ART 
           READ(1,*) ET
           READ(1,*) AOT
           READ(1,*) ALOT
           READ(1,*) CDOT
           READ(1,*) CMOT
           READ(1,*) IT
           READ(1,*) CT
           READ(1,*) CFC
           READ(1,*) LBAT
           READ(1,*) TSPAN
           READ(1,*) SELEV


C          DISTANCES AND FUSELAGE DATA
C          ***************************
           READ(1,*) XA
           READ(1,*) ZA
           READ(1,*) LT 
           READ(1,*) HT
           READ(1,*) ZJ
           READ(1,*) LB
           READ(1,*) WFUS
           READ(1,*) XFUS
           READ(1,*) ZT
           READ(1,*) CDF
           READ(1,*) IYY
           CLOSE(1)


C   B E G I N    C A L C U L A T I O N S
C******************************************

      PI = 3.14159



 
C WING CALCULATIONS
C************************************************


C  AW..............3-D LIFT CURVE SLOPE FOR WING

             AW = AOW/(1+(57.3*AOW/(PI*EW*ARW)))


C  CLO.............LIFT AT ZERO ANGLE OF ATTACK
 
               CLO = (0-ALO)*AW



   
C TAIL CALCULATIONS
C************************************************  

C  AT..............3-D LIFT CURVE SLOPE FOR TAIL
  
             AT = AOT/(1+(57.3*AOT/(PI*ET*ART)))



C  CLOT............LIFT AT ZERO ANGLE OF ATTACK FOR TAIL

            CLOT =(0-ALOT)*AT


 
C  CMDE............MOMENT COEFFICIENT FROM ELEVATOR DEFLECTION
 

           CMDE=-(AT/PI)*((1-CFC)*(CFC*(1-CFC))**0.5)    


   
C  CLDE............LIFT COEFICIENT FROM ELEVATOR DEFLECTION


           CLDE = (AT/PI)*(acos(1-2*CFC)+2*(CFC*(1-CFC))**0.5) 



    
 
C********************************************************
C          SET VALUES BEFORE ROUTINE BEGINS
C********************************************************

           ALOW = -20.0
           P = 2.3E-03
           A = 8.0
 

C AT THIS POINT IN THE CODE THE PROGRAM HAS CALCULATED
C ALL THE NECESSARY 3-D COEFFICIENTS AND IS READY TO 
C BALANCE THE MOMENTS TO FIND TRIM.  THE CODE CHOOSES 
C AN ANGLE OF ATTACK AND SOLVES FOR THE REQUIRED
C ELEVATOR DEFLECTION. IT THEN COMPARES THE LIFT PRODUCED
C TO THE WEIGHT.  IT ITERATES OVER ANGLE OF ATTACK UNTIL
C THE LIFT EQUALS THE WEIGHT. THIS IS THE FINAL TRIM CONDITION.
C THE CODE OUTPUT WILL CONTAIN THE FINAL ELEVATOR DEFLECTION 
C AND THE FINAL ANGLE OF ATTACK. 



C START ITERATION WITH ANGLE OF ATTACK
C**********************************

 10           Q =  0.5*P*V**2




C CALCULATE THE WING LIFT, DRAG AND MOMENT
C************************************

           WL = Q*SW*(AW*A +CLO)

           MW = Q*SW*CW*CMW
  
           CDIW =(AW*A +CLO)**2/(PI*EW*ARW) 

           DW = Q*SW*(CDO + CDIW)



C CALCULATE DOWNWASH ANGLE 
C******************************
C 
       TMP = (((1/LBA)**0.3)/(ARW**0.725))*((3*CW)/(XA+LT))**0.25  
       E = 20*(AW*A +CLO)*TMP                      
       OLDE = 35*(CLO+AW*(A+IW))/ARW 
 

C CALCUATE THE FINAL MOMENT. (Mcg) ASSUMING 
C THAT THIS VALUE IS ZERO, WHAT ELEVATOR TRIM 
C WILL BE NECESSARY TO MAKE THE BALANCE?

C THIS MOMENT EQUATION IS TAKEN DIRECTLY FROM SMETANAS BOOK
C ON PAGE 52 

C IT IS AS FOLLOWS: Mcg = Lxa + Dza + Ma.c. -Ltlt +Dtht + Ma.c.t 
C THE CODE ASSUMES THAT THE THREE WING TERMS WILL BE BALANCED BY
C THE REMAINING TAIL TERMS. HENCE THE TAIL VARIABLE. 
C THE THRUST IS OMMITED TO SIMPLIFY THE EQUATIONS.  

          TAIL = -(WL*XA + DW*ZA + MW)

       
C ALL THREE TAIL EXPRESSIONS CONTAIN BOTH ANGLE OF ATTACK AND 
C ELEVATOR DEFLECTION ANGLE.  THE INCLUSION OF INDUCED DRAG MAKES 
C THE SOLUTION QUADRATIC.

C THEREFORE:

C THE FOLLOWING TEMPORARY VARIABLES WILL
C BE USED TO SIMPLIFY THE EXPRESSION WHICH
C SOLVES FOR DE


        AA = Q*ST*LT*(CLOT + AT*(A+IT-E))

        BB = Q*ST*LT*CLDE

        CC = Q*ST*CT*CMT

        DD = Q*ST*CT*CMDE

        EE = Q*ST*HT*CDOT

        FF = AA/(Q*ST*LT)     
      
        GG = CLDE


        QA = GG**2/(PI*ET*ART) 

     
        QB = 2*FF*GG/(PI*ET*ART) +DD - BB


        QC = FF**2/(PI*ET*ART) -AA +EE -TAIL 



C    ELEVATOR DEFLECTION

        DE = (-QB - (QB**2 -4*QA*QC)**.5)/(2*QA) 

 

C     TAIL LIFT 

        TL = (AA+BB*DE)/LT  

C    LIFT OF THE TOTAL AIRPLANE

        L = TL + WL   




C THE FOLLOWING ROUTINE WILL CHOOSE THE NEXT 
C ANGLE OF ATTACK TO BE TRIED


        IF (ABS(WEIGHT-L) .GT. .1) THEN 

         IF(L .GT. WEIGHT) THEN

            AHIGH = A
            A = A - (ABS(ALOW - AHIGH))/2 


            GOTO 10 


         ELSE  

            ALOW  = A
            A = A + (ABS(ALOW - AHIGH))/2 
            GOTO 10  

        ENDIF          
        ENDIF       


C ONCE THE CODE REACHES THIS POINT THE ITERATION IS FINISHED AND 
C ANGLE OF ATTACK AND ELEVATOR DEFLECTION HAVE BEEN SOLVED. 


C  F I N A L   C O M P U T A T I O N S 




C      CDIT............INDUCED DRAG COEF. FROM TAIL

       CDIT = (AT*(A+IT-E) +CLOT)**2/(PI*ET*ART) 

C  DT..............TAIL DRAG (lbs) 

       DT = Q*ST*(CDOT+CDIT)

C  FD.........  FUSELAGE DRAG 

       FD = Q*SW*CDF
       
C      TOTAL DRAG

       D = DT + DW + FD


C      OVERALL LIFT AND DRAG COEFFICIENTS

       CD = D/(Q*SW)
C* *****************************************************
C        PRINT PERTINANT DATA TO THE SCREEN  
C******************************************************
       CL = L/(Q*SW)





        print *, '                     ********************************'
        print *, '                     *    AIRCRAFT TRIM ANALYSIS    *'      
        print *, '                     *         CALCULATIONS         *'
        print *, '                     ********************************'

       print *, ' '
       print *, ' '
       print *, ' '
       print *, ' '
       print *, ' D A T A   F R O M   I N P U T  F I L E'   
       print *, '*****************************************************'

      
        PRINT*,'W I N G   D A T A'
        PRINT*,'*****************'                  
        PRINT*,'WING AREA.................................',SW,'ft^2' 
        PRINT*,'ASPECT RATIO OF THE WING..................',ARW
        PRINT*,'OSWALD EFFICIENCY.........................',EW  
        PRINT*,'2-D LIFT CURVE SLOPE......................',AOW  
        PRINT*,'ALPHA AT WHICH LIFT IS ZERO FOR WING......',ALO,'deg' 
        PRINT*,'DRAG COEFFICENT...........................',CDO
        PRINT*,'WING MOMENT COEFFICIENT...................',CMW
        PRINT*,'WING INCIDENCE IN DEGREES ................',IW
        PRINT*,'MEAN AERODYNAMIC CHORD FOR THE WING.......',CW,'ft'
        PRINT*,'LAMBDA OF SMETANA. (TIP CHORD/ROOT CHORD).',LBA
        PRINT*,'WING SPAN (ft) ...........................',WSPAN
        PRINT*,' '


        PRINT*,'T A I L   D A T A'
        PRINT*,'*****************'
        PRINT*,'TAIL AREA.................................',ST,'ft^2'
        PRINT*,'ASPECT RATIO FOR THE TAIL.................',ART
        PRINT*,'OSWALD EFFICIENCY FOR THE TAIL............',ET
        PRINT*,'2-D LIFT CURVE SLOPE......................',AOT
        PRINT*,'ALPHA AT WHICH LIFT IS ZERO FOR THE WING..',ALOT,'deg' 
        PRINT*,'DRAG COEFFICIENT FOR TAIL (PARASITE)......',CDOT
        PRINT*,'MOMENT COEFFICIENT FOR TAIL...............',CMOT
        PRINT*,'INCIDENCE OF TAIL (DEGREES)...............',IT,'deg'
        PRINT*,'MEAN AERODYNAMIC CHORD FOR TAIL...........',CT,'ft'
        PRINT*,'FLAP TO CHORD RATIO.......................',CFC 
        PRINT*,'TAIL TAPER RATIO..........................',LBAT
        PRINT*,'TAIL SPAN (SMETANA USES BT)...............',TSPAN
        PRINT*,'AREA OF ELEVATOR (ft^2) ..................',SELEV 
        PRINT*,' '


        PRINT*,'DISTANCES AND FUSELAGE DATA'
        PRINT*,'*************************'
        PRINT*,'XA.......................',XA 
        PRINT*,'ZA.......................',ZA
        PRINT*,'LT.......................',LT
        PRINT*,'HT.......................',HT
        PRINT*,'ZJ.......................',ZJ
        PRINT*,'LB.......................',LB
        PRINT*,'WFUS.....................',WFUS
        PRINT*,'XFUS.....................',XFUS
        PRINT*,'ZT.......................',ZT
        PRINT*,'FUSELAGE DRAG COEFFICIENT.',CDF
        PRINT *,'  '




       PRINT *,'T R I M   C O N D I T I O N   I N P U T S'
       PRINT *,'******************************************'
       PRINT *,'AIRCRAFT VELOCITY......................',V,'ft/sec'
       PRINT *,'AIR DENSITY.........................',P,'slug/ft^3'
       PRINT *,'AIRCRAFT WEIGHT.....................',WEIGHT,'lbs'
       PRINT *,' '
        

       PRINT *,' '
       PRINT *,' ' 
       PRINT *,'C A L C U L A T E D   V A L U E S'
       PRINT *,'*********************************'
       PRINT *,' '
       PRINT *, 'AIRCRAFT ANGLE OF ATTACK ...',A,'deg'
       PRINT *, 'TOTAL LIFT PRODUCED ........',L,'lbs'
       PRINT *, 'TOTAL DRAG PRODUCED .........',D,'lbs'
       PRINT *, 'OVERALL LIFT COEFFICIENT.....',CL 
       PRINT *, 'OVERALL DRAG COEFFICIENT.....',CD
       PRINT *, ' ' 
       PRINT *,'S M E T A N A S   G E O M E T R Y'
       PRINT *,'**********************************'
       PRINT *,'XA..............',XA,'ft'
       PRINT *,'LT..............',LT,'ft'
       PRINT *,'ZA..............',ZA,'ft'
       PRINT *,'HT..............',HT,'ft'
       PRINT *,' '
       PRINT *,'W I N G  D A T A'
       PRINT *,'********************************'
       PRINT *,'3-D LIFT CURVE SLOPE FOR WING.........',AW
       PRINT *,'LIFT COEF.  AT ZERO ANGLE OF ATTACK...',CLO
       PRINT *,'INDUCED DRAG COEF FOR WING............',CDIW
       PRINT *,'WING LIFT (lbs).......................',WL,'lbs'
       PRINT *,'WING DRAG (lbs).......................',DW,'lbs'
       PRINT *,' ' 
       PRINT *,'T A I L  D A T A'
       PRINT *,'*********************************'        
      PRINT *,'3-D LIFT CURVE SLOPE FOR TAIL.........',AT
      PRINT *,'LIFT COEF. AT ZERO ANGLE OF ATTACK....',CLOT 
      PRINT *,'INDUCED DRAG COEF. FROM TAIL.........',CDIT
      PRINT *,'DOWNWASH ANGLE FROM WING.............',E,'deg' 
      PRINT *,'MOMENT COEF. ELEVATOR DEFLECTION.....',CMDE
      PRINT *,'LIFT COEF. FROM ELEVATOR DEFLECTION...',CLDE 
      PRINT *,'ELEVATOR DEFLECTION...................',DE,'deg'
      PRINT *,'TAIL LIFT............................. ',TL,'lbs'
      PRINT *,'TAIL DRAG.............................',DT,'lbs' 
      PRINT *,' '
 

       
       PRINT *,'OLDE.................',OLDE,'deg' 


      



C******************************************************
C             SMETANA DATA FILE FOR LOPROG 
C******************************************************

C        THIS IS CHOICE1


         IF (CHOICE1 .EQ. 'y') THEN 


          print *, 'ENTER THE NAME FOR SMETANAS DATA FILE'
          read  *, DATAFILE


        OPEN(UNIT=1, FILE=DATAFILE)
  
          WRITE(UNIT=1, *) '0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0'    

          WRITE(UNIT=1, *) P

          WRITE(UNIT=1, *) V
  
          WRITE(UNIT=1, *) WEIGHT/32.2 

          WRITE(UNIT=1, *) IYY

C         CALCULATE THE REQUIRED THRUST
C         THRUST = DRAG   
C         THRUST IS ASSUMED TO BE PARALLEL WITH FUSE
C         REFERENCE LINE 
          THRUST = D/COS((A/57.3))

          WRITE(UNIT=1, *) THRUST 

          WRITE(UNIT=1, *) ZJ

C         THIS PROGRAM DOES ONLY LEVEL FLIGHT 
C         THEREFORE:
          GCG = 32.2
          GSG = 0.0

          WRITE(UNIT=1, *) GCG

          WRITE(UNIT=1, *) GSG

C         THESE VARIABLES DESCRIBE THE ANGLE BETWEEN THE 
C         THRUST AXIS AND THE RELATIVE WIND.  THEREFORE
C         WITH THIS PROGRAM AND ITS PREVIOUS THRUST ASSUMPTION

          CXZ  = COS(A/57.3)   
          SXZ  = SIN(A/57.3)


          WRITE(UNIT=1, *) CXZ

          WRITE(UNIT=1, *) SXZ

          WRITE(UNIT=1, *) SW

          WRITE(UNIT=1, *) ST

          WRITE(UNIT=1, *) WSPAN

          WRITE(UNIT=1, *) TSPAN

          WRITE(UNIT=1, *) EW  

          WRITE(UNIT=1, *) EC

          WRITE(UNIT=1, *) LBA

          WRITE(UNIT=1, *) LBAT 

C         SMETANAS CL2DW:  THIS VARIABLE IS THE 2-D LIFT 
C         COEFFICIENT OF THE WING.

          CL2DW = (0-ALO)*AOW  + (A+IW)*AOW 

          WRITE(UNIT=1, *) CL2DW

          WRITE(UNIT=1, *) IT
          
          WRITE(UNIT=1, *) AW

          WRITE(UNIT=1, *) AT


C         PARASITE DRAG COEFFICIENT 
C         
          CDPIE = CDO+ CDOT*(ST/SW) + CDF

          WRITE(UNIT=1, *) CDPIE


C         2-D CDA IS ALWAYS ZERO.

          WRITE(UNIT=1, *) 0.0    

C         SMETANA DEFINES ANGLE OF ATTACK AS THE ANGLE
C         OF ATTACK AS THE WING.
     
          WRITE(UNIT=1, *) A+IW

          WRITE(UNIT=1, *) IW

          WRITE(UNIT=1, *) LT

          WRITE(UNIT=1, *) LT+XA
 
          WRITE(UNIT=1, *) XA  

          WRITE(UNIT=1, *) ZA
      
          WRITE(UNIT=1, *) LB

          WRITE(UNIT=1, *) WFUS

          WRITE(UNIT=1, *) XFUS

          WRITE(UNIT=1, *) -XA

          WRITE(UNIT=1, *) ZT

          WRITE(UNIT=1, *) SELEV

          WRITE(UNIT=1, *) 1,1,0.0001
         
          WRITE(UNIT=1, *) 0,0,0


            CLOSE(1)     

             ENDIF  


C************************************************************
C  MOMENT COEFFICIENT VARIATIONS  
C************************************************************

C  THIS IS CHOICE NUMBER 2 

         IF(CHOICE2 .EQ. 'y' ) THEN

        PRINT *, 'ENTER FILENAME FOR Mcg VECTOR'
        READ  *, vector

        OPEN(UNIT=1, FILE=vector)


     

C THE PURPOSE OF THIS SECTION IS TO PERTURBATE THE AIRCRAFT 
C FROM EQUILIBRIUM AND CALCULATE THE RESTORING MOMENT.

C************************************************
C       VARY ALPHA ABOUT EQUILIBRIUM POINT
C************************************************


        DO 50 A=-10,10
        
C CALCULATE THE WING LIFT, DRAG AND MOMENT
C************************************

           WL = Q*SW*(AW*A +CLO)

           MW = Q*SW*CW*CMW
  
           CDIW =(AW*A +CLO)**2/(PI*EW*ARW) 

           DW = Q*SW*(CDO + CDIW)

       


C CALCULATE DOWNWASH ANGLE 
C******************************
C 
       TMP = (((1/LBA)**0.3)/(ARW**0.725))*((3*CW)/(XA+LT))**0.25  
       E = 20*(AW*A +CLO)*TMP                      
       OLDE = 35*(CLO+AW*(A+IW))/ARW 






C     TAIL LIFT AND MOMENT  

        AA = Q*ST*LT*(CLOT + AT*(A+IT-E))

        BB = Q*ST*LT*CLDE

        TL = (AA+BB*DE)/LT

        CDIT = (AT*(A+IT-E) +CLOT)**2/(PI*ET*ART) 

        DT = Q*ST*(CDOT+CDIT)

        MT =Q*ST*CT*(CMT + CMDE*DE)

 
        MCG = WL*XA + DW*ZA + MW - TL*LT + DT*HT + MT      
C       MCG = WL*XA  +MW  -TL*LT +MT (No Drag Version) 
 
       write(unit=1,*) MCG,A 

 50   continue


       close(1)


         ENDIF



C *********************************************
C  THIS ROUTINE VARIES VELOCITY AND CREATES THE 
C  FOLLOWING VECTORS
C**********************************************

          IF(CHOICE3 .EQ. 'r') then


         print *,'ENTER A FILENAME FOR THE VECTOR MATRIX'
         READ *,  mat

         open(UNIT=2, FILE=mat)


        DO 150 VV=1,100

        V =VLOW + (VV/100)*(VHIGH-VLOW) 


C   SET VALUES BEFORE ROUTINE BEGINS
           ALOW = -20.0
           P = 2.3E-03
           A = 8.0


C START ITERATION WITH ANGLE OF ATTACK
C**********************************

 100           Q =  0.5*P*V**2




C CALCULATE THE WING LIFT, DRAG AND MOMENT
C************************************

           WL = Q*SW*(AW*A +CLO)

           MW = Q*SW*CW*CMW

           CDIW =(AW*A +CLO)**2/(PI*EW*ARW)

           DW = Q*SW*(CDO + CDIW)



C CALCULATE DOWNWASH ANGLE
C******************************
C
       TMP = (((1/LBA)**0.3)/(ARW**0.725))*((3*CW)/(XA+LT))**0.25
       E = 20*(AW*A +CLO)*TMP
       OLDE = 35*(CLO+AW*(A+IW))/ARW


           TAIL = -(WL*XA + DW*ZA + MW)


        AA = Q*ST*LT*(CLOT + AT*(A+IT-E))

        BB = Q*ST*LT*CLDE

        CC = Q*ST*CT*CMT

        DD = Q*ST*CT*CMDE

        EE = Q*ST*HT*CDOT

        FF = AA/(Q*ST*LT)

        GG = CLDE


        QA = GG**2/(PI*ET*ART)


        QB = 2*FF*GG/(PI*ET*ART) +DD - BB


        QC = FF**2/(PI*ET*ART) -AA +EE -TAIL



C    ELEVATOR DEFLECTION

        DE = (-QB - (QB**2 -4*QA*QC)**.5)/(2*QA)



C     TAIL LIFT

        TL = (AA+BB*DE)/LT

C    LIFT OF THE TOTAL AIRPLANE

        L = TL + WL



C THE FOLLOWING ROUTINE WILL CHOOSE THE NEXT
C ANGLE OF ATTACK TO BE TRIED


        IF (ABS(WEIGHT-L) .GT. .1) THEN

         IF(L .GT. WEIGHT) THEN

            AHIGH = A
            A = A - (ABS(ALOW - AHIGH))/2


            GOTO 100


         ELSE

            ALOW  = A
            A = A + (ABS(ALOW - AHIGH))/2
            GOTO 100

        ENDIF
        ENDIF


C ONCE THE CODE REACHES THIS POINT THE ITERATION IS FINISHED AND
C ANGLE OF ATTACK AND ELEVATOR DEFLECTION HAVE BEEN SOLVED.


C  F I N A L   C O M P U T A T I O N S
C*************************************

       CDIT = (AT*(A+IT-E) +CLOT)**2/(PI*ET*ART)

       DT = Q*ST*(CDOT+CDIT)

       FD = Q*SW*CDF

       D = DT + DW + FD

       CD = D/(Q*SW)

       CL = L/(Q*SW)


          write(unit=2,*) V,DE,A,D 


 150    CONTINUE

          CLOSE(2)


         ENDIF

           END
