      PROGRAM MASTER

c     file JKayVLM.f - last original verison

c     a program to estimate the subsonic aerodynamics and stability
c     and control derivatives for an assessment of the control power
c     required in conceptual aircraft designs

c     the basic method is the vortex ring variation of the vlm method

c     conceived and written by Jacob Kay
c     VPI Aerospace and Ocean Engineering Department, 1992

c     mods by W.H. Mason and Alex Benoliel, March 1993
c     last mod, May 15, 1993
c     mod to correct small bug, Nov. 16, 1996 by W.H. Mason
c	mods by L.T Leifsson, Feb. 2004:
c		- program writes to screen and output file a header
c		  with information about the program and mods
c		- program pauses when run is finished, asks to press
c		  enter to close program
c		

      REAL  VCref(5,10,6,3),DCL(5,2),    DCM(5,2),      QBAR,        
     +      DELslat(5),     DELflap(5),  VC(5,10,6,3),  GAMMA(5,8,10),
     +      FSECT,          Xcg,         Ycg,           Zcg,
     +      LIFT0,          PITCH0,      PI,            LIFT,
     +      PITCH,          ROLL,        YAW,           p,
     +      q,              V,           S,             C,
     +      ALPHA,          LCSDW,       DRAG,          YAW0,
     +      CLalpha,        CmCL,        CLq,           Cmq,
     +      CROLLp,         B,           ASYM,          ASYMDEFL,
     +      DCroll(5,2),    CLATAIL1,    CLATAIL2,      CYAWp,
     +      Cdi,            CdiCL2,      CL0,           Cm0,
     +      Cmalpha,        cref,        Xcgref,        CL,
     +      CH,             Cybeta,      Cnbeta,        Clbeta,
     +      Cyr,            Cnr,         Clr,           BETA,
     +      DCy(5,2),       DCn(5,2),    LIFTTAIL,      LIFTAIL0,
     +      MACH,           bplon(5,4,3),bplat(5,4,3),  slatlon(5),
     +      slatlat(5),     flaplon(5),  flaplat(5),    GAMMAref(5,8,10) 

      INTEGER I,      J,      K,     LAT,   IW1,  IW2,  IWING1,  IWING2,
     +        ITAIL1, ITAIL2, ICASE, ISECT, jsect
      character*1     input,lotus
      character*12    loninf, latinf,in
      character*8     vlmout


c     ===================================================================
c     This master program calls CENTRAL to determine the stability
c     and control derivatives of a design configuration.
c
c     Prepared by Jacob Kay,    Sept. 1992.
c     ===================================================================

      PI   = 4.0*ATAN(1.0)
      V    = 1.0
      QBAR = .5*V*V
c
c     *** Get basic geometric information ***

      write(6,3000)
	write(6,3001)
	write(6,3002)
	write(6,3003)
      write(*,*) 'Enter output file name (max. of 8 characters):'
      read (*,'(a8)') vlmout
      open(12,file=vlmout)
C      open(UNIT=12,FILE=vlmout,STATUS='NEW',ACTION='WRITE')
      rewind(12)
      write(*,*)  'Input from file (y/n)?'
      read(*,'(a1)') input
      if(input.eq.'y'.or.input.eq.'Y') then
	  write(*,*)'Enter file name:'
	  read(*,1100)in
	  write(*,*)in
	  open(9,file=in)
        read(9,1005)mach,s,c,b,ch,icase,xcg,zcg
        read(9,1100) loninf
	  if(icase.eq.1.or.icase.eq.4)then
	    read(9,*)
        else
          read(9,1100) latinf
        endif
      else
        write(6,*)  'ENTER MACH NUMBER: (0 - 0.8)'
        READ*, MACH

        write(6,*) 'ENTER REFERENCE AREA.'
        READ*,S

        write(6,*) 'ENTER REFERENCE CHORD.'
        READ*,C

        write(6,*) 'ENTER REFERENCE SPAN.'
        READ*,B

        write(6,*) 'ENTER HEIGHT ABOVE GROUND.'
        READ*,CH

c     *** Ask for longitudinal or lateral case ***

        write(6,*) 'ENTER: 1 - Longitudinal derivatives only'
        write(6,*) '       2 - Lateral/directional derivatives only'
        write(6,*) '       3 - Both'
        write(6,*) '       4 - Basic calculations only'
        READ*, ICASE

	    write(*,*)'Enter Longitudinal Geometry File:'
        read(*,1100) loninf
        if(icase.ne.1.and.icase.ne.4)then
	          write(*,*)'Enter Lateral Geometry File:'
              read(*,1100) latinf
              endif
		
        write(6,*) 'ENTER X-cg.'
        READ*,Xcg 

        write(6,*)'ENTER Z-cg (+ UP)'
        READ*, Zcg
      endif
C	LOTUS is an obsolte program!
C      write(*,*) 'Ouput LOTUS 1-2-3 file (y/n)?'
C      read(*,'(a1)') lotus
	lotus = 'n'

c     *** Decide to turn on or off ground effect ***

      IF (CH.GT.5*B) CH = 99999.99999
      
c 50   IF (ICASE.EQ.2) GOTO 60
      if(input.eq.'y'.or.input.eq.'Y') then
        read(9,1010)itail1,itail2,iw1,iw2
      else
        write(6,*)'ENTER 1st SECTION NO. DESIGNATED AS H. TAIL SURFACE.'
        READ*, ITAIL1
        write(6,*)'ENTER 2nd SECTION NO. DESIGNATED AS H. TAIL SURFACE.'
        READ*, ITAIL2
        write(6,*)'ENTER 1st SECTION NO. DESIGNATED AS WING SURFACE.'
        READ*, IW1
        write(6,*)'ENTER 2nd SECTION NO. DESIGNATED AS WING SURFACE.'
        READ*, IW2
      endif
c     *** Write conditions at top of output ***

 60   BETA = SQRT(1.0-MACH*MACH)
      Xcg = Xcg/BETA    
      C = C/BETA
      xcgref = xcg*beta
      cref   = c*beta
      write(6,1400)  icase, Mach, S, cref, b, ch, xcgref, zcg
	write(12,3000)
	write(12,3001)
	write(12,3002)
	write(12,3003)
      write(12,1400) icase, Mach, S, cref, b, ch, xcgref, zcg
 3000 format('  **************************************************'/
     1        '    VT Aerospace Aircraft Design Software Series'/
     2		'  **************************************************'/)
 3001 format('                   Program JKayVLM'/
     1		/'  A program to estimate the subsonic aerodynamics and'/
     2        '  stability and control derivatives for an assessment'/
     3		'  of the control power required in conceptual aircraft'/
     4		'  designs. The basic method is the vortex ring'/
     5		'  variation of the vlm method.'/)
 3002 format('  Conceived and written by Jacob Kay.'/
     1			'  This program modified by'/
     2        '    - W.H. Mason and Alex Benoliel, March 1993'/
     3        '    - W.H. Mason, Nov. 16 1996'/
     4		'    - L.T. Leifsson, Feb 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'/)
 1400 format(/3x,'VLM Based Stability & Control Evaluation Code'/
     1       /3x,'icase               = ', i12,
     2       /3x,'MACH NUMBER         = ', f12.5,
     3       /3x,'Sref                = ', f12.5, 
     4       /3x,'c ref               = ', f12.5, 
     5       /3x,'b ref               = ', f12.5,     
     6       /3x,'HEIGHT ABOVE GROUND = ', f12.5,
     7       /3x,'X-cg                = ', f12.5,
     8       /3x,'Z-cg                = ', f12.5)

c     *** Get geometry information ***

      OPEN(UNIT = 10, FILE = loninf, STATUS = 'OLD')
      REWIND(UNIT = 10)

* In reading lat/dir data, y coordinates are entered in as
*    z coordinates and z coordinates are entered in as y coordinates
*    to conform to the codes coordinate common coordinate system

      if(icase.ne.1.and.icase.ne.4) then
        OPEN(UNIT = 11, FILE = latinf, STATUS = 'OLD')
        REWIND(UNIT = 11)
        read(11,*)
        READ(11,1000) FSECT
        jsect = fsect
        write(6,1530) jsect
        DO 110, K = 1, jSECT
	     write(6,1505)k,k,k
           DO 100, I = 1, 4
              READ(11,1000) BPlat(k,I,1),BPlat(k,I,3),BPlat(k,I,2)
              write(6,1510) i,BPlat(k,I,1),BPlat(k,I,3),BPlat(k,I,2)
	        bplat(k,i,1) = bplat(k,i,1)/beta
 100       CONTINUE
           READ(11,1000) slatlat(K), flaplat(K)
           write(6,1520) k, slatlat(K),k, flaplat(K)
 110    continue
      endif

      read(10,*)
      READ(10,1000) FSECT
      ISECT     = FSECT
      write(6,1500) isect
      DO 130, I = 1, ISECT
         write(6,1505)i,i,i  
         DO 120, J = 1, 4
            READ(10,1000) BPlon(i,J,1),BPlon(i,J,2),BPlon(i,J,3)
            write(6,1510) j,BPlon(i,J,1),BPlon(i,J,2),BPlon(i,J,3) 
	      bplon(i,j,1) = bplon(i,j,1)/beta
 120     CONTINUE
         READ(10,1000) slatlon(I), flaplon(I)
         write(6,1520) i, slatlon(I),i, flaplon(I)
 130  CONTINUE


 1000 FORMAT(3F10.5)
 1005 format(/,f4.1,/,4(f10.2,/),i1,/,f10.3,/,f10.3)
 1010 format(3(i1,/),i1)
 1100 format(a12)
 1500 format(/3x,'reading longitudinal planform data, isect = ', i2)
 1505 format(/,3x,' Point No.     x(',i1,')',8x,'y(',i1,')',
     &        8x,'z(',i1,')')
 1510 format( 5x,i3,3x,3f12.5)
 1520 format(3x,2x,'slat(',i1,') = ',f9.5,2x'flap(',i1,') = ',f9.5)

 1530 format(/3x,'reading sideview/v tail data, isect = ',i2)


c     *** Longitudinal derivatives first ***

      Ycg      = 0.0
     
c     *** Call 'CENTRAL' to get aerodynamic properties at AOA = 0 ***

      ALPHA    = 0.0
      q        = 0.0
      p        = 0.0
      ASYM     = 1.0
      ASYMDEFL = 1.0
      LAT      =   0
      IWING1   =  99
      IWING2   = -99

c     *** Constant ASYM is a indicator to VLM that it's to calculate
c         for asymmetric span load of equal but opposite gamma 
c         distribution ***

c     *** Constant IWING1 & IWING2 are the range of section numbers
c         that corresponds to the wing's sections (IWING1<=IWING2)***

      write(6,*)'GETTING NOMINAL AERODYNAMICS.'

      CALL DELZERO(ISECT,DELslat,DELflap)
      CALL CENTRAL(ISECT,DELslat,DELflap,ALPHA,CH,q,p,
     +             Xcg,Ycg,Zcg,V,ASYM,ASYMDEFL,LAT,VCref,GAMMAref,
     +             ITAIL1,ITAIL2,IWING1,IWING2,
     +             LIFT0,DRAG,PITCH0,ROLL,YAW0,LIFTAIL0,VCref,GAMMAref,
     +             bplon,slatlon,flaplon)

      Cdi            = 2.*DRAG/(QBAR*S)
      CL0            = 2.*LIFT0/(QBAR*S)
      CM0            = 2.*PITCH0/(QBAR*S*C)
      write(6,1405)  CL0,Cm0,Cdi
      write(12,1405) CL0,Cm0,Cdi
1405  format( 3x,'CL0                 = ',f12.5,
     1       /3x,'Cm0                 = ',f12.5,
     2       /3x,'Cdi(at alpha = 0)   = ',f12.5)
	
c     *** Skip if lateral derivatives only ***
      IF (ICASE .EQ. 2 ) GOTO 500

c     *** Get aerodynamic data at AOA = 5 ***

      ALPHA    = 5.0*PI/180.0
      q        = 0.0
      p        = 0.0
      ASYM     = 1.0
      ASYMDEFL = 1.0
      LAT      =   0

      write(6,*)'AT ALPHA = ', ALPHA,' rad.'

      CALL DELZERO(ISECT,DELslat,DELflap)
c
c     *** CH changed to reflect that of CG's height above ground ***

      CALL CENTRAL(ISECT,DELslat,DELflap,ALPHA,CH+Xcg*sin(alpha),q,
     +             p,Xcg,Ycg,Zcg,V,ASYM,ASYMDEFL,LAT,VCref,GAMMAref,
     +             ITAIL1,ITAIL2,IWING1,IWING2,
     +             LIFT,DRAG,PITCH,ROLL,YAW,LIFTTAIL,VC,GAMMA,
     +             bplon,slatlon,flaplon)
c
      CLalpha        = 2.*(LIFT-LIFT0)/(QBAR*S)/ALPHA
      Cmalpha        = 2.*(PITCH-PITCH0)/(QBAR*S*c)/alpha
      CmCL           =    (PITCH-PITCH0)/C/(LIFT-LIFT0)
      CDiCL2         =    2*DRAG*(QBAR*S)/(2*LIFT)**2
      CL             = 2.*LIFT/(QBAR*S)
      Cdi            = 2.*DRAG/(QBAR*S)

      write(6,1410)  CLalpha,Cmalpha, cmcl,CdiCL2,CL,Cdi 
      write(12,1410) CLalpha,Cmalpha, cmcl,CdiCL2,CL,Cdi
 1410 format( 3x,'CL-alpha            = ', f12.5,
     1       /3x,'Cm-alpha            = ', f12.5,
     2       /3x,'Cm/CL               = ', f12.5,
     3       /3x,'Cdi/CL^2            = ', f12.5,
     4       /3x,'CL (at alpha = 5)   = ', f12.5,
     5       /3x,'Cdi(at alpha = 5)   = ', f12.5)

      if(icase .eq. 4) goto 70
      
      if(itail1.eq.6.and.itail2.eq.6) goto 70
c     
c     *** Calc. lift-curve slope of tail in downwash ***

      CLATAIL1 = 2.0*(LIFTTAIL-LIFTAIL0)/(QBAR*S)/ALPHA

      write(6,*)'CLATAIL1 = ',CLATAIL1
c
c     *** Now run VLM without wings and determine the lift-curve
c         slope of the tail ***
c
      write(6,*) 'Performing wingless calculations:'

      IWING1   = IW1
      IWING2   = IW2
      ALPHA    = 0.0
      q        = 0.0
      p        = 0.0
      ASYM     = 1.0
      ASYMDEFL = 1.0
      LAT      =   0

      CALL CENTRAL(ISECT,DELslat,DELflap,ALPHA,CH,q,p,
     +             Xcg,Ycg,Zcg,V,ASYM,ASYMDEFL,LAT,VCref,GAMMAref,
     +             ITAIL1,ITAIL2,IWING1,IWING2,
     +             LIFT,DRAG,PITCH,ROLL,YAW,LIFTAIL0,VC,GAMMA,
     +             bplon,slatlon,flaplon)

      ALPHA = 1.0*PI/180.0
      CALL CENTRAL(ISECT,DELslat,DELflap,ALPHA,CH+Xcg*sin(alpha),q,
     +             p,Xcg,Ycg,Zcg,V,ASYM,ASYMDEFL,LAT,VCref,GAMMAref,
     +             ITAIL1,ITAIL2,IWING1,IWING2,
     +             LIFT,DRAG,PITCH,ROLL,YAW,LIFTTAIL,VC,GAMMA,
     +             bplon,slatlon,flaplon)
      write(6,*)'LIFT OF TAIL',LIFTTAIL

c     *** Determine lift-curve slope of tail w/o downwash ***
      CLATAIL2 = 2.0*(LIFTTAIL-LIFTAIL0)/(QBAR*S)/ALPHA
      write(6,*)'CLATAIL2 = ',CLATAIL2

c     *** Lift-curve slope of tial due to downwash ***
      LCSDW = CLATAIL2-CLATAIL1

      WRITE (6,1420) LCSDW
      WRITE(12,1420) LCSDW

 1420 format( 3x,'Lift-curve slope of',
     1       /3x,'tail due to downwash'
     2       /3x,'of wing             = ', f12.5)


 70   IWING1 = 99
      IWING2 = -99
c
c
c
c     *** Get pitch rate damping coefficients ***
c
      write(6,*)'PITCH RATE DAMPING COEFFICIENTS:'
      ALPHA    = 0.0
      q        = 0.1
      p        = 0.0
      ASYM     = 1.0
      ASYMDEFL = 1.0
      LAT      =   0
      CALL DELZERO(ISECT,DELslat,DELflap)
      CALL CENTRAL(ISECT,DELslat,DELflap,ALPHA,CH,q,p,
     +             Xcg,Ycg,Zcg,V,ASYM,ASYMDEFL,LAT,VCref,GAMMAref,
     +             ITAIL1,ITAIL2,IWING1,IWING2,
     +             LIFT,DRAG,PITCH,ROLL,YAW,LIFTTAIL,VC,GAMMA,
     +             bplon,slatlon,flaplon)
c
c
      CLq = 2*(LIFT-LIFT0)/(QBAR*S)/(q*(C/2/V))
      CMq = 2*(PITCH-PITCH0)/(QBAR*S*C)/(q*(C/2/V))

      WRITE (6,1430) CLq,CMq
      WRITE(12,1430) CLq,CMq
      if(lotus.eq.'y'.or.lotus.eq.'Y') then
	   open (14,file='lotus.out')
	   write(14,1431)Mach,Cl0,0.0,CM0,S,CLalpha,cref,Cmalpha,b,
     &        CmCl,xcgref,CDiCL2,zcg,CLq,ch,CMq
      endif
 
 1430 format( 3x,'CL-q                = ', f12.5,
     1       /3x,'CM-q                = ', f12.5)
 1431 format(1x,8(f12.5,3x,f12.5,/,1x))

      if(icase .eq. 4) GOTO 4000

c     *** Get longitudinal control effectiveness(symmetric deflection) ***

      if(icase.ne.4) then
      write(6,*)'PITCH CONTROL EFFECTIVENESS:'
      ALPHA    = 0.0
      q        = 0.0
      p        = 0.0
      ASYM     = 1.0
      ASYMDEFL = 1.0
      LAT      =   0
      DO 200, I = 1, ISECT
         DO 210, J = 1, 2
c
c           *** Catch those cases w/o control surfaces ***
            IF ((J .EQ. 1) .AND. (slatlon(I) .EQ. 0.0)) GOTO 210
            IF ((J .EQ. 2) .AND. (flaplon(I) .EQ. 0.0)) GOTO 210
c           *** Skip TE deflection if all moving surface ***
            IF ((J .EQ. 2) .AND. (slatlon(I) .EQ. 1.0)) GOTO 210
c
            CALL DELZERO(ISECT,DELslat,DELflap)
            IF (J .EQ. 1) DELslat(I) = 5.0*PI/180.0
            IF (J .EQ. 2) DELflap(I) = 5.0*PI/180.0
            CALL CENTRAL(ISECT,DELslat,DELflap,ALPHA,CH,q,p,
     +                   Xcg,Ycg,Zcg,V,ASYM,ASYMDEFL,LAT,VCref,GAMMAref,
     +                   ITAIL1,ITAIL2,IWING1,IWING2,
     +                   LIFT,DRAG,PITCH,ROLL,YAW,LIFTTAIL,VC,GAMMA,      
     +                   bplon,slatlon,flaplon)
            DCL(I,J) = 2*(LIFT-LIFT0)/(QBAR*S)/(5.0*PI/180.0)
            DCM(I,J) = 2*(PITCH-PITCH0)/(QBAR*S*C)/(5.0*PI/180.0)   
            WRITE(12,2000) I,J, DCL(I,J)
            WRITE(12,2001) I,J, DCm(I,J)
            WRITE (6,2000) I,J, DCL(I,J)
            WRITE (6,2001) I,J, DCm(I,J)  
 210     CONTINUE
 200  CONTINUE
      endif

 2000 FORMAT (3x,'CL-delta [',I1,' ',I1,']      = ',f12.5)
 2001 FORMAT (3x,'Cm-delta [',I1,' ',I1,']      = ',f12.5)

c
c
c     *** Stop if Longitudinal derivatives only ***
      IF (ICASE .EQ. 1) GOTO 4000
c
c
c     *** LATERAL DERIVATIVES FROM HERE DOWN ***
c
c
c     *** Note that when non-dimensionalizinging yaw-rate & yaw moment
c         must account for compressibility distortion by dividing
c         the b (span) by BETA for each yaw-rate and yaw moment.
c         This does not apply to roll (moment and rate) nor
c         sideforce. ***

c     *** Get roll rate damping coefficient (w/o Vertical Tail) ***
 500  write(6,*)'ROLL RATE DAMPING COEFFICIENT:'
      ALPHA    =  0.0
      q        =  0.0
      p        =  .10
      ASYM     = -1.0
      ASYMDEFL =  1.0
      LAT      =    0
      CALL DELZERO(ISECT,DELslat,DELflap)
      CALL CENTRAL(ISECT,DELslat,DELflap,ALPHA,CH,q,p,
     +             Xcg,Ycg,Zcg,V,ASYM,ASYMDEFL,LAT,VCref,GAMMAref,
     +             ITAIL1,ITAIL2,IWING1,IWING2,
     +             LIFT,DRAG,PITCH,ROLL,YAW,LIFTTAIL,VC,GAMMA,
     +             bplon,slatlon,flaplon)

c     *** Multiplied by 2 for the other half of configuration ***
      CROLLp = 2.0*ROLL/(p*(B/2/V))/(QBAR*S*B)
      CYAWp  = 2.0*YAW/(p*(B/BETA/2/V))/(QBAR*S*B)


c     *** Hold off printing until counting the V. tail contribution ***

c     *** Determine roll control power (of each control surface
c         when deflected antisymmetrically about the x-axis ***

      if(icase.ne.4) then
      write(6,*)'Determining roll control EFFECTIVENESS:'
      ALPHA = 0.0
      p        =  0.0
      q        =  0.0
      ASYM     = -1.0
      ASYMDEFL = -1.0
      LAT      =    0
c
c
      DO 220, I = 1, ISECT
         DO 230, J = 1, 2
c
c           *** Catch those cases w/o control surfaces ***
            IF ((J .EQ. 1) .AND. (slatlon(I) .EQ. 0.0)) GOTO 230
            IF ((J .EQ. 2) .AND. (flaplon(I) .EQ. 0.0)) GOTO 230
c           *** Skip TE deflection if all moving surface ***
            IF ((J .EQ. 2) .AND. (slatlon(I) .EQ. 1.0)) GOTO 230
c
            CALL DELZERO(ISECT,DELslat,DELflap)
            IF (J .EQ. 1) DELslat(I) = 5.0*PI/180.0
            IF (J .EQ. 2) DELflap(I) = 5.0*PI/180.0
C
            CALL CENTRAL(ISECT,DELslat,DELflap,ALPHA,CH,q,p,
     +             Xcg,Ycg,Zcg,V,ASYM,ASYMDEFL,LAT,VCref,GAMMAref,
     +             ITAIL1,ITAIL2,IWING1,IWING2,
     +             LIFT,DRAG,PITCH,ROLL,YAW,LIFTTAIL,VC,GAMMA,
     +             bplon,slatlon,flaplon)
            write(6,*) 'yaw = ', yaw
            DCn(I,J) = 2.0*YAW/(5.0*PI/180.0)/(QBAR*S*B/BETA)      
            DCroll(I,J) = -2.0*ROLL/(5.0*PI/180.0)/(QBAR*S*B)
            WRITE(12,2002) I,J, DCroll(I,J)
            WRITE(12,2003) I,J, DCn(I,J)
            WRITE (6,2002) I,J, DCroll(I,J)
            WRITE (6,2003) I,J, DCn(I,J)
 230     CONTINUE
 220  CONTINUE
      endif

 2002 FORMAT (3x,'Cl-delta [',I1,' ',I1,']      = ',f12.5)
 2003 FORMAT (3x,'Cn-delta [',I1,' ',I1,']      = ',f12.5)


c     *** Estimating directional derivatives ***

      write(6,*)'DETERMINING C-y-beta & C-n-beta:'

c     *** Exchange y- & z-coordinates of CG for lat. stab. derivatives ***

      Ycg = Zcg
      Zcg = 0.0

c
c     *** Estimating sideforce, yaw & roll moment due to side slip ***

      ALPHA = 5.0*PI/180.0
c
c     *** Get out of Ground effect for Lat/directional derivatives ***

      CH       = 9E9
      q        = 0.0
      p        = 0.0
      ASYM     = 1.0
      ASYMDEFL = 1.0
      LAT      =   1
      CALL DELZERO(JSECT,DELslat,DELflap)
      CALL CENTRAL(JSECT,DELslat,DELflap,ALPHA,CH,q,p,
     +             Xcg,Ycg,Zcg,V,ASYM,ASYMDEFL,LAT,VCref,GAMMAref,
     +             ITAIL1,ITAIL2,IWING1,IWING2,
     +             LIFT,DRAG,PITCH,ROLL,YAW,LIFTTAIL,VC,GAMMA,
     +             bplat,slatlat,flaplat)

c
      Cybeta = -LIFT/(QBAR*S)/ALPHA
      Cnbeta = -PITCH/(QBAR*S*B/BETA)/ALPHA
      Clbeta = -ROLL/(QBAR*S*B)/ALPHA

      WRITE (6,1435) Cybeta,Cnbeta,Clbeta
      WRITE(12,1435) Cybeta,Cnbeta,Clbeta
 
 1435 format( 3x,'Cy-beta             = ', f12.5,
     1       /3x,'Cn-beta             = ', f12.5,
     2       /3x,'Cl-beta             = ', f12.5)
    

c     *** Estimating yaw rate coefficients ***

      write(6,*)'YAW RATE DAMPING COEFFICIENTS:'
      ALPHA    = 0.0
      CH       = 9E9
      q        = .01
      p        = 0.0
      ASYM     = 1.0
      ASYMDEFL = 1.0
      LAT      =   1
      CALL DELZERO(JSECT,DELslat,DELflap)
      CALL CENTRAL(JSECT,DELslat,DELflap,ALPHA,CH,q,p,
     +             Xcg,Ycg,Zcg,V,ASYM,ASYMDEFL,LAT,VCref,GAMMAref,
     +             ITAIL1,ITAIL2,IWING1,IWING2,
     +             LIFT,DRAG,PITCH,ROLL,YAW,LIFTTAIL,VC,GAMMA,
     +             bplat,slatlat,flaplat)

      Cyr = LIFT/(QBAR*S)/(q*(B/BETA/2/V))
      Cnr = PITCH/(QBAR*S*B/BETA)/(q*(B/BETA/2/V))
      Clr = ROLL/(QBAR*S*B)/(q*(B/BETA/2/V))

      WRITE (6,1440) Cyr,Cnr,Clr
      WRITE(12,1440) Cyr,Cnr,Clr

 1440 format( 3x,'Cy-r                = ', f12.5,
     1       /3x,'Cn-r                = ', f12.5,
     2       /3x,'Cl-r                = ', f12.5)


c     *** Roll rate coefficient due to vertical tail ***

      write(6,*)'ROLL RATE DAMPING COEFFICIENT (due to V. tail):'
      ALPHA    = 0.0
      CH       = 9E9
      q        = 0.0
      p        = 1.0
      ASYM     = 1.0
      ASYMDEFL = 1.0
      LAT      =   1
      CALL DELZERO(JSECT,DELslat,DELflap)
      CALL CENTRAL(JSECT,DELslat,DELflap,ALPHA,CH,q,p,
     +             Xcg,Ycg,Zcg,V,ASYM,ASYMDEFL,LAT,VCref,GAMMAref,
     +             ITAIL1,ITAIL2,IWING1,IWING2,
     +             LIFT,DRAG,PITCH,ROLL,YAW,LIFTTAIL,VC,GAMMA,
     +             bplat,slatlat,flaplat)
c
c     *** Add to wing and H tail's contribution ***
      CROLLp = CROLLp+ROLL/(p*(B/2/V))/(QBAR*S*B)
      CYAWp = CYAWp+PITCH/(p*(B/2/v))/(QBAR*S*B/BETA)

      WRITE (6,1450) CROLLp,CYAWp
      WRITE(12,1450) CROLLp,CYAWp

 1450 format( 3x,'Cl-p                = ', f12.5,
     1       /3x,'Cn-p                = ', f12.5)


c     *** Yaw effector coefficients ***

      if(icase.ne.4) then
      write(6,*)'RUDDER EFFECTIVENESS:'
      ALPHA    = 0.0
      CH       = 9E9
      q        = 0.0
      p        = 0.0
      ASYM     = 1.0
      ASYMDEFL = 1.0
      LAT      =   1
      DO 260, I = 1, JSECT
         DO 270, J = 1, 2

c           *** Catch those cases w/o control surfaces ***
            IF ((J .EQ. 1) .AND. (slatlat(I) .EQ. 0.0)) GOTO 270
            IF ((J .EQ. 2) .AND. (flaplat(I) .EQ. 0.0)) GOTO 270
c
            CALL DELZERO(JSECT,DELslat,DELflap)
            IF (J .EQ. 1) DELslat(I) = -5.0*PI/180.0
            IF (J .EQ. 2) DELflap(I) = -5.0*PI/180.0
            CALL CENTRAL(JSECT,DELslat,DELflap,ALPHA,CH,q,p,
     +                   Xcg,Ycg,Zcg,V,ASYM,ASYMDEFL,LAT,VCref,GAMMAref,
     +                   ITAIL1,ITAIL2,IWING1,IWING2,
     +                   LIFT,DRAG,PITCH,ROLL,YAW,LIFTTAIL,VC,GAMMA,      
     +                   bplat,slatlat,flaplat)
            DCy(I,J) = -LIFT/(QBAR*S)/(5.0*PI/180.0)
            DCn(I,J) = -PITCH/(QBAR*S*B/BETA)/(5.0*PI/180.0)
            DCl(I,J) = -ROLL/(QBAR*S*B)/(5.0*PI/180.0)   
            WRITE(12,2004) I,J, DCy(I,J)
            WRITE(12,2005) I,J, DCl(I,J)  
            WRITE(12,2006) I,J, DCN(I,J)
            WRITE (6,2004) I,J, DCy(I,J)
            WRITE (6,2005) I,J, DCl(I,J)  
            WRITE (6,2006) I,J, DCN(I,J)
 270     CONTINUE
 260  CONTINUE
      endif

 2004 FORMAT (3x,'Cy-delta [',I1,' ',I1,']      = ',f12.5)
 2005 FORMAT (3x,'Cl-delta [',I1,' ',I1,']      = ',f12.5)
 2006 FORMAT (3x,'Cn-delta [',I1,' ',I1,']      = ',f12.5)

 4000	WRITE(6,3100)
	WRITE(12,3100)
 3100	FORMAT(/'  Run completed successfully.'/)

	CLOSE(6)
	CLOSE(12)

	PAUSE '  Press enter to close program.'

      STOP
      END

c     *******************************************************************
      SUBROUTINE DELZERO(ISECT,DELslat,DELflap)
      REAL DELslat(5), DELflap(5)
      INTEGER ISECT
c
c     ===================================================================
c     This subroutine sets all surface deflection to zero.
c     ===================================================================
c
      DO 100, I = 1, ISECT
         DELslat(I) = 0.0
         DELflap(I) = 0.0
 100  CONTINUE
      RETURN
      END

      SUBROUTINE CENTRAL(ISECT,DELslat,DELflap,ALPHA,CH,q,p,
     +                   Xcg,Ycg,Zcg,V,ASYM,ASYMDEFL,LAT,VCref,GAMMAref,
     +                   ITAIL1,ITAIL2,IWING1,IWING2,
     +                   LIFT,DRAG,PITCH,ROLL,YAW,LIFTTAIL,VC,GAMMA,
     +                   bpin,cslat,cflap)

      REAL VC(5,10,6,3), BP(4,3), OUT(10,6,3), Cslat(5), Cflap(5)
      REAL Hslat(5,2,3), Hflap(5,2,3), DHslat(2,3), DHflap(2,3)
      REAL DELslat(5), DELflap(5), CP(5,8,5,3),VNOR(5,8,5,4)
      REAL VCref(5,10,6,3),GAMMA(5,8,10),GAMMAref(5,8,10)
      REAL LIFT,PITCH,ROLL,LIFTTAIL,bpin(5,4,3)

      INTEGER  ISECT, IP(5), JP(5), IPslat(5), IPflap(5), LAT


      OPEN (UNIT = 15, FILE = 'RESULTS')
      REWIND (UNIT = 15)

c     ===================================================================
c     This program controls the input of VLM and its grid geometry.
c     It repeatedly calls 'GEOMETRY' to determine the grid points of 
c     each section of the configuration.  Calling 'DEFLECTION' rotates
c     the points on the flap/slats about the hinge lines of the sections.
c     'NORMAL' determines the normal vectors at each control point.
c     After the geometry is set up, VLM is called.

c     Prepared by Jacob Kay,   Sept, 1992.
c     =====================================================================
c
      DO 90, I = 1, ISECT
         IP(I) = 8
         JP(I) = 5
 90   CONTINUE

c     *** For each section ***

      DO 100, K = 1, ISECT


 1000    FORMAT (3F10.5)
         do 30 i = 1,4
	      do 20 j = 1,3
               bp(i,j)=bpin(k,i,j)
   20       continue
   30    continue
c        *** Calculate vortex ring corner point coordinates ***
         CALL GEOMETRY(BP,OUT,ALPHA,IP(K),JP(K),Cslat(K),Cflap(K),
     +                 DHslat,DHflap,IPslat(K),IPflap(K))
c
c        *** Fill vortex ring corner pt coord. matrix ***
         DO 120, I = 1, IP(K)+2
            DO 130, J = 1, JP(K)+1
               VC(K,I,J,1) = OUT(I,J,1)
               VC(K,I,J,2) = OUT(I,J,2)
               VC(K,I,J,3) = OUT(I,J,3)
 130        CONTINUE
 120     CONTINUE
c
c        *** Fill hinge line end point coord. matrix ***
         DO 140, I = 1, 2
            DO 150, J = 1, 3
               Hslat(K,I,J) = DHslat(I,J)
               Hflap(K,I,J) = DHflap(I,J)
 150        CONTINUE
 140     CONTINUE
c
c
 100  CONTINUE
      write(6,*)'Finished defining vortex ring corners.'
c
c
c     *** Reconfigure geometry when deflecting ***
c
      write(15,1011)
c     *** For each section ***
      DO 200, K = 1, ISECT
c
c        *** Slat first ***
c
c        *** Skip if no deflection ***
         IF (DELslat(K) .EQ. 0.0) GOTO 210
c        
c        *** For all-moving surface ***
         DO 230, J = 1, JP(K)+1
            DO 240, I = 1,IPslat(K)+1
               CALL DEFLECT(Hslat(K,1,1),Hslat(K,1,2),Hslat(K,1,3),
     +                      Hslat(K,2,1),Hslat(K,2,2),Hslat(K,2,3),
     +                      VC(K,I,J,1),VC(K,I,J,2),VC(K,I,J,3),
     +                      DELslat(K))
 240        CONTINUE
 230     CONTINUE
      write(6,*)'Slat deflected.'
c
c        *** Flap ***
 210     IF (DELflap(K) .EQ. 0.0) GOTO 250
         DO 260, J = 1, JP(K)+1
            DO 270, I = IP(K)-IPflap(K)+2,IP(K)+1
               CALL DEFLECT(Hflap(K,1,1),Hflap(K,1,2),Hflap(K,1,3),
     +                      Hflap(K,2,1),Hflap(K,2,2),Hflap(K,2,3),
     +                      VC(K,I,J,1),VC(K,I,J,2),VC(K,I,J,3),
     +                      DELflap(K))
 270        CONTINUE
 260     CONTINUE
      write(6,*)'Flap deflected.'
c
c
c        *** Print 'em out to check ***
 250     WRITE(15,1001)k
         DO 280, J = 1, JP(K)+1
            DO 290, I = 1, IP(K)+1
               WRITE(15,1000) VC(K,I,J,1),VC(K,I,J,2),VC(K,I,J,3)
 290        CONTINUE
 280     CONTINUE
c
c
 200  CONTINUE
 1001 FORMAT ('zone t="Section ',i1,'",i=9,j=6,f=point')
 1011 format ('Title="VLM Planform Output"',/,'variables=x,y,z')

c       
c
c     *** Determine control points and their normal vectors ***
      CALL CONPT(VC,CP,VNOR,IP,JP,ISECT)
      write(6,*)'Defining control points and normal vectors.'
c
c     *** PRINT 'EM OUT TO CHECK ***
      WRITE(15,1002)
      DO 300, K = 1, ISECT
         DO 310, J = 1, JP(K)
            DO 320, I = 1, IP(K)
               WRITE(15,1000) CP(K,I,J,1), CP(K,I,J,2), CP(K,I,J,3)
 320        CONTINUE
 310     CONTINUE
 300  CONTINUE
 1002 FORMAT ('zone t="CONTROL POINTS"')
c
c     *** CHECK NORMAL VECTORS ***
      WRITE(15,1003)
      DO 400, K = 1, ISECT
         DO 410, J = 1, JP(K)
            DO 420, I = 1, IP(K)
               WRITE(15,1000) VNOR(K,I,J,1),VNOR(K,I,J,2), VNOR(K,I,J,3)
 420        CONTINUE
 410     CONTINUE
 400  CONTINUE
 1003 FORMAT('zone t="NORMAL VECTORS"')
 
c
c 
c     *** Perform VLM calculation ***
      write(6,*) 'Entering VLM.'
      CALL VLM(VC,CP,VNOR,IP,JP,ISECT,ALPHA,CH,V,Xcg,Ycg,Zcg,q,p,ASYM,
     +         VCref,GAMMAref,ASYMDEFL,DELslat,DELflap,Hslat,Hflap,
     +         IPslat,IPflap,LAT,ITAIL1,ITAIL2,IWING1,IWING2,
     +         LIFT,DRAG,PITCH,ROLL,YAW,LIFTTAIL,GAMMA)
C
C
      END


c     *******************************************************************
      SUBROUTINE GEOMETRY(BP,VC,ALPHA,IP,JP,Cslat,Cflap,
     +                    Hslat,Hflap,IPslat,IPflap)

      REAL BP(4,3), Cslat, Cflap, ALPHA, VC(10,6,3), LAMBDA
      REAL Hslat(2,3), Hflap(2,3)

      INTEGER IP, JP, IPslat, IPflap, IPunflap

c     ===================================================================
c     This subroutine will use the corner points of sections of planforms
c     and determine the corner points of vortex rings.
c
c     Prepared by Jacob Kay, VPI&SU  Sept, 1992.
c     ===================================================================

c     *** Begin Loop ***


c     *** Some constants for the section ***

      B         = BP(2,2)-BP(1,2)
      DY        = B/JP
      Xfract   = 1.0/IP
      IPslat   = NINT(Cslat/Xfract)
      IPflap    = NINT(Cflap/Xfract)
      Cr        = BP(4,1)-BP(1,1)
      Ct        = BP(3,1)-BP(2,1)
      LAMBDA    = Ct/Cr
      THETAr    = -ASIN((BP(4,3)-BP(1,3))/(BP(4,1)-BP(1,1)))
      THETAt    = -ASIN((BP(3,3)-BP(2,3))/(BP(3,1)-BP(2,1)))
c
c     *** At each spanwise station: ***
      DO 110, K = 1, JP+1
         Y      = DY*(K-1)
         CHORD  = Cr+(Ct-Cr)*Y/B
         THETA  = THETAr+ATAN((LAMBDA*TAN(THETAt-THETAr)*Y/B)/
     +            (1-(1-LAMBDA)*Y/B))
         XLE    = BP(1,1)+(BP(2,1)-BP(1,1))*Y/B
         ZLE    = BP(1,3)+(BP(2,3)-BP(1,3))*Y/B
         XDIS   = .25*(CHORD/IP)
c
c        *** Define vortex rings' corners on slat ***
c    
c        *** Skip over if no slat ***
         IF (IPslat .EQ. 0) THEN
            DX = 0.0
         else 
            DX           = Cslat*CHORD/IPslat
         endif
 50      DO 120, J    = 1, IPslat+1
            XINC      = DX*(J-1)
            X         = XLE+XINC+XDIS
            VC(J,K,1) = X*COS(ALPHA)-(ZLE-(X-XLE)*
     +                  TAN(THETA))*SIN(ALPHA)
            VC(J,K,2) = Y+BP(1,2)
            VC(J,K,3) = -X*SIN(ALPHA)+(ZLE-(X-XLE)*
     +                   TAN(THETA))*COS(ALPHA)
 120     CONTINUE
c        *** Skip to find hinge line end pts if all-moving surface ***
         IF (Cslat.EQ.1.0) GO TO 110
c
c
c        *** Define vortex rings'corners on flap ***
         XHINGE = BP(4,1)+(BP(3,1)-BP(4,1))*Y/B-CfLAP*CHORD
c        *** If no flap ***
         IF (IPflap .EQ. 0) THEN
            DX = 0.0
         else 
            DX = Cflap*CHORD/IPflap
         endif
 60      DO 130, J = IP+1-IPflap, IP+1
            XINC = DX*(J-IP+IPflap-1)
            X = XHINGE+XDIS+XINC
            VC(J,K,1) = X*COS(ALPHA)-(ZLE-(X-XLE)*
     +                  TAN(THETA))*SIN(ALPHA)
            VC(J,K,2) = Y+BP(1,2)
            VC(J,K,3) = -X*SIN(ALPHA)+(ZLE-(X-
     +                  XLE)*TAN(THETA))*COS(ALPHA)
 130     CONTINUE
c
c        *** Define vortex rings' corners on unflapped section ***
         IPunflap = IP-IPslat-IPflap
         DX = (VC(IP+1-IPflap,K,1)-VC(IPslat+1,K,1))/IPunflap
         DZ = (VC(IP+1-IPflap,K,3)-VC(IPslat+1,K,3))/IPunflap
         DO 140, J = IPslat+2, IP-IPflap
            X = VC(IPslat+1,K,1)+(J-IPslat-1)*DX
            Z = VC(IPslat+1,K,3)+(J-IPslat-1)*DZ
            VC(J,K,1) = X
            VC(J,K,2) = Y+BP(1,2)
            VC(J,K,3) = Z
 140     CONTINUE
c
c
c        *** Define downstream wake points ***
         VC(IP+2,K,1) = VC(IP+1,K,1)+100*B
         VC(IP+2,K,2) = VC(IP+1,K,2)
         VC(IP+2,K,3) = VC(IP+1,K,3)
c
 110  Continue
c
c     *** Hinge lines ***
      IF (Cslat .EQ. 1.0) GOTO 300
c
c     *** End points of flap and slat hinge lines ***
      DO 200, I = 1, 3
         Hslat(1,I) = VC(IPslat+1,1,I)
         Hslat(2,I) = VC(IPslat+1,JP+1,I)
         Hflap(1,I) = VC(IP-IPflap+1,1,I)
         Hflap(2,I) = VC(IP-IPflap+1,JP+1,I)
 200  CONTINUE
c
c
c
      RETURN
c
c
c     *** Find hinge line end points if all-moving surface ***
c
c     *** One pt going thru c/4 of mean chord ***
 300  Ymac       = B*(1.0+2.0*LAMBDA)/(3.0*(1.0+LAMBDA))+BP(1,2)
      Cmac       = Cr+(Ct-Cr)*Ymac/B
      XLE        = BP(1,1)+(BP(2,1)-BP(1,1))*Ymac/B
      ZLE        = BP(1,3)+(BP(2,3)-BP(1,3))*Ymac/B
      X          = XLE+.25*Cmac
      THETA      = THETAr+ATAN((LAMBDA*TAN(THETAt-THETAr)*Y/B)/
     +             (1-(1-LAMBDA)*Y/B))
      Hslat(2,1) = X*COS(ALPHA)-(ZLE-(X-XLE)*TAN(THETA))*SIN(ALPHA)
      Hslat(2,2) = Ymac
      Hslat(2,3) = -X*SIN(ALPHA)+(ZLE-(X-XLE)*TAN(THETA))*COS(ALPHA) 
c
c     *** The other pt forms a line perpendicular to A/C centerline ***
      XLE          = BP(1,1)
      ZLE          = BP(1,3)
      Hslat(1,1)   = Hslat(2,1)
      Hslat(1,2)   = BP(1,2)
      THETA        = THETAr
      X            = (ZLE*SIN(ALPHA)+XLE*TAN(THETA)*SIN(ALPHA))/
     +               (COS(ALPHA)+SIN(ALPHA)*TAN(THETA))
      Hslat(1,3)   = -X*SIN(ALPHA)+(ZLE-(X-XLE)*TAN(THETA))*COS(ALPHA)

      RETURN                                             

      END



c     ******************************************************************
      SUBROUTINE DEFLECT(A1,A2,A3,B1,B2,B3,P1,P2,P3,DEL)
      REAL C(3), MAG, T(3,3), TINV(3,3)
c
c     ===================================================================
c     This subroutine calculates the new coordinates of a point
c     after been rotated certain degrees about an arbitrary axis of
c     rotation (described by two points).
c
c     Prepared by Jacob Kay, Sept, 1992.
c     ===================================================================
c
c
C     PI = 4.0*ATAN(1.0)
c
c     *** Find point C on a line-segment AB ***
      MAG = -((P1-A1)*(A1-B1)+(P2-A2)*(A2-B2)+
     +        (P3-A3)*(A3-B3))/((A1-B1)**2+(A2-B2)**2+
     +        (A3-B3)**2)
      C(1) = A1+MAG*(B1-A1)
      C(2) = A2+MAG*(B2-A2)
      C(3) = A3+MAG*(B3-A3)
c
c     *** Vector from A to B ***
      X0 = B1-A1
      Y0 = B2-A2
      Z0 = B3-A3
c     *** Vector from C to P ***
      X1 = P1-C(1)
      Y1 = P2-C(2)
      Z1 = P3-C(3)
c
c     *** Calculate transformation to rotate the inertial x-axis to 
c         be parallel to the vector from A to B ***
      PSI   =  ATAN2(Y0,X0)
      THETA = -ASIN(Z0/SQRT(X0**2+Y0**2+Z0**2))
c
c     *** Get transformation matrix (including deflection) ***
      CALL TRANSFORM(PSI,THETA,-DEL,T)
c
c     *** Vector from C to P' (P after deflection) expressed in axis 
c         system where x-axis is parallel to the vector from A to B ***
      X2P = T(1,1)*X1+T(1,2)*Y1+T(1,3)*Z1
      Y2P = T(2,1)*X1+T(2,2)*Y1+T(2,3)*Z1
      Z2P = T(3,1)*X1+T(3,2)*Y1+T(3,3)*Z1
c
c     *** Getting back to the inertial axis system ***

      CALL TRANSFORM(PSI,THETA,0.0,T)
c
c     *** Transpose of T ***
      DO 100, I       = 1,3
         DO 110, J    = 1,3
            TINV(I,J) = T(J,I)
 110     CONTINUE
 100  CONTINUE
c
c     *** Vector from C to P' expressed in inertial axis system ***
      X2 = TINV(1,1)*X2P+TINV(1,2)*Y2P+TINV(1,3)*Z2P
      Y2 = TINV(2,1)*X2P+TINV(2,2)*Y2P+TINV(2,3)*Z2P
      Z2 = TINV(3,1)*X2P+TINV(3,2)*Y2P+TINV(3,3)*Z2P
c     *** Position of P' (P after deflection) in inertial axis system ***
      P1 = C(1)+X2
      P2 = C(2)+Y2
      P3 = C(3)+Z2
c
c
c
      RETURN
      END

c     ******************************************************************
      SUBROUTINE TRANSFORM(PSI, THETA, PHI, A)
      REAL A(3,3), PSI, THETA, PHI
c
c     ==================================================================
c     This subroutine determines the transformation matrix for a axis
c     system rotation psi about the z-axis, then theta about the y'-axis
c     and finally phi about the x"-axis.
c     ==================================================================
c
c     *** CALC. MATRIX ELEMENT ***
      A(1,1) = COS(THETA)*COS(PSI)
      A(1,2) = COS(THETA)*SIN(PSI)
      A(1,3) = -SIN(THETA)
      A(2,1) = SIN(PHI)*SIN(THETA)*COS(PSI)-COS(PHI)*SIN(PSI)
      A(2,2) = SIN(PHI)*SIN(THETA)*SIN(PSI)+COS(PHI)*COS(PSI)
      A(2,3) = SIN(PHI)*COS(THETA)
      A(3,1) = COS(PHI)*SIN(THETA)*COS(PSI)+SIN(PHI)*SIN(PSI)
      A(3,2) = COS(PHI)*SIN(THETA)*SIN(PSI)-SIN(PHI)*COS(PSI)
      A(3,3) = COS(PHI)*COS(THETA)
      RETURN	
      END


c     *******************************************************************
      SUBROUTINE CONPT(QF,QC,DS,IP,JP,ISECT)
      REAL QF(5,10,6,3),QC(5,8,5,3),DS(5,8,5,4)
      INTEGER IP(5), JP(5)
c
c
c     ==================================================================
c     This subroutine determines the locations and the normal vectors
c     of control points of each section.
c     ==================================================================
c     *** Control Points ***
      DO 100, K = 1, ISECT
         DO 120, J = 1, JP(K)
            DO 130, I = 1, IP(K)
               QC(K,I,J,1) = (QF(K,I,J,1)+QF(K,I,J+1,1)+
     +                        QF(K,I+1,J+1,1)+QF(K,I+1,J,1))/4.0
               QC(K,I,J,2) = (QF(K,I,J,2)+QF(K,I,J+1,2)+
     +                        QF(K,I+1,J+1,2)+QF(K,I+1,J,2))/4.0
               QC(K,I,J,3) = (QF(K,I,J,3)+QF(K,I,J+1,3)+
     +                        QF(K,I+1,J+1,3)+QF(K,I+1,J,3))/4.0
c
c           *** Calulate panels' normal vectors ***
            CALL NORMAL (QF(K,I,J,1),QF(K,I,J,2),QF(K,I,J,3),
     +	                  QF(K,I+1,J,1),QF(K,I+1,J,2),QF(K,I+1,J,3),
     +                   QF(K,I,J+1,1),QF(K,I,J+1,2),QF(K,I,J+1,3),
     +                   QF(K,I+1,J+1,1),QF(K,I+1,J+1,2),
     +                   QF(K,I+1,J+1,3),DS(K,I,J,1),
     +                   DS(K,I,J,2),DS(K,I,J,3),DS(K,I,J,4))
 130        CONTINUE
 120     CONTINUE
 100  CONTINUE
c
c
      RETURN
      END


c     ********************************************************************
      SUBROUTINE NORMAL(x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,c1,c2,c3,s)
c
c     ====================================================================
c     This subroutine calculates the average normal vector of a 4-corner
c     panel and its area.
c     ====================================================================
c
      xa   = x2-x3
      ya   = y2-y3
      za   = z2-z3
      xb   = x4-x1
      yb   = y4-y1
      zb   = z4-z1
c     *** Normal vector components ***
      x    = ya*zb-za*yb
      y    = xb*za-xa*zb
      z    = xa*yb-ya*xb
      amag = SQRT(x**2+y**2+z**2)
c     *** Normalizing the normal vector ***
      c1   = x/amag
      c2   = y/amag
      c3   = z/amag
c     *** Normal area ***
      e1   = x3-x1
      e2   = y3-y1
      e3   = z3-z1
      f1   = x2-x1
      f2   = y2-y1
      f3   = z2-z1
c   
      s11  = f2*zb-f3*yb
      s12  = xb*f3-f1*zb
      s13  = f1*yb-f2*xb
      s21  = yb*e3-zb*e2
      s22  = e1*zb-xb*e3
      s23  = xb*e2-yb*e1
      s    = .5*(SQRT(s11**2+s12**2+s13**2)+SQRT(s21**2+s22**2+s23**2))
c
c
      RETURN
      END


c     ******************************************************************
      SUBROUTINE VLM(QF,QC,DS,IP,JP,ISECT,ALPHA,CH,VF,Xcg,Ycg,Zcg,
     +               q,p,ASYM,
     +               QFref,GAMMAref,ASYMDEFL,DELslat,DELflap,Hslat,
     +               Hflap,IPslat,IPflap,LAT,ITAIL1,ITAIL2,
     +               IWING1,IWING2,
     +               LIFT,DRAG,PITCH,roll,YAW,LIFTTAIL,GAMMA)
      REAL LIFT,PITCH,ROLL, LIFTTAIL,xmact,zmact

      common /consts/ PI
      DIMENSION QF(5,10,6,3), QC(5,8,5,3), DS(5,8,5,4)
      DIMENSION GAMMA(5,8,10),GAMMAref(5,8,10),QFref(5,10,6,3)
      DIMENSION A(200,200), GAMMA1(200), DW(200), SYM(200)
      DIMENSION GAMMAnew(5,8,10),DELslat(5),DELflap(5)
      DIMENSION Hslat(5,2,3),Hflap(5,2,3),IPslat(5),IPflap(5)
      DIMENSION A1(5,9,5), JP(5), IP(5)
c    
c     ==================================================================
c     This program uses a lattice of closed vortex rings to estimate
c     aerodynamic properties of a planform.  
c
c     Prepared by Jacob Kay		Sept. 1992.
c     ==================================================================
c
c
      RHO = 1.0
c
      pi  = 4.0*ATAN(1.0)
c
c     *** Set gamma to 1 for influence coefficient matrix calulation. ***
      DO 102, K = 1, 5
         DO 100, I = 1, 8
            DO 101, J = 1, 10
               GAMMA(K,I,J) = 1.0
 101        CONTINUE
 100     CONTINUE
 102  CONTINUE
      sn1 = SIN(alpha)
      cs1 = COS(alpha)
c
c     *** For antisymmetric surface deflections, the reference gamma
c         distribution has to be reconfigured to account for the change
c         in geometry. ***
      IF (ASYM .EQ. -1.0) THEN
      DO 50, M = 1, ISECT
         DO 55, I = 1, 1
            DO 55, J = 1, JP(M)
               GAMMAnew(M,I,J) = GAMMAref(M,I,J)*(QFref(M,I,J+1,2)-
     +                           QFref(M,I,J,2))/
     +                           (QF(M,I,J+1,2)-QF(M,I,J,2))
 55      CONTINUE
         DO 50, I = 2, IP(M)
            DO 50, J = 1, JP(M)
               GAMMAnew(M,I,J) = (GAMMAref(M,I,J)-GAMMAref(M,I-1,J))*
     +                           (QFref(M,I,J+1,2)-QFref(M,I,J,2))/
     +                           (QF(M,I,J+1,2)-QF(M,I,J,2))+
     +                           GAMMAnew(M,I-1,J)
 50   CONTINUE
      ENDIF 
c
c
c     *** Influence Coefficients ***
c     *** Sum up the induced velocity due to all vortex rings at each
c         control point. ***
c
      K = 0
c     *** For each control pt. ***
      DO 105, M = 1, ISECT
      IF ((M.GE.IWING1).AND.(M.LE.IWING2)) GOTO 105
      write(6,*)'CALCULATING INFLUENCE COEFFICIENTS FOR SECTION',M
      DO 110, I = 1, IP(M)
         DO 120, J = 1, JP(M)
            K = K+1

 2001 format(2x,'m=',i2,2x,'i=',i2,2x,'ip(m)=',i2,2x,'j=',i2,2x,
     1          'jp(m)=',i2,2x,'k=',i3)


c           *** The original half of wing ***
            CASE = 1.0
            CALL WING(QC(M,I,J,1),QC(M,I,J,2),QC(M,I,J,3),GAMMA,u,v,w,
     +                1.0,M,I,J,DS,ip,jp,CASE,A1,QF,ISECT,
     +                IWING1,IWING2)
c
            L = 0
            DO 125, M1 = 1, ISECT
               IF ((M1.GE.IWING1).AND.(M1.LE.IWING2)) GOTO 125
               DO 130, I1 = 1, IP(M1)
                  DO 140, J1 = 1, JP(M1)
                     L = L+1
                     A(K,L) = A1(M1,I1,J1)
 140              CONTINUE
 130           CONTINUE
 125        CONTINUE
c
c           *** The other half of wing ***
            CASE = 2.0
c
c           *** Skip if lateral case ***
            IF (LAT .EQ. 1) GOTO 170
c
c           *** Reverse flap deflection for antisymmetric deflection ***
            IF (ASYMDEFL .EQ. -1.0) CALL REVERSE(DELslat,DELflap,
     +                              QF,Hslat,Hflap,IPslat,IPflap,
     +                              ISECT,IP,JP)
            CALL WING(QC(M,I,J,1),-QC(M,I,J,2),QC(M,I,J,3),GAMMA,u,v,w,
     +                1.0,M,I,J,DS,ip,jp,CASE,A1,QF,ISECT,
     +                IWING1,IWING2)
c
c           *** Return vortex corners back to its original deflected
c               position. ***
            IF (ASYMDEFL .EQ. -1.0) CALL REVERSE(DELslat,DELflap,
     +                              QF,Hslat,Hflap,IPslat,IPflap,
     +                              ISECT,IP,JP) 
            L = 0
            DO 145, M1 = 1, ISECT
               IF ((M1.GE.IWING1).AND.(M1.LE.IWING2)) GOTO 145
               DO 150, I1 = 1, IP(M1)
                  DO 160, J1 = 1, JP(M1)
                     L = L+1
                     A(K,L) = A(K,L)+ASYM*A1(M1,I1,J1)
 160              CONTINUE
 150           CONTINUE
 145        CONTINUE
            IF (CH.GT.9.0E4) GOTO 170
c           
c           *** Ground effect (mirror image) ***
            CASE = 3.0
            CALL WING(QC(M,I,J,1),QC(M,I,J,2),-2*CH-QC(M,I,J,3),GAMMA,
     +                u,v,w,1.0,M,I,J,DS,ip,jp,CASE,A1,QF,ISECT,
     +                IWING1,IWING2)
            L = 0
            DO 175, M1 = 1, ISECT
               IF ((M1.GE.IWING1).AND.(M1.LE.IWING2)) GOTO 175
               DO 180, I1 = 1, IP(M1)
                  DO 190, J1 = 1, JP(M1)
                     L = L+1
                     A(K,L) = A(K,L)+A1(M1,I1,J1)
 190              CONTINUE
 180           CONTINUE
 175        CONTINUE
c
c           *** Ground effect due to other half of wing ***
            CASE = 4.0
c           *** Reverse flap deflection for antisymmetric deflection ***
            IF (ASYMDEFL .EQ. -1.0) CALL REVERSE(DELslat,DELflap,
     +                              QF,Hslat,Hflap,IPslat,IPflap,
     +                              ISECT,IP,JP)
c
            CALL WING(QC(M,I,J,1),-QC(M,I,J,2),-2*CH-QC(M,I,J,3),GAMMA,
     +                u,v,w,1.0,M,I,J,DS,ip,jp,CASE,A1,QF,ISECT,
     +                IWING1,IWING2)
c
c           *** Return vortex corners back to its original deflected
c               position. ***
            IF (ASYMDEFL .EQ. -1.0) CALL REVERSE(DELslat,DELflap,
     +                              QF,Hslat,Hflap,IPslat,IPflap,
     +                              ISECT,IP,JP) 
            L = 0
            DO 195, M1 = 1, ISECT
               IF ((M1.GE.IWING1).AND.(M1.LE.IWING2)) GOTO 195
               DO 200, I1 = 1, IP(M1)
                  DO 210, J1 = 1, JP(M1)
                     L = L+1
                     A(K,L) = A(K,L)+ASYM*A1(M1,I1,J1)
 210              CONTINUE
 200           CONTINUE
 195        CONTINUE
c
c
c           *** Downwash due to free stream & roll/pitch rate induced
c               velocity components ***
c
c
c           *** If flow is antisymmetric as in the cases of roll rate
c               damping calculations and antisymmetric deflections,
c               the problem becomes solving the antisymmetric vortex
c               distribution.  The symmetric gamma distribution is
c               designated as GAMMAnew and is known.  It will be
c               accounted for on the RHS. ***
c           
  170       SYM(K) = 0.0
            IF (ASYM.EQ.-1.0) THEN

c 
c           *** The original half of wing ***
            CASE = 1.0
            CALL WING(QC(M,I,J,1),QC(M,I,J,2),QC(M,I,J,3),GAMMAnew,
     +                u,v,w,
     +                1.0,M,I,J,DS,ip,jp,CASE,A1,QF,ISECT,
     +                IWING1,IWING2)
            SYM(K) = u*DS(M,I,J,1)+v*DS(M,I,J,2)+w*DS(M,I,J,3)
c
c
c           *** The other half of wing ***
            CASE = 2.0
c           *** Reverse flap deflection for antisymmetric deflection ***
            IF (ASYMDEFL .EQ. -1.0) CALL REVERSE(DELslat,DELflap,
     +                              QF,Hslat,Hflap,IPslat,IPflap,
     +                              ISECT,IP,JP)
c 
           CALL WING(QC(M,I,J,1),-QC(M,I,J,2),QC(M,I,J,3),GAMMAnew,
     +                u,v,w,
     +                1.0,M,I,J,DS,ip,jp,CASE,A1,QF,ISECT,
     +                IWING1,IWING2)
c           *** Return vortex corners back to its original deflected
c               position. ***
            IF (ASYMDEFL .EQ. -1.0) CALL REVERSE(DELslat,DELflap,
     +                              QF,Hslat,Hflap,IPslat,IPflap,
     +                              ISECT,IP,JP) 
c
            SYM(K) = SYM(K)+u*DS(M,I,J,1)-v*DS(M,I,J,2)+w*DS(M,I,J,3)
            IF (CH.GT.100.0) GOTO 215
c           
c           *** Ground effect (mirror image) ***
            CASE = 3.0
            CALL WING(QC(M,I,J,1),QC(M,I,J,2),-2*CH-QC(M,I,J,3),
     +                GAMMAnew,
     +                u,v,w,1.0,M,I,J,DS,ip,jp,CASE,A1,QF,ISECT,
     +                IWING1,IWING2)
            SYM(K) = SYM(K)+u*DS(M,I,J,1)+v*DS(M,I,J,2)-w*DS(M,I,J,3)
c
c           *** Ground effect due to other half of wing ***
            CASE = 4.0
c           *** Reverse flap deflection for antisymmetric deflection ***
            IF (ASYMDEFL .EQ. -1.0) CALL REVERSE(DELslat,DELflap,
     +                              QF,Hslat,Hflap,IPslat,IPflap,
     +                              ISECT,IP,JP)
c
            CALL WING(QC(M,I,J,1),-QC(M,I,J,2),-2*CH-QC(M,I,J,3),
     +                GAMMAnew,
     +                u,v,w,1.0,M,I,J,DS,ip,jp,CASE,A1,QF,ISECT,
     +                IWING1,IWING2)
c
c           *** Return vortex corners back to its original deflected
c               position. ***
            IF (ASYMDEFL .EQ. -1.0) CALL REVERSE(DELslat,DELflap,
     +                              QF,Hslat,Hflap,IPslat,IPflap,
     +                              ISECT,IP,JP) 
c
            SYM(K) = SYM(K)+u*DS(M,I,J,1)-v*DS(M,I,J,2)-w*DS(M,I,J,3)
c
c           *** Done computing symmetrical gamma's contribution to be
c               added to the RHS ***
            ENDIF
c
c
 215        uinf = VF
            vinf = 0.0
            winf = q*(QC(M,I,J,1)-Xcg)-p*(QC(M,I,J,2)-Ycg)
            DW(K) = -(uinf*DS(M,i,j,1)+vinf*DS(M,i,j,2)+winf*
     +              DS(M,i,j,3))-SYM(K)
c
c           *** Total number of panels
            ISUM = K         
 120     CONTINUE
 110  CONTINUE
 105  CONTINUE
C
C     *** check influence coefficient matrix ***
      DO 1001, I = 1,8
         WRITE(15,1002) A(I,1),A(I,2),A(I,3),A(I,4),A(I,5),A(I,6),
     +                  A(I,7),A(I,8),
     +                  A(I,9),DW(I)
 1001 CONTINUE
 1002 FORMAT(10(F8.5,2X))
C
c
c     *** Solve for vortex strengths ***
      write(6,*)'CALLING MATRIX'
      CALL MATRIX(ISUM,A,DW,GAMMA1)
c
c
c     *** Reassign system solution to GAMMA array ***
      K = 0
      DO 225 M = 1, ISECT
         IF ((M.GE.IWING1).AND.(M.LE.IWING2)) GOTO 225
         DO 220, I = 1, IP(M)
            DO 230, J = 1, JP(M)
               K = K+1
               GAMMA(M,I,J) = GAMMA1(K)
 230        CONTINUE
 220     CONTINUE
 225  CONTINUE
c
c
c     *** Determine forces and moment ***
      FL = 0.0
      FD = 0.0
      FM = 0.0
      RM = 0.0
      YM = 0.0

      DO 245 M=1,ISECT
         IF ((M.GE.IWING1).AND.(M.LE.IWING2)) GOTO 245
         DO 240, I = 1, IP(M)
            DO 250, J = 1, JP(M)
               IF (I.EQ.1) GAMMAIJ = GAMMA(M,I,J)
               IF (I.GT.1) GAMMAIJ = GAMMA(M,I,J)-GAMMA(M,I-1,J)
               SPANP = QF(M,I,J+1,2)-QF(M,I,J,2)
               DLIFT = RHO*VF*GAMMAIJ*SPANP
               DSIDEFORCE = -RHO*VF*GAMMAIJ*(QF(M,I,J+1,3)-
     +                      QF(M,I,J,3))
               FL = FL+DLIFT
c
c              *** Induced Drag ***
               CALL WING(QC(M,I,J,1),QC(M,I,J,2),QC(M,I,J,3),GAMMA,
     +                u1,v1,w1,
     +                0.0,M,I,J,DS,ip,jp,CASE,A1,QF,ISECT,
     +                IWING1,IWING2)
               CALL WING(QC(M,I,J,1),-QC(M,I,J,2),QC(M,I,J,3),GAMMA,
     +                u2,v2,w2,
     +                0.0,M,I,J,DS,ip,jp,CASE,A1,QF,ISECT,
     +                IWING1,IWING2)
               IF (CH .GT. 9.0E4) GOTO 260
               CALL WING(QC(M,I,J,1),QC(M,I,J,2),-2*CH-QC(M,I,J,3),
     +                GAMMA,u3,v3,w3,
     +                0.0,M,I,J,DS,ip,jp,CASE,A1,QF,ISECT,
     +                IWING1,IWING2)
               CALL WING(QC(M,I,J,1),-QC(M,I,J,2),-2*CH-QC(M,I,J,3),
     +                GAMMA,u4,v4,w4,
     +                0.0,M,I,J,DS,ip,jp,CASE,A1,QF,ISECT,
     +                IWING1,IWING2)
               GOTO 265
 260           W3 = 0.0
               W4 = 0.0
 265           WIND = W1+W2-W3-W4
               ALPHAI = -WIND/VF
               DDRAG = RHO*SPANP*VF*GAMMAIJ*ALPHAI
               FD = FD+DDRAG
c
c 
c              *** Moments w/r C.G. ***
               xmact   = -((QF(M,I,J,1)+QF(M,I,J+1,1))/2-Xcg*cos(alpha))
               zmact   = -((QF(M,I,J,3)+QF(M,I,J+1,3))/2+Xcg*sin(alpha))
               DMOMENT =  xmact*DLIFT+zmact*ddrag
c              *** Moments for antisymmetric cases ***
c              *** To locate the aerodyanmic center of panel ***
               CLAMBDA = (QF(M,I+1,J+1,1)-QF(M,I,J+1,1))/
     +                   (QF(M,I+1,J,1)-QF(M,I,J,1))
               YARM    = QF(M,I,J,2)+(QF(M,I,J+1,2)-QF(M,I,J,2))/3.0*
     +                   ((1.0+2.0*CLAMBDA)/(1.0+CLAMBDA))-Ycg
               ZARM    = (QF(M,I,J+1,3)+QF(M,I,J,3))/2.0-Zcg
               Rmoment = YARM*DLIFT-ZARM*DSIDEFORCE
               XARM    = (QF(M,I,J+1,1)+QF(M,I,J,1))/2.0-Xcg
               Ymoment = XARM*DSIDEFORCE-YARM*DDRAG         
c
c
               FM = FM+DMOMENT
               RM = RM+RMOMENT
               YM = YM+YMOMENT
 250        CONTINUE
            WRITE(15,1002)GAMMA(M,I,1),GAMMA(M,I,2),GAMMA(M,I,3),
     +      GAMMA(M,I,4),GAMMA(M,I,5),GAMMA(M,I,6),GAMMA(M,I,7),
     +      GAMMA(M,I,8),GAMMA(M,I,9),GAMMA(M,I,10)
 240     CONTINUE
 245  CONTINUE
 2000 format(/4x,'m',4x,'i',4x,'j',3x,'dlift',5x,'xmact',5x,'dmoment',
     1         5x,'xcg',5x,'QF(M,I,J,1)',2x,'QF(M,I,J+1,1)')
 2005 format(3i5,6f12.5)
c
c
      LIFT = FL
      DRAG = FD
      PITCH = FM
c     *** 'roll' is the dimensional roll moment of half of plane (used
c         only for asymmetric loading only. ***
      roll = Rm
c
c     *** 'YAW' is the dimensional yaw moment of half of plane used 
c          only for asymmetric loading only. ***
      YAW = YM
c
c     *** This section is used only to calculate the lift due
c         tail section(s) ***
c
      LIFTTAIL = 0.0
      DO 305, K = 1, ISECT
      IF ((K.GE.IWING1).AND.(K.LE.IWING2)) GOTO 305
c
      DO 300, I = 1, IP(K)
c	   write(*,*) QC(K,I,1,1) 
c	   write(*,*) K,I,ITAIL1,ITAIL2
c	   pause 'press enter'
         IF ((QC(K,I,1,1).GE.QC(ITAIL1,1,1,1)).OR.
     +		(QC(K,I,1,1).GE.QC(ITAIL2,1,1,1))) THEN
         DO 310, J = 1, JP(K)
            IF(I.EQ.1) GAMMAIJ = GAMMA(K,I,J)
            IF(I.GT.1) GAMMAIJ = GAMMA(K,I,J)-GAMMA(K,I-1,J)
            SPANP = QF(K,I,J+1,2)-QF(K,I,J,2)
            LIFTTAIL = LIFTTAIL+RHO*VF*GAMMAIJ*SPANP
 310     CONTINUE
         ENDIF
 300  CONTINUE          
 305  CONTINUE     
c
c
c
      RETURN
      END

c     ********************************************************************
      SUBROUTINE WING(X,Y,Z,GAMMA,U,V,W,ONOFF,M1,I1,J1,
     +                DS,IP,JP,CASE,A1,QF,ISECT,IWING1,IWING2)
      DIMENSION GAMMA(5,8,10),QF(5,10,6,3),A1(5,9,5)
      DIMENSION DS(5,8,5,4),IP(5),JP(5)
c
c     ====================================================================
c     This subroutine determines the induced velocity at (x,y,z) due to 
c     the (semi-span) vortex lattice system.
c     ====================================================================
c
c     *** Initialize some variables ***
      u = 0.0
      v = 0.0
      w = 0.0
c
      DO 50, M = 1, ISECT
         IF ((M.GE.IWING1).AND.(M.LE.IWING2)) GOTO 50      
         DO 100, I = 1, IP(M)+1
            DO 110, J = 1, JP(M)
c              *** For trailing vortices ***
               IF (I.EQ.IP(M)+1) THEN
                  VORTIC = GAMMA(M,IP(M),J)
                  GOTO 120
               ENDIF
c
c
               VORTIC = GAMMA(M,I,J)
               U0 = 0.0
               V0 = 0.0
               W0 = 0.0
c
c              *** Skip for induced drag calc. ***
 120           continue
               U1 = 0.0
               V1 = 0.0
               W1 = 0.0
               U3 = 0.0
               V3 = 0.0
               W3 = 0.0

               IF (ONOFF.EQ.0.0) GOTO 130
c
c              *** Leading segment of vortex ring ***
               CALL VORTEX (X,Y,Z,QF(M,I,J,1),QF(M,I,J,2),QF(M,I,J,3),
     +                      QF(M,I,J+1,1),QF(M,I,J+1,2),QF(M,I,J+1,3),
     +                      VORTIC,U1,V1,W1)
c
c              *** Trailing segment of vortex ring ***  
               CALL VORTEX(X,Y,Z,QF(M,I+1,J+1,1),QF(M,I+1,J+1,2),
     +                     QF(M,I+1,J+1,3),QF(M,I+1,J,1),
     +                     QF(M,I+1,J,2),QF(M,I+1,J,3),
     +                     VORTIC,U3,V3,W3)
c
c              *** Two side segments ***
 130           CALL VORTEX (X,Y,Z,QF(M,I,J+1,1),QF(M,I,J+1,2),
     +                      QF(M,I,J+1,3),QF(M,I+1,J+1,1),
     +                      QF(M,I+1,J+1,2),QF(M,I+1,J+1,3),
     +                      VORTIC,U2,V2,W2)
               CALL VORTEX (X,Y,Z,QF(M,I+1,J,1),QF(M,I+1,J,2),
     +                      QF(M,I+1,J,3),QF(M,I,J,1),
     +                      QF(M,I,J,2),QF(M,I,J,3),
     +                      VORTIC,U4,V4,W4)
c
c           *** Sum vortex-induced velocity components from secgments ***
               U0 = U2+U4+(U1+U3)*ONOFF
               V0 = V2+V4+(V1+V3)*ONOFF
               W0 = W2+W4+(W1+W3)*ONOFF
c
c              *** Skip if Calculating Induced Drag ***
               IF (ONOFF .EQ. 0.0) GOTO 140
c
c              *** For different semi-span vortice lattice systems ***
               IF (CASE .EQ. 1.0)
     +            A1(M,I,J) = U0*DS(M1,I1,J1,1)+V0*DS(M1,I1,J1,2)+
     +                        W0*DS(M1,I1,J1,3)
               IF (CASE .EQ. 2.0)
     +            A1(M,I,J) = U0*DS(M1,I1,J1,1)-V0*DS(M1,I1,J1,2)+
     +                        W0*DS(M1,I1,J1,3)
               IF (CASE .EQ. 3.0)
     +            A1(M,I,J) = U0*DS(M1,I1,J1,1)+V0*DS(M1,I1,J1,2)-
     +                        W0*DS(M1,I1,J1,3)
               IF (CASE .EQ. 4.0)
     +            A1(M,I,J) = U0*DS(M1,I1,J1,1)-V0*DS(M1,I1,J1,2)-
     +                        W0*DS(M1,I1,J1,3)

c              *** Adding trailing (free) vortex influence to T.E. vortex 
c                  ring since they have the same strength ***
 140           IF (I.EQ.IP(M)+1) A1(M,IP(M),J) = A1(M,IP(M),J)+
     +                                           A1(M,IP(M)+1,J)
               U = U+U0
               V = V+V0
               W = W+W0
c
c              
 110        CONTINUE
 100     CONTINUE
 50   CONTINUE
c
c
c
      RETURN
      END


c     ********************************************************************
      SUBROUTINE VORTEX(X,Y,Z,X1,Y1,Z1,X2,Y2,Z2,GAMMA,U,V,W)
c
c     ====================================================================
c     This subroutine determines the induced velocity at (X,Y,Z) due to 
c     the vortex segment from (X1,Y1,Z1) to (X2,Y2,Z2).
c     ====================================================================
c
      common /consts/ PI

C     pi = 4.0*ATAN(1.0)
c     *** R1 is the vector from (X1,Y1,Z1) to (X,Y,Z) ***
c     *** R2 is the vector from (X2,Y2,Z2) to (X,Y,Z) ***
c     *** R1 cross R2 ***
      R1R2X  =   (Y-Y1)*(Z-Z2) - (Z-Z1)*(Y-Y2)
      R1R2Y  = -((X-X1)*(Z-Z2) - (X-X2)*(Z-Z1))
      R1R2Z  =   (X-X1)*(Y-Y2) - (X-X2)*(Y-Y1)
c
c     *** Magnitude (squared) of vectors ***
      SQUARE = (R1R2X**2) + (R1R2Y**2) + (R1R2Z**2)
      R1     = SQRT((X-X1)**2 + (Y-Y1)**2 + (Z-Z1)**2)
      R2     = SQRT((X-X2)**2 + (Y-Y2)**2 + (Z-Z2)**2)
c
c     *** Dot products ***
      R0R1   = (X2-X1)*(X-X1)+(Y2-Y1)*(Y-Y1)+(Z2-Z1)*(Z-Z1)
      R0R2   = (X2-X1)*(X-X2)+(Y2-Y1)*(Y-Y2)+(Z2-Z1)*(Z-Z2) 
c
      COEFF  = GAMMA/(4.0*PI*SQUARE)*(R0R1/R1 - R0R2/R2)
      U      = R1R2X*COEFF
      V      = R1R2Y*COEFF
      W      = R1R2Z*COEFF    
C
      RETURN
      END

c     *******************************************************************
      SUBROUTINE MATRIX (N,A,D,C)
      INTEGER N
      DIMENSION A(200,200),C(200),D(200)

C      NOTE:  N=SIZE OF MATRIX
C      USER MUST UPDATE THE VALUE OF N IN THE ABOVE PARAMETER STATEMENT
C      INPUT MATRIX [A] & [D]; OUTPUT MATRIX [C] SUCH THAT [A][C]=[D]

C      STARTING MATRIX INVERSION

      DO 10, I = 1,N
 10   C(I)     = 0.0

      DO 1100 K= 1,N - 1
      M = K + 1
      L = K
 1130 Q = ABS(A(M,K)) - ABS(A(L,K))
      IF (Q.GT.0.0) L = M
      IF (M.LT.N)   THEN
                    M=M+1
                    GO TO 1130
                    ENDIF
      IF (L.EQ.K) GO TO 1110

      DO 1200 J = K, N
      DU        = A(K,J)
      A(K,J)    = A(L,J)
      A(L,J)    = DU
 1200 CONTINUE

      DD        = D(K)
      D(K)      = D(L)
      D(L)      = DD

 1110 M         = K + 1
 1120 Q         = A(M,K)/A(K,K)

      A(M,K)    = 0
      DO 1300 J = K+1, N
      A(M,J)    = A(M,J)-Q*A(K,J)
 1300 CONTINUE

      D(M)      = D(M) - Q*D(K)
      IF (M.LT.N) THEN
                  M = M+1
                  GO TO 1120
                  ENDIF
 1100 CONTINUE
      C(N)      = D(N)/A(N,N)
      DO 1500 M = N-1,1,-1
      Q         = 0
      DO 1400 J = M+1,N
      Q         = Q + A(M,J)*C(J)
      C(M)      = (D(M)-Q)/A(M,M)
 1400 CONTINUE
 1500 CONTINUE

       END


c     *****************************************************************
      SUBROUTINE REVERSE(DELslat,DELflap,VC,Hslat,Hflap,IPslat,
     +                   IPflap,ISECT,IP,JP)
      REAL DELslat(5),DELflap(5),VC(5,10,6,3),Hslat(5,2,3),Hflap(5,2,3)
      INTEGER IPslat(5),IPflap(5),ISECT,IP(5),JP(5)
c
c     ==================================================================
c     This subroutine is to reverse the control surfaces' deflection
c     in the opposite direction for antisymmetric deflection cases.
c     Calling this subroutine the second time will deflect the affected
c     control surface(s) back to the original deflection angle.
c     ==================================================================
c
c     *** Find the affected vortex corner points ***
      DO 100, M = 1, ISECT
c
c        *** If the deflected is (are) slat(s) ***
         IF (DELslat(M) .NE. 0.0) THEN
            DO 110, J = 1, JP(M)+1
               DO 120, I = 1, IPslat(M)
                  CALL DEFLECT(Hslat(M,1,1),Hslat(M,1,2),Hslat(M,1,3),
     +                         Hslat(M,2,1),Hslat(M,2,2),Hslat(M,2,3),
     +                         VC(M,I,J,1),VC(M,I,J,2),VC(M,I,J,3),
     +                         -2.0*DELslat(M))
 120           CONTINUE
 110        CONTINUE
            DELslat(M) = -DELslat(M)
         ENDIF
c
c        *** For deflected flaps ***
         IF (DELflap(M) .NE. 0.0) THEN
            DO 130, J = 1, JP(M)+1
               DO 140, I = IP(M)+2-IPFLAP(M), IP(M)+1
                  CALL DEFLECT(Hflap(M,1,1),Hflap(M,1,2),Hflap(M,1,3),
     +                         Hflap(M,2,1),Hflap(M,2,2),Hflap(M,2,3),
     +                         VC(M,I,J,1),VC(M,I,J,2),VC(M,I,J,3),
     +                         -2.0*DELflap(M))
 140           CONTINUE
 130        CONTINUE
            DELflap(M) = -DELflap(M)
         ENDIF
 100  CONTINUE
      RETURN
      END
