PROGRAM WingPlanAnal

c    REM ***************************************************************
C        *                                                             *
c    REM * wing geomtry program, from:                                 *
c        *                                                             *
c    REM *          GEOMETRY FOR AERODYNAMICS, BY W.H. MASON           *
c    REM *           AEROCAL, BOX 799, HUNTINGTON, NY  11743           *
c    REM *        now: Box 11850, Blacksburg, VA 24062-1850 and        *
c    REM *        Department of Aerospace and Ocean Engineering        *
c    REM *        Virginia Tech, Blacksburg, VA, 24062                 *
c    REM *               contact: whmason@vt.edu                       *
c    REM *            see also: http://www.aoe.vt.edu/                 *
c    REM *                                                             *
c        *                calculator program: 1982                     *
c    REM *     APPLESOFT VERSION FOR APPLE ][ SERIES COMPUTERS         *
c    REM *    RELEASE 1.2            LAST MOD: JANUARY 4, 1987         *
c    REM *                                                             *
c    REM *  PORTED WITH MINIMUM MODS TO MAC COMPUTER FOR USE           *
c    REM * WITH QuickBASIC and no print out capability - FEBRUARY 1990 *
c    REM *     last mod: June 20, 1993 (minor for web, Oct 16, 1997    *
c        *                                                             *
c        *      fortran, December 2005                                 *
c        *                                                             *
c    REM ***************************************************************

      DIMENSION YL(101),XL(101),YT(101),XT(101)
      DIMENSION TITLE(15)
      character*1  ans
      character*4  TITLE
      character*15 datset
      IREAD  = 1

      write(6,100)
  100 format(/5x,'WingPlanAnal'//5x,'Wing Planform Geometry Analysis'/
     1  5x,'Virginia Tech Aircraft Design Software Series'/
     2  5x,'W.H. Mason, Department of Aerospace and Ocean Engineering'/
     3  5x,'Virginia Tech, Blacksburg, VA 24061, email: whmason@vt.edu'/
     3  5x,'version: January 22, 2006'/) 

      Write(6,101)
  101 Format(5x,'Planform Properties')

      write(6,5)
      read (5,7) datset
      OPEN (unit = IREAD,file=datset)

    5 FORMAT(' Enter name of data set:')
    7 FORMAT(a15)
   10 FORMAT(15A4)

      READ (IREAD,10) TITLE
      write(6,102) TITLE
  102 format(/5x,'Input Case Title: ',15A4/)

      read(iread,15) FILM
   15 FORMAT(4F10.5,I2)
      iLM = FILM
      read(iread,7)

      write(6,103)
  103 format(5x,'Planform Points')

      write(6,104)  iLM
  104 format(/5x,'Leading Edge',5x,'iLM =',i3//
     1        7x,'i',8x,'YLE',12x,'XLE')

      do 110 ii = 1, iLM
      read(iread,15) yl(ii),xl(ii)
      write(6,106) ii,yl(ii),xl(ii)
  110 continue
  106 format(5x,i3,4x,f10.4,4x,f10.4)

      read(iread,15) FITM
      iTM = FITM
      read(iread,7)

      write(6,108)  iTM
  108 format(/5x,'Trailing Edge',5x,'iTM =',i3//
     1        7x,'i',8x,'YTE',12x,'XTE')

      do 120 ii = 1,iTM
      read(iread,15) yt(ii),xt(ii)
      write(6,106) ii,yt(ii),xt(ii)
  120 continue

      BB = yt(iTM)
      IF (yl(iLM) .gt. BB ) BB = yl(iLM)

      xi1 = 0
      xI2 = 0
      xI3 = 0
      xi4 = 0
      xi5 = 0

      N  = 20
      DS = BB / N
      IL = 1
      IT = 1


      write(6,140)
  140 format(/5x,'Interpolated LE and TE points and sweep'//
     1 2x,' i',8x,'eta',8x,'y',7x,'XLE',2x,'LE sweep(deg)',2x,'XTE',
     2 2x,'TE sweep(deg)')

      Ncount = N + 1
      do 200 icount = 1,Ncount
      ii = icount - 1
      IL = 1
      IT = 1
      Q = 1
      IF (ii .eq. 0 .OR. ii .eq. N) Q = .5

      y   = ii * DS
      eta = y/bb

      call XLEPTS(yl,xl,y,sl,xf,ilm,il)
      call XTEPTS(yt,xt,y,st,xb,itm,it)

      sldeg = 57.2958*atan(sl)
      stdeg = 57.2958*atan(st)

      write(6,404) ii,eta,y,xf,sldeg,xb,stdeg
  404 format(2x,i2,2x,6f10.3)

      C = xb - xf

      IF (ii .eq. 0) CR = C
      IF (ii .eq. N) CT = C

      xi1 = xi1 + Q * C
      xI2 = xI2 + Q * C**2
      xI3 = xI3 + Q * (C * xf + C**2 / 2)
      xi4 = xi4 + Q * y * C
      xi5 = xi5 + Q * C * xf
  200 continue

      xi1 = 2 * DS * xi1
      Q   = 2 * DS / xi1
      xI2 = Q * xI2
      xI3 = Q * xI3
      xi4 = Q * xi4
      xi5 = Q * xi5
      xi6 = xi5 + xi2/4.0
      AR  = (2 * y)**2 / xi1
      CA  = xi1 / (2 * y)
      TP  = CT / CR

      write(6,142)
  142 format(/5x,'Integral Quantities')

      write(6,160) xi1,xi2,xi3,xi4,xi5,xi6,ar,ca,tp
  160 format(/5x,'Planform Area            =  ',f10.5/
     1        5x,'Mean Aerodynamic Chord   =  ',f10.5/
     2        5x,'X-Centroid               =  ',f10.5/
     3        5x,'Spanwise position of MAC =  ',f10.5/
     4        5x,'X-Leading Edge of MAC    =  ',f10.5/
     5        5x,'Quarter Chord of MAC     =  ',f10.5/
     6        5x,'Aspect Ratio             =  ',f10.5/
     7        5x,'Average Chord            =  ',f10.5/
     8        5x,'Taper Ratio              =  ',f10.5/
     9       //)
 
c     REM   * POINT BY POINT INTERPOLATION, WITH SWEEP
   
      write(6,502)
  502 format(/5x,'Do you want LE/TE and Sweep? (Y/N) ')
      read(5,503) ans
  503 format(a1)
      if( ans .eq. 'N' .OR. ans .eq. 'n') go to 999

c       PRINT: NA$ = " Local Planform Properties":GOSUB 5000

  504 continue
      write(6,506) 
  506 format(/5x,'Input span location, Y= ')
      read(5,508) y
  508 format(f10.5)

      eta = y/bb

      call XLEPTS(yl,xl,y,sl,xf,ilm,il)
      call XTEPTS(yt,xt,y,st,xb,itm,it)

      sldeg = 57.2958*atan(sl)
      stdeg = 57.2958*atan(st)

      chord = xb - xf

      write(6,510) y,eta,xf,sldeg,xb,stdeg,chord
  510 format(/5x,'y   = 'f10.5,5x,'eta   = ',f10.5/
     1 5x,'X-le= ',f10.5,5x,'sweep = ',f10.5/
     2 5x,'X-te= ',f10.5,5x,'sweep = ',f10.5/
     3 5x,'c   = ',f10.5)

      write(6,512)
  512 format(/5x,'Do you want another station? (Y/N) ')
      read(5,503) ans

      IF (ANS .eq. 'Y' .or. ANS .eq. 'y') GO TO 504

  999 continue 

      stop
      end

      subroutine XLEPTS(yl,xl,y,sl,xf,ilm,il)
c     REM  * LINEAR INTERPOLATION FOR LEADING EDGE

c     yl,xl - input array of pts defining the LE
c     ilm   - the number of yl,xl pairs
c     y     - input value for interpolation
c     xf    - output value of Xle
c     sl    - output: tan of the sweep of the LE
c     il    - interval integer input (=1)

      dimension yl(101), xl(101)

  720 D2 = y - yl(IL)
      IF (D2 .gt. 0.0) go to 736
      IF (D2 .eq. 0.0) go to 742
  726 IF (IL .eq. 1)   go to 742
      IL = IL - 1
      D2 = y - yl(IL)
      IF (D2 .lt. 0.0) go to 726
      IF (D2 .ge. 0.0) go to 742
  734 IL = IL + 1
      D2 = D1
  736 IF (IL .eq. iLM - 1) go to 742
      D1 = y - yl(IL + 1)
      IF (D1 .ge. 0.0) go to 734

  742 SL = (xl(IL + 1) - xl(IL)) / (yl(IL + 1) - yl(IL))
      xf = xl(IL) + SL * (y - yl(IL))

      RETURN
      end

      subroutine XTEPTS(yt,xt,y,st,xb,itm,it) 
c     REM  * LINEAR INTERPOLATION FOR TRAILING EDGE

c     yt,xt - input array of pts defining the TE
c     itm   - the number of yt,xt pairs
c     y     - input value for interpolation
c     xb    - output value of Xte
c     st    - output: tan of the sweep of the TE
c     it    - interval integer input (=1)
 
      dimension yt(101), xt(101)
      eps = 1.0e-4

      if (y .lt. yt(1)) write(6,100) y,yt(1)
  100 format(/5x,'y lt yt(1)',3x,2e12.7)

      if (y .gt. yt(itm)) write(6,102) y,yt(itm)
  102 format(/5x,'y gt yt(itm)',3x,2e12.7)

  760 D4 = y - yt(IT)
      IF (D4 .gt. 0.0) go to 776
      IF (abs(D4) .lt. eps) go to 782
  766 IF (IT .eq. 1)   go to 782
      IT = IT - 1
      D4 = y - yt(IT)
      IF (D4 .lt. 0.0) go to 766
      IF (D4 .ge. 0.0) go to 782
  774 IT = IT + 1
      D4 = D3
  776 IF (IT .eq. iTM - 1) go to 782
      D3 = y - yt(IT + 1)
      IF (D3 .ge. 0) go to 774

  782 ST = (xt(IT + 1) - xt(IT)) / (yt(IT + 1) - yt(IT))
      xb = xt(IT) + ST * (y - yt(IT))

      RETURN
      end