aboutsummaryrefslogtreecommitdiff
path: root/src/xfoil.f
diff options
context:
space:
mode:
Diffstat (limited to 'src/xfoil.f')
-rw-r--r--src/xfoil.f2580
1 files changed, 2580 insertions, 0 deletions
diff --git a/src/xfoil.f b/src/xfoil.f
new file mode 100644
index 0000000..73cf11e
--- /dev/null
+++ b/src/xfoil.f
@@ -0,0 +1,2580 @@
+C***********************************************************************
+C Module: xfoil.f
+C
+C Copyright (C) 2000 Mark Drela
+C
+C This program is free software; you can redistribute it and/or modify
+C it under the terms of the GNU General Public License as published by
+C the Free Software Foundation; either version 2 of the License, or
+C (at your option) any later version.
+C
+C This program is distributed in the hope that it will be useful,
+C but WITHOUT ANY WARRANTY; without even the implied warranty of
+C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+C GNU General Public License for more details.
+C
+C You should have received a copy of the GNU General Public License
+C along with this program; if not, write to the Free Software
+C Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+C***********************************************************************
+C
+ PROGRAM XFOIL
+C--- Uncomment for Win32/Compaq Visual Fortran compiler (needed for GETARG)
+ccc USE DFLIB
+C
+ INCLUDE 'XFOIL.INC'
+ CHARACTER*4 COMAND
+ CHARACTER*128 COMARG, PROMPT
+ CHARACTER*1 ANS
+C
+ DIMENSION IINPUT(20)
+ DIMENSION RINPUT(20)
+ LOGICAL ERROR
+C
+C---- max panel angle threshold for warning
+ DATA ANGTOL / 40.0 /
+C
+ VERSION = 6.97
+ WRITE(*,1005) VERSION
+ 1005 FORMAT(
+ & /' ==================================================='
+ & /' XFOIL Version', F5.2
+ & /' Copyright (C) 2000 Mark Drela, Harold Youngren'
+ & //' This software comes with ABSOLUTELY NO WARRANTY,'
+ & /' subject to the GNU General Public License.'
+ & //' Caveat computor'
+ & /' ===================================================')
+C
+ CALL INIT
+ LU = 8
+ CALL GETDEF(LU,'xfoil.def', .TRUE.)
+C
+C---- try to read airfoil from command line argument, if any
+ FNAME = ' '
+ NARG = IARGC()
+ IF(NARG.GT.0) CALL GETARG(NARG,FNAME)
+C
+ IF(FNAME(1:1) .NE. ' ') THEN
+ CALL LOAD(FNAME,ITYPE)
+C
+ IF(ITYPE.GT.0 .AND. NB.GT.0) THEN
+ccc CALL PANGEN(.TRUE.)
+ CALL ABCOPY(.TRUE.)
+C
+ CALL CANG(X,Y,N,0, IMAX,AMAX)
+ IF(ABS(AMAX).GT.ANGTOL) THEN
+ WRITE(*,1081) AMAX, IMAX
+ 1081 FORMAT(
+ & /' WARNING: Poor input coordinate distribution'
+ & /' Excessive panel angle', F7.1,' at i =', I4
+ & /' Repaneling with PANE and/or PPAR suggested'
+ & /' (doing GDES,CADD before repaneling _may_'
+ & /' improve excessively coarse LE spacing' )
+ CALL PANPLT
+ ENDIF
+ ENDIF
+ ENDIF
+C
+ WRITE(*,1100) XCMREF,YCMREF,NPAN
+ 1100 FORMAT(
+ & /' QUIT Exit program'
+ & //' .OPER Direct operating point(s)'
+ & /' .MDES Complex mapping design routine'
+ & /' .QDES Surface speed design routine'
+ & /' .GDES Geometry design routine'
+ & //' SAVE f Write airfoil to labeled coordinate file'
+ & /' PSAV f Write airfoil to plain coordinate file'
+ & /' ISAV f Write airfoil to ISES coordinate file'
+ & /' MSAV f Write airfoil to MSES coordinate file'
+ & /' REVE Reverse written-airfoil node ordering'
+ & /' DELI i Change written-airfoil file delimiters'
+ & //' LOAD f Read buffer airfoil from coordinate file'
+ & /' NACA i Set NACA 4,5-digit airfoil and buffer airfoil'
+ & /' INTE Set buffer airfoil by interpolating two airfoils'
+ & /' NORM Buffer airfoil normalization toggle'
+ & /' HALF Halve the number of points in buffer airfoil'
+ & /' XYCM rr Change CM reference location, currently ',2F8.5
+ & //' BEND Display structural properties of current airfoil'
+ & //' PCOP Set current-airfoil panel nodes directly',
+ & ' from buffer airfoil points'
+ & /' PANE Set current-airfoil panel nodes (',I4,' )',
+ & ' based on curvature'
+ & /' .PPAR Show/change paneling'
+ & //' .PLOP Plotting options'
+ & //' WDEF f Write current-settings file'
+ & /' RDEF f Reread current-settings file'
+ & /' NAME s Specify new airfoil name'
+ & /' NINC Increment name version number'
+ & //' Z Zoom | (available in all menus)'
+ & /' U Unzoom | ')
+C
+C---- start of menu loop
+ 500 CONTINUE
+ CALL ASKC(' XFOIL^',COMAND,COMARG)
+C
+C---- get command line numeric arguments, if any
+ DO I=1, 20
+ IINPUT(I) = 0
+ RINPUT(I) = 0.0
+ ENDDO
+ NINPUT = 0
+ CALL GETINT(COMARG,IINPUT,NINPUT,ERROR)
+ NINPUT = 0
+ CALL GETFLT(COMARG,RINPUT,NINPUT,ERROR)
+C
+C===============================================
+ IF(COMAND.EQ.' ') THEN
+ GO TO 500
+C
+C===============================================
+ ELSEIF(COMAND.EQ.'? ') THEN
+ WRITE(*,1100) XCMREF, YCMREF, NPAN
+C
+C===============================================
+ ELSEIF(COMAND.EQ.'QUIT') THEN
+ CALL PLCLOSE
+ STOP
+C
+C===============================================
+ ELSEIF(COMAND.EQ.'OPER') THEN
+ CALL OPER
+C
+C===============================================
+ ELSEIF(COMAND.EQ.'MDES') THEN
+ CALL MDES
+C
+C===============================================
+ ELSEIF(COMAND.EQ.'QDES') THEN
+ CALL QDES
+C
+C===============================================
+ ELSEIF(COMAND.EQ.'GDES') THEN
+ CALL GDES
+C
+C===============================================
+ ELSEIF(COMAND.EQ.'SAVE') THEN
+ CALL SAVE(1,COMARG)
+C
+C===============================================
+ ELSEIF(COMAND.EQ.'PSAV') THEN
+ CALL SAVE(0,COMARG)
+C
+C===============================================
+ ELSEIF(COMAND.EQ.'USAV') THEN
+ CALL SAVE(-1,COMARG)
+C
+C===============================================
+ ELSEIF(COMAND.EQ.'ISAV') THEN
+ CALL SAVE(2,COMARG)
+C
+C===============================================
+ ELSEIF(COMAND.EQ.'MSAV') THEN
+ CALL MSAVE(COMARG)
+C
+C===============================================
+ ELSEIF(COMAND.EQ.'REVE') THEN
+ LCLOCK = .NOT.LCLOCK
+ IF(LCLOCK) THEN
+ WRITE(*,*) 'Airfoil will be written in clockwise order'
+ ELSE
+ WRITE(*,*) 'Airfoil will be written in counterclockwise order'
+ ENDIF
+C
+C===============================================
+ ELSEIF(COMAND.EQ.'DELI') THEN
+ 40 CONTINUE
+ IF(NINPUT.GE.1) THEN
+ KDNEW = IINPUT(1)
+ ELSE
+ WRITE(*,2100) KDELIM
+ 2100 FORMAT(/' --------------------------'
+ & /' 0 blank'
+ & /' 1 comma'
+ & /' 2 tab',
+ & //' currently, delimiter =', I2)
+ CALL ASKI('Enter new delimiter',KDNEW)
+ ENDIF
+C
+ IF(KDNEW.LT.0 .OR. KDNEW.GT.2) THEN
+ NINPUT = 0
+ GO TO 40
+ ELSE
+ KDELIM = KDNEW
+ ENDIF
+C
+C===============================================
+ ELSEIF(COMAND.EQ.'LOAD') THEN
+ CALL LOAD(COMARG,ITYPE)
+ IF(ITYPE.GT.0 .AND. NB.GT.0) THEN
+ccc CALL PANGEN(.TRUE.)
+ CALL ABCOPY(.TRUE.)
+C
+ CALL CANG(X,Y,N,0, IMAX,AMAX)
+ IF(ABS(AMAX).GT.ANGTOL) THEN
+ WRITE(*,1081) AMAX, IMAX
+ CALL PANPLT
+ ENDIF
+ ENDIF
+C
+C===============================================
+ ELSEIF(COMAND.EQ.'NACA') THEN
+ CALL NACA(IINPUT(1))
+C
+C===============================================
+ ELSEIF(COMAND.EQ.'INTE') THEN
+ CALL INTE
+C
+C===============================================
+ ELSEIF(COMAND.EQ.'INTX') THEN
+ CALL INTX
+C
+C===============================================
+ ELSEIF(COMAND.EQ.'NORM') THEN
+ LNORM = .NOT.LNORM
+ IF(LNORM) THEN
+ WRITE(*,*) 'Loaded airfoil will be normalized'
+ ELSE
+ WRITE(*,*) 'Loaded airfoil won''t be normalized'
+ ENDIF
+C
+C===============================================
+ ELSEIF(COMAND.EQ.'HALF') THEN
+ CALL HALF(XB,YB,SB,NB)
+ CALL SCALC(XB,YB,SB,NB)
+ CALL SEGSPL(XB,XBP,SB,NB)
+ CALL SEGSPL(YB,YBP,SB,NB)
+C
+ CALL GEOPAR(XB,XBP,YB,YBP,SB,NB, W1,
+ & SBLE,CHORDB,AREAB,RADBLE,ANGBTE,
+ & EI11BA,EI22BA,APX1BA,APX2BA,
+ & EI11BT,EI22BT,APX1BT,APX2BT,
+ & THICKB,CAMBRB )
+C
+C==========================================
+ ELSEIF(COMAND.EQ.'XYCM') THEN
+ IF(NINPUT.GE.2) THEN
+ XCMREF = RINPUT(1)
+ YCMREF = RINPUT(2)
+ ELSE
+ CALL ASKR('Enter new CM reference X^',XCMREF)
+ CALL ASKR('Enter new CM reference Y^',YCMREF)
+ ENDIF
+C
+C===============================================
+ ELSEIF(COMAND.EQ.'BEND') THEN
+ IF(N.EQ.0) THEN
+ WRITE(*,*)
+ WRITE(*,*) ' No airfoil available'
+ GO TO 500
+ ENDIF
+C
+ CALL BENDUMP(N,X,Y)
+C
+C===============================================
+ ELSEIF(COMAND.EQ.'BENP') THEN
+ IF(N.EQ.0) THEN
+ WRITE(*,*)
+ WRITE(*,*) ' No airfoil available'
+ GO TO 500
+ ENDIF
+C
+ DO I = 1, N
+ W1(I) = 1.0
+ ENDDO
+ CALL BENDUMP2(N,X,Y,W1)
+C
+C===============================================
+ ELSEIF(COMAND.EQ.'PCOP') THEN
+ CALL ABCOPY(.TRUE.)
+ccc CALL PANPLT
+C
+C===============================================
+ ELSEIF(COMAND.EQ.'PANE') THEN
+ CALL PANGEN(.TRUE.)
+ccc CALL PANPLT
+C
+C===============================================
+ ELSEIF(COMAND.EQ.'PPAR') THEN
+ CALL GETPAN
+C
+C===============================================
+ ELSEIF(COMAND.EQ.'PLOP') THEN
+ CALL OPLSET(IDEV,IDEVRP,IPSLU,
+ & SIZE,PLOTAR,
+ & XMARG,YMARG,XPAGE,YPAGE,
+ & CH,SCRNFR,LCURS,LLAND)
+C
+C===============================================
+ ELSEIF(COMAND.EQ.'WDEF') THEN
+ LU = 8
+ IF(COMARG(1:1).EQ.' ') THEN
+ FNAME = 'xfoil.def'
+ ELSE
+ FNAME = COMARG
+ ENDIF
+ CALL STRIP(FNAME,NFN)
+ OPEN(LU,FILE=FNAME,STATUS='OLD',ERR=703)
+ WRITE(*,701) FNAME(1:NFN)
+ 701 FORMAT(/' File ', A, ' exists. Overwrite? Y')
+ READ(*,1000) ANS
+ IF(INDEX('Nn',ANS).EQ.0) GO TO 706
+ WRITE(*,*)
+ WRITE(*,*) 'No action taken'
+ CLOSE(LU)
+C
+ 703 OPEN(LU,FILE=FNAME,STATUS='UNKNOWN')
+ 706 CALL WRTDEF(LU)
+ WRITE(*,708) FNAME(1:NFN)
+ 708 FORMAT(/' File ', A, ' written')
+ CLOSE(LU)
+C
+C===============================================
+ ELSEIF(COMAND.EQ.'RDEF') THEN
+ IF(COMARG(1:1).EQ.' ') THEN
+ FNAME = 'xfoil.def'
+ ELSE
+ FNAME = COMARG
+ ENDIF
+C
+ LU = 8
+ CALL GETDEF(LU,FNAME, .FALSE.)
+C
+C===============================================
+ ELSEIF(COMAND.EQ.'NAME') THEN
+ IF(COMARG.EQ.' ') THEN
+ CALL NAMMOD(NAME,0,-1)
+ ELSE
+ NAME = COMARG
+ ENDIF
+ CALL STRIP(NAME,NNAME)
+C
+C===============================================
+ ELSEIF(COMAND.EQ.'NINC') THEN
+ CALL NAMMOD(NAME,1,1)
+ CALL STRIP(NAME,NNAME)
+C
+C===============================================
+ ELSEIF(COMAND.EQ.'Z ') THEN
+ IF(LPLOT) THEN
+ CALL USETZOOM(.TRUE.,.TRUE.)
+ CALL REPLOT(IDEV)
+ ENDIF
+C
+C===============================================
+ ELSEIF(COMAND.EQ.'U ') THEN
+ IF(LPLOT) THEN
+ CALL CLRZOOM
+ CALL REPLOT(IDEV)
+ ENDIF
+C
+C===============================================
+ ELSE
+ WRITE(*,1050) COMAND
+ 1050 FORMAT(1X,A4,' command not recognized. Type a "?" for list')
+C
+ ENDIF
+C
+C===============================================
+ GO TO 500
+C
+ 1000 FORMAT(A)
+ END ! XFOIL
+
+
+ SUBROUTINE INIT
+C---------------------------------------------------
+C Variable initialization/default routine.
+C See file XFOIL.INC for variable description.
+C---------------------------------------------------
+ INCLUDE 'XFOIL.INC'
+C
+ PI = 4.0*ATAN(1.0)
+ HOPI = 0.50/PI
+ QOPI = 0.25/PI
+ DTOR = PI/180.0
+C
+C---- default Cp/Cv (air)
+ GAMMA = 1.4
+ GAMM1 = GAMMA - 1.0
+C
+C---- set unity freestream speed
+ QINF = 1.0
+C
+C---- initialize freestream Mach number to zero
+ MATYP = 1
+ MINF1 = 0.
+C
+ ALFA = 0.0
+ COSA = 1.0
+ SINA = 0.0
+C
+ DO 10 I=1, IQX
+ GAMU(I,1) = 0.
+ GAMU(I,2) = 0.
+ GAM(I) = 0.
+ GAM_A(I) = 0.
+ 10 CONTINUE
+ PSIO = 0.
+C
+ CL = 0.
+ CM = 0.
+ CD = 0.
+C
+ SIGTE = 0.0
+ GAMTE = 0.0
+ SIGTE_A = 0.
+ GAMTE_A = 0.
+C
+ DO 20 I=1, IZX
+ SIG(I) = 0.
+ 20 CONTINUE
+C
+ NQSP = 0
+ DO 30 K=1, IPX
+ ALQSP(K) = 0.
+ CLQSP(K) = 0.
+ CMQSP(K) = 0.
+ DO 302 I=1, IBX
+ QSPEC(I,K) = 0.
+ 302 CONTINUE
+ 30 CONTINUE
+C
+ AWAKE = 0.0
+ AVISC = 0.0
+C
+ KIMAGE = 1
+ YIMAGE = -10.0
+ LIMAGE = .FALSE.
+C
+C---- output coordinate file delimiters
+C- KDELIM = 0 blanks
+C- = 1 commas
+C- = 2 tabs
+ KDELIM = 0
+C
+ LGAMU = .FALSE.
+ LQINU = .FALSE.
+ LVISC = .FALSE.
+ LWAKE = .FALSE.
+ LPACC = .FALSE.
+ LBLINI = .FALSE.
+ LIPAN = .FALSE.
+ LQAIJ = .FALSE.
+ LADIJ = .FALSE.
+ LWDIJ = .FALSE.
+ LCPXX = .FALSE.
+ LQVDES = .FALSE.
+ LQSPEC = .FALSE.
+ LQREFL = .FALSE.
+ LVCONV = .FALSE.
+ LCPREF = .FALSE.
+ LFOREF = .FALSE.
+ LPFILE = .FALSE.
+ LPFILX = .FALSE.
+ LPPSHO = .FALSE.
+ LBFLAP = .FALSE.
+ LFLAP = .FALSE.
+ LEIW = .FALSE.
+ LSCINI = .FALSE.
+ LPLOT = .FALSE.
+ LCLIP = .FALSE.
+ LVLAB = .TRUE.
+ LCMINP = .FALSE.
+ LHMOMP = .FALSE.
+ LFREQP = .TRUE.
+C
+ LCURS = .TRUE.
+ LLAND = .TRUE.
+ LGSAME = .FALSE.
+C
+ LGPARM = .TRUE.
+ LPLCAM = .FALSE.
+C
+C---- input airfoil will not be normalized
+ LNORM = .FALSE.
+C
+C---- airfoil will not be forced symmetric
+ LQSYM = .FALSE.
+ LGSYM = .FALSE.
+C
+C---- endpoint slopes will be matched
+ LQSLOP = .TRUE.
+ LGSLOP = .TRUE.
+ LCSLOP = .TRUE.
+C
+C---- grids on Qspec(s) and buffer airfoil geometry plots will be plotted
+ LQGRID = .TRUE.
+ LGGRID = .TRUE.
+ LGTICK = .TRUE.
+C
+C---- no grid on Cp plots
+ LCPGRD = .FALSE.
+C
+C---- grid and no symbols are to be used on BL variable plots
+ LBLGRD = .TRUE.
+ LBLSYM = .FALSE.
+C
+C---- buffer and current airfoil flap hinge coordinates
+ XBF = 0.0
+ YBF = 0.0
+ XOF = 0.0
+ YOF = 0.0
+C
+ NCPREF = 0
+C n
+C---- circle plane array size (257, or largest 2 + 1 that will fit array size)
+ ANN = LOG(FLOAT((2*IQX)-1))/LOG(2.0)
+ NN = INT( ANN + 0.00001 )
+ NC1 = 2**NN + 1
+ NC1 = MIN( NC1 , 257 )
+C
+C---- default paneling parameters
+ NPAN = 160
+ CVPAR = 1.0
+ CTERAT = 0.15
+ CTRRAT = 0.2
+C
+C---- default paneling refinement zone x/c endpoints
+ XSREF1 = 1.0
+ XSREF2 = 1.0
+ XPREF1 = 1.0
+ XPREF2 = 1.0
+C
+C---- no polars present to begin with
+ NPOL = 0
+ IPACT = 0
+ DO IP = 1, NPX
+ PFNAME(IP) = ' '
+ PFNAMX(IP) = ' '
+ ENDDO
+C
+C---- no reference polars
+ NPOLREF = 0
+C
+C---- plot aspect ratio, character size
+ PLOTAR = 0.55
+ CH = 0.015
+C
+C---- airfoil node tick-mark size (as fraction of arc length)
+ GTICK = 0.0005
+C
+C---- Cp limits in Cp vs x plot
+ CPMAX = 1.0
+ CPMIN = -2.0
+ CPDEL = -0.5
+ PFAC = PLOTAR/(CPMAX-CPMIN)
+C
+C---- Ue limits in Ue vs x plot
+ UEMAX = 1.8
+ UEMIN = -1.0
+ UEDEL = 0.2
+ UFAC = PLOTAR/(UEMAX-UEMIN)
+C
+C---- DCp limits in CAMB loading plot
+ YPMIN = -0.6
+ YPMAX = 0.6
+C
+C---- scaling factor for Cp vector plot
+ VFAC = 0.25
+C
+C---- offsets and scale factor for airfoil in Cp vs x plot
+ XOFAIR = 0.09
+ YOFAIR = -.01
+ FACAIR = 0.70
+C
+C---- u/Qinf scale factor for profile plotting
+ UPRWT = 0.02
+C
+C---- polar plot options, grid, list, legend, no CDW
+ LPGRID = .TRUE.
+ LPCDW = .FALSE.
+ LPLIST = .TRUE.
+ LPLEGN = .TRUE.
+ LAECEN = .FALSE.
+ LPCDH = .FALSE.
+ LPCMDOT = .FALSE.
+C
+C---- axis limits and annotation deltas for polar plot
+ CPOLPLF(1,ICD) = 0.0
+ CPOLPLF(2,ICD) = 0.04
+ CPOLPLF(3,ICD) = 0.01
+C
+ CPOLPLF(1,ICL) = 0.
+ CPOLPLF(2,ICL) = 1.5
+ CPOLPLF(3,ICL) = 0.5
+C
+ CPOLPLF(1,ICM) = -0.25
+ CPOLPLF(2,ICM) = 0.0
+ CPOLPLF(3,ICM) = 0.05
+C
+ CPOLPLF(1,IAL) = -4.0
+ CPOLPLF(2,IAL) = 10.0
+ CPOLPLF(3,IAL) = 2.0
+C
+C---- widths of plot boxes in polar plot page
+ XCDWID = 0.45
+ XALWID = 0.25
+ XOCWID = 0.20
+C
+C---- line style and color index for each polar
+C
+C 1 ***************************** SOLID
+C 2 **** **** **** **** **** **** LONG DASHED
+C 3 ** ** ** ** ** ** ** ** ** ** SHORT DASHED
+C 4 * * * * * * * * * * * * * * * DOTTED
+C 5 ***** * ***** * ***** * ***** DASH-DOT
+C 6 ***** * * ***** * * ***** * * DASH-DOT-DOT
+C 7 ***** * * * ***** * * * ***** DASH-DOT-DOT-DOT
+C 8 **** **** * * **** **** * * DASH-DASH-DOT-DOT
+C
+C 3 red
+C 4 orange
+C 5 yellow
+C 6 green
+C 7 cyan
+C 8 blue
+C 9 violet
+C 10 magenta
+C
+ DO IP=1, NPX
+cc ILINP(IP) = 1 + MOD(IP-1,8)
+cc ICOLP(IP) = 3 + MOD(IP-1,8)
+C
+C------ normally solid, going to dashed after IP=7
+ ILINP(IP) = 1 + (IP-1)/7
+C
+C------ skip yellow (hard to see on white background)
+ ICOLP(IP) = 3 + MOD(IP-1,7)
+ IF(ICOLP(IP) .GE. 5) ICOLP(IP) = ICOLP(IP) + 1
+ ENDDO
+C
+C---- polar variables to be written to polar save file
+ IPOL(1) = IAL
+ IPOL(2) = ICL
+ IPOL(3) = ICD
+ IPOL(4) = ICP
+ IPOL(5) = ICM
+ NIPOL = 5
+ NIPOL0 = 5
+C
+ JPOL(1) = JTN
+ NJPOL = 1
+C
+C---- default Cm reference location
+ XCMREF = 0.25
+ YCMREF = 0.
+C
+C---- default viscous parameters
+ RETYP = 1
+ REINF1 = 0.
+ ACRIT = 9.0
+ XSTRIP(1) = 1.0
+ XSTRIP(2) = 1.0
+ XOCTR(1) = 1.0
+ XOCTR(2) = 1.0
+ YOCTR(1) = 0.
+ YOCTR(2) = 0.
+ WAKLEN = 1.0
+C
+ IDAMP = 0
+C
+C---- set BL calibration parameters
+ CALL BLPINI
+C
+C---- Newton iteration limit
+ ITMAX = 20
+C
+C---- max number of unconverged sequence points for early exit
+ NSEQEX = 4
+C
+C---- drop tolerance for BL system solver
+ VACCEL = 0.01
+C
+C---- inverse-mapping auto-filter level
+ FFILT = 0.0
+C
+C---- default overlay airfoil filename
+ ONAME = ' '
+C
+C---- default filename prefix
+ PREFIX = ' '
+C
+C---- Plotting flag
+ IDEV = 1 ! X11 window only
+c IDEV = 2 ! B&W PostScript output file only (no color)
+c IDEV = 3 ! both X11 and B&W PostScript file
+c IDEV = 4 ! Color PostScript output file only
+c IDEV = 5 ! both X11 and Color PostScript file
+C
+C---- Re-plotting flag (for hardcopy)
+ IDEVRP = 2 ! B&W PostScript
+c IDEVRP = 4 ! Color PostScript
+C
+C---- PostScript output logical unit and file specification
+ IPSLU = 0 ! output to file plot.ps on LU 4 (default case)
+c IPSLU = ? ! output to file plot?.ps on LU 10+?
+C
+C---- screen fraction taken up by plot window upon opening
+ SCRNFR = 0.80
+C
+C---- Default plot size in inches
+C- (Default plot window is 11.0 x 8.5)
+C- (Must be smaller than XPAGE if objects are to fit on paper page)
+ SIZE = 10.0
+
+C---- plot-window dimensions in inches for plot blowup calculations
+C- currently, 11.0 x 8.5 default window is hard-wired in libPlt
+ XPAGE = 11.0
+ YPAGE = 8.5
+C
+C---- page margins in inches
+ XMARG = 0.0
+ YMARG = 0.0
+C
+C---- set top and bottom-side colors
+cc ICOLS(1) = 5
+cc ICOLS(2) = 7
+ ICOLS(1) = 8
+ ICOLS(2) = 3
+C
+C 3 red
+C 4 orange
+C 5 yellow
+C 6 green
+C 7 cyan
+C 8 blue
+C 9 violet
+C 10 magenta
+C
+C
+ CALL PLINITIALIZE
+C
+C---- set up color spectrum
+ NCOLOR = 64
+ CALL COLORSPECTRUMHUES(NCOLOR,'RYGCBM')
+C
+C
+ NNAME = 32
+ NAME = ' '
+CCC 12345678901234567890123456789012
+C
+C---- MSES domain parameters (not used in XFOIL)
+ ISPARS = ' -2.0 3.0 -2.5 3.5'
+C
+C---- set MINF, REINF, based on current CL-dependence
+ CALL MRCL(1.0,MINF_CL,REINF_CL)
+C
+C---- set various compressibility parameters from MINF
+ CALL COMSET
+C
+ RETURN
+ END ! INIT
+
+
+ SUBROUTINE MRCL(CLS,M_CLS,R_CLS)
+C-------------------------------------------
+C Sets actual Mach, Reynolds numbers
+C from unit-CL values and specified CLS
+C depending on MATYP,RETYP flags.
+C-------------------------------------------
+ INCLUDE 'XFOIL.INC'
+ REAL M_CLS
+C
+ CLA = MAX( CLS , 0.000001 )
+C
+ IF(RETYP.LT.1 .OR. RETYP.GT.3) THEN
+ WRITE(*,*) 'MRCL: Illegal Re(CL) dependence trigger.'
+ WRITE(*,*) ' Setting fixed Re.'
+ RETYP = 1
+ ENDIF
+ IF(MATYP.LT.1 .OR. MATYP.GT.3) THEN
+ WRITE(*,*) 'MRCL: Illegal Mach(CL) dependence trigger.'
+ WRITE(*,*) ' Setting fixed Mach.'
+ MATYP = 1
+ ENDIF
+C
+C
+ IF(MATYP.EQ.1) THEN
+C
+ MINF = MINF1
+ M_CLS = 0.
+C
+ ELSE IF(MATYP.EQ.2) THEN
+C
+ MINF = MINF1/SQRT(CLA)
+ M_CLS = -0.5*MINF/CLA
+C
+ ELSE IF(MATYP.EQ.3) THEN
+C
+ MINF = MINF1
+ M_CLS = 0.
+C
+ ENDIF
+C
+C
+ IF(RETYP.EQ.1) THEN
+C
+ REINF = REINF1
+ R_CLS = 0.
+C
+ ELSE IF(RETYP.EQ.2) THEN
+C
+ REINF = REINF1/SQRT(CLA)
+ R_CLS = -0.5*REINF/CLA
+C
+ ELSE IF(RETYP.EQ.3) THEN
+C
+ REINF = REINF1/CLA
+ R_CLS = -REINF /CLA
+C
+ ENDIF
+C
+C
+ IF(MINF .GE. 0.99) THEN
+ WRITE(*,*)
+ WRITE(*,*) 'MRCL: CL too low for chosen Mach(CL) dependence'
+ WRITE(*,*) ' Aritificially limiting Mach to 0.99'
+ MINF = 0.99
+ M_CLS = 0.
+ ENDIF
+C
+ RRAT = 1.0
+ IF(REINF1 .GT. 0.0) RRAT = REINF/REINF1
+C
+ IF(RRAT .GT. 100.0) THEN
+ WRITE(*,*)
+ WRITE(*,*) 'MRCL: CL too low for chosen Re(CL) dependence'
+ WRITE(*,*) ' Aritificially limiting Re to ',REINF1*100.0
+ REINF = REINF1*100.0
+ R_CLS = 0.
+ ENDIF
+C
+ RETURN
+ END ! MRCL
+
+
+
+ SUBROUTINE GETDEF(LU,FILNAM,LASK)
+ CHARACTER*(*) FILNAM
+ LOGICAL LASK
+C-----------------------------------------------------
+C Reads in default parameters from file xfoil.def
+C If LASK=t, ask user if file is to be read.
+C-----------------------------------------------------
+ INCLUDE 'XFOIL.INC'
+ LOGICAL LCOLOR
+ CHARACTER*1 ANS
+C
+ 1000 FORMAT(A)
+C
+ OPEN(LU,FILE=FILNAM,STATUS='OLD',ERR=90)
+ IF(LASK) THEN
+ WRITE(*,1050) FILNAM
+ 1050 FORMAT(/' Read settings from file ', A, ' ? Y')
+ READ(*,1000) ANS
+ IF(INDEX('Nn',ANS).NE.0) THEN
+ CLOSE(LU)
+ RETURN
+ ENDIF
+ ENDIF
+C
+ CLMIN = CPOLPLF(1,ICL)
+ CLMAX = CPOLPLF(2,ICL)
+ CLDEL = CPOLPLF(3,ICL)
+C
+ CDMIN = CPOLPLF(1,ICD)
+ CDMAX = CPOLPLF(2,ICD)
+ CDDEL = CPOLPLF(3,ICD)
+C
+ ALMIN = CPOLPLF(1,IAL)
+ ALMAX = CPOLPLF(2,IAL)
+ ALDEL = CPOLPLF(3,IAL)
+C
+ CMMIN = CPOLPLF(1,ICM)
+ CMMAX = CPOLPLF(2,ICM)
+ CMDEL = CPOLPLF(3,ICM)
+C
+C---- default paneling parameters (viscous)
+ READ(LU,*,ERR=80) NPAN, CVPAR, CTERAT, CTRRAT
+ READ(LU,*,ERR=80) XSREF1, XSREF2, XPREF1, XPREF2
+C
+C---- plotting parameters
+ READ(LU,*,ERR=80) SIZE, PLOTAR, CH, SCRNFR
+C
+C---- plot sizes
+ READ(LU,*,ERR=80) XPAGE, YPAGE, XMARG, YMARG
+C
+C---- plot flags
+ READ(LU,*,ERR=80) LCOLOR, LCURS
+C
+C---- Cp limits in Cp vs x plot
+ READ(LU,*,ERR=80) CPMAX, CPMIN, CPDEL
+ PFAC = PLOTAR/(CPMAX-CPMIN)
+C
+C---- airfoil x-offset and scale factor in Cp vs x plot, BL profile weight
+ READ(LU,*,ERR=80) XOFAIR, FACAIR, UPRWT
+C
+C---- polar plot CL,CD,alpha,CM min,max,delta
+ READ(LU,*,ERR=80) (CPOLPLF(K,ICL), K=1, 3)
+ READ(LU,*,ERR=80) (CPOLPLF(K,ICD), K=1, 3)
+ READ(LU,*,ERR=80) (CPOLPLF(K,IAL), K=1, 3)
+ READ(LU,*,ERR=80) (CPOLPLF(K,ICM), K=1, 3)
+C
+C---- default Mach and viscous parameters
+ READ(LU,*,ERR=80) MATYP, MINF1, VACCEL
+ READ(LU,*,ERR=80) RETYP, RMILL, ACRIT
+ READ(LU,*,ERR=80) XSTRIP(1), XSTRIP(2)
+C
+ IF( LCOLOR) IDEVRP = 4
+ IF(.NOT.LCOLOR) IDEVRP = 2
+C
+ REINF1 = RMILL * 1.0E6
+C
+C---- set MINF, REINF
+ CALL MRCL(1.0,MINF_CL,REINF_CL)
+C
+C---- set various compressibility parameters from new MINF
+ CALL COMSET
+C
+ CLOSE(LU)
+ WRITE(*,1600) FILNAM
+ 1600 FORMAT(/' Default parameters read in from file ', A,':' /)
+ CALL WRTDEF(6)
+ RETURN
+C
+ 80 CONTINUE
+ CLOSE(LU)
+ WRITE(*,1800) FILNAM
+ 1800 FORMAT(/' File ', A,' read error'
+ & /' Settings may have been changed')
+ RETURN
+C
+ 90 CONTINUE
+ WRITE(*,1900) FILNAM
+ 1900 FORMAT(/' File ', A,' not found')
+ RETURN
+C
+ END ! GETDEF
+
+
+
+ SUBROUTINE WRTDEF(LU)
+C------------------------------------------
+C Writes default parameters to unit LU
+C------------------------------------------
+ INCLUDE 'XFOIL.INC'
+ LOGICAL LCOLOR
+C
+ LCOLOR = IDEVRP.EQ.4
+C
+C---- default paneling parameters (viscous)
+ WRITE(LU,1010) NPAN , CVPAR , CTERAT, CTRRAT
+ WRITE(LU,1020) XSREF1, XSREF2, XPREF1, XPREF2
+C
+C---- plotting parameters
+ WRITE(LU,1030) SIZE, PLOTAR, CH, SCRNFR
+C
+C---- plot sizes
+ WRITE(LU,1032) XPAGE, YPAGE, XMARG, YMARG
+C
+C---- plot flags
+ WRITE(LU,1034) LCOLOR, LCURS
+C
+C---- Cp limits in Cp vs x plot
+ WRITE(LU,1040) CPMAX, CPMIN, CPDEL
+C
+C---- x-offset and scale factor for airfoil on Cp vs x plot
+ WRITE(LU,1050) XOFAIR, FACAIR, UPRWT
+C
+C---- polar plot CL,CD,alpha,CM min,max,delta
+ WRITE(LU,1061) (CPOLPLF(K,ICL), K=1, 3)
+ WRITE(LU,1062) (CPOLPLF(K,ICD), K=1, 3)
+ WRITE(LU,1063) (CPOLPLF(K,IAL), K=1, 3)
+ WRITE(LU,1064) (CPOLPLF(K,ICM), K=1, 3)
+C
+C---- default viscous parameters
+ WRITE(LU,1071) MATYP , MINF1 , VACCEL
+ WRITE(LU,1072) RETYP , REINF1/1.0E6 , ACRIT
+ WRITE(LU,1080) XSTRIP(1), XSTRIP(2)
+C
+ RETURN
+C...............................................
+ 1010 FORMAT(1X,I5,4X,F9.4,F9.4,F9.4,' | Npan PPanel TErat REFrat')
+ 1020 FORMAT(1X,F9.4 ,F9.4,F9.4,F9.4,' | XrefS1 XrefS2 XrefP1 XrefP2')
+ 1030 FORMAT(1X,F9.4 ,F9.4,F9.4,F9.4,' | Size plotAR CHsize ScrnFr')
+ 1032 FORMAT(1X,F9.4 ,F9.4,F9.4,F9.4,' | Xpage Ypage Xmargn Ymargn')
+ 1034 FORMAT(1X,L2,7X,L2,7X,9X , 9X ,' | Lcolor Lcursor' )
+ 1040 FORMAT(1X,F9.4 ,F9.4,F9.4, 9X ,' | CPmax CPmin CPdel' )
+ 1050 FORMAT(1X,F9.4 ,F9.4,F9.4, 9X ,' | XoffAir ScalAir BLUwt' )
+ 1061 FORMAT(1X,F9.4 ,F9.4,F9.4, 9X ,' | CLmin CLmax CLdel' )
+ 1062 FORMAT(1X,F9.4 ,F9.4,F9.4, 9X ,' | CDmin CDmax CDdel' )
+ 1063 FORMAT(1X,F9.4 ,F9.4,F9.4, 9X ,' | ALmin ALmax ALdel' )
+ 1064 FORMAT(1X,F9.4 ,F9.4,F9.4, 9X ,' | CMmin CMmax CMdel' )
+ 1071 FORMAT(1X,I3,6X,F9.4,F9.4, 9X ,' | MAtype Mach Vaccel' )
+ 1072 FORMAT(1X,I3,6X,F9.4,F9.4, 9X ,' | REtype Re/10^6 Ncrit' )
+ 1080 FORMAT(1X,F9.4 ,F9.4, 9X , 9X ,' | XtripT XtripB' )
+ END ! WRTDEF
+
+
+ SUBROUTINE COMSET
+ INCLUDE 'XFOIL.INC'
+C
+C---- set Karman-Tsien parameter TKLAM
+ BETA = SQRT(1.0 - MINF**2)
+ BETA_MSQ = -0.5/BETA
+C
+ TKLAM = MINF**2 / (1.0 + BETA)**2
+ TKL_MSQ = 1.0 / (1.0 + BETA)**2
+ & - 2.0*TKLAM/ (1.0 + BETA) * BETA_MSQ
+C
+C---- set sonic Pressure coefficient and speed
+ IF(MINF.EQ.0.0) THEN
+ CPSTAR = -999.0
+ QSTAR = 999.0
+ ELSE
+ CPSTAR = 2.0 / (GAMMA*MINF**2)
+ & * (( (1.0 + 0.5*GAMM1*MINF**2)
+ & /(1.0 + 0.5*GAMM1 ))**(GAMMA/GAMM1) - 1.0)
+ QSTAR = QINF/MINF
+ & * SQRT( (1.0 + 0.5*GAMM1*MINF**2)
+ & /(1.0 + 0.5*GAMM1 ) )
+ ENDIF
+C
+ RETURN
+ END ! COMSET
+
+
+ SUBROUTINE CPCALC(N,Q,QINF,MINF,CP)
+C---------------------------------------------
+C Sets compressible Cp from speed.
+C---------------------------------------------
+ DIMENSION Q(N),CP(N)
+ REAL MINF
+C
+ LOGICAL DENNEG
+C
+ BETA = SQRT(1.0 - MINF**2)
+ BFAC = 0.5*MINF**2 / (1.0 + BETA)
+C
+ DENNEG = .FALSE.
+C
+ DO 20 I=1, N
+ CPINC = 1.0 - (Q(I)/QINF)**2
+ DEN = BETA + BFAC*CPINC
+ CP(I) = CPINC / DEN
+ IF(DEN .LE. 0.0) DENNEG = .TRUE.
+ 20 CONTINUE
+C
+ IF(DENNEG) THEN
+ WRITE(*,*)
+ WRITE(*,*) 'CPCALC: Local speed too large. ',
+ & 'Compressibility corrections invalid.'
+ ENDIF
+C
+ RETURN
+ END ! CPCALC
+
+
+ SUBROUTINE CLCALC(N,X,Y,GAM,GAM_A,ALFA,MINF,QINF,
+ & XREF,YREF,
+ & CL,CM,CDP, CL_ALF,CL_MSQ)
+C-----------------------------------------------------------
+C Integrates surface pressures to get CL and CM.
+C Integrates skin friction to get CDF.
+C Calculates dCL/dAlpha for prescribed-CL routines.
+C-----------------------------------------------------------
+ DIMENSION X(N),Y(N), GAM(N), GAM_A(N)
+ REAL MINF
+C
+ccC---- moment-reference coordinates
+cc XREF = 0.25
+cc YREF = 0.
+C
+ SA = SIN(ALFA)
+ CA = COS(ALFA)
+C
+ BETA = SQRT(1.0 - MINF**2)
+ BETA_MSQ = -0.5/BETA
+C
+ BFAC = 0.5*MINF**2 / (1.0 + BETA)
+ BFAC_MSQ = 0.5 / (1.0 + BETA)
+ & - BFAC / (1.0 + BETA) * BETA_MSQ
+C
+ CL = 0.0
+ CM = 0.0
+
+ CDP = 0.0
+C
+ CL_ALF = 0.
+ CL_MSQ = 0.
+C
+ I = 1
+ CGINC = 1.0 - (GAM(I)/QINF)**2
+ CPG1 = CGINC/(BETA + BFAC*CGINC)
+ CPG1_MSQ = -CPG1/(BETA + BFAC*CGINC)*(BETA_MSQ + BFAC_MSQ*CGINC)
+C
+ CPI_GAM = -2.0*GAM(I)/QINF**2
+ CPC_CPI = (1.0 - BFAC*CPG1)/ (BETA + BFAC*CGINC)
+ CPG1_ALF = CPC_CPI*CPI_GAM*GAM_A(I)
+C
+ DO 10 I=1, N
+ IP = I+1
+ IF(I.EQ.N) IP = 1
+C
+ CGINC = 1.0 - (GAM(IP)/QINF)**2
+ CPG2 = CGINC/(BETA + BFAC*CGINC)
+ CPG2_MSQ = -CPG2/(BETA + BFAC*CGINC)*(BETA_MSQ + BFAC_MSQ*CGINC)
+C
+ CPI_GAM = -2.0*GAM(IP)/QINF**2
+ CPC_CPI = (1.0 - BFAC*CPG2)/ (BETA + BFAC*CGINC)
+ CPG2_ALF = CPC_CPI*CPI_GAM*GAM_A(IP)
+C
+ DX = (X(IP) - X(I))*CA + (Y(IP) - Y(I))*SA
+ DY = (Y(IP) - Y(I))*CA - (X(IP) - X(I))*SA
+ DG = CPG2 - CPG1
+C
+ AX = (0.5*(X(IP)+X(I))-XREF)*CA + (0.5*(Y(IP)+Y(I))-YREF)*SA
+ AY = (0.5*(Y(IP)+Y(I))-YREF)*CA - (0.5*(X(IP)+X(I))-XREF)*SA
+ AG = 0.5*(CPG2 + CPG1)
+C
+ DX_ALF = -(X(IP) - X(I))*SA + (Y(IP) - Y(I))*CA
+ AG_ALF = 0.5*(CPG2_ALF + CPG1_ALF)
+ AG_MSQ = 0.5*(CPG2_MSQ + CPG1_MSQ)
+C
+ CL = CL + DX* AG
+ CDP = CDP - DY* AG
+ CM = CM - DX*(AG*AX + DG*DX/12.0)
+ & - DY*(AG*AY + DG*DY/12.0)
+C
+ CL_ALF = CL_ALF + DX*AG_ALF + AG*DX_ALF
+ CL_MSQ = CL_MSQ + DX*AG_MSQ
+C
+ CPG1 = CPG2
+ CPG1_ALF = CPG2_ALF
+ CPG1_MSQ = CPG2_MSQ
+ 10 CONTINUE
+C
+ RETURN
+ END ! CLCALC
+
+
+
+ SUBROUTINE CDCALC
+ INCLUDE 'XFOIL.INC'
+C
+ SA = SIN(ALFA)
+ CA = COS(ALFA)
+C
+ IF(LVISC .AND. LBLINI) THEN
+C
+C----- set variables at the end of the wake
+ THWAKE = THET(NBL(2),2)
+ URAT = UEDG(NBL(2),2)/QINF
+ UEWAKE = UEDG(NBL(2),2) * (1.0-TKLAM) / (1.0 - TKLAM*URAT**2)
+ SHWAKE = DSTR(NBL(2),2)/THET(NBL(2),2)
+C
+C----- extrapolate wake to downstream infinity using Squire-Young relation
+C (reduces errors of the wake not being long enough)
+ CD = 2.0*THWAKE * (UEWAKE/QINF)**(0.5*(5.0+SHWAKE))
+C
+ ELSE
+C
+ CD = 0.0
+C
+ ENDIF
+C
+C---- calculate friction drag coefficient
+ CDF = 0.0
+ DO 20 IS=1, 2
+ DO 205 IBL=3, IBLTE(IS)
+ I = IPAN(IBL ,IS)
+ IM = IPAN(IBL-1,IS)
+ DX = (X(I) - X(IM))*CA + (Y(I) - Y(IM))*SA
+ CDF = CDF + 0.5*(TAU(IBL,IS)+TAU(IBL-1,IS))*DX * 2.0/QINF**2
+ 205 CONTINUE
+ 20 CONTINUE
+C
+ RETURN
+ END ! CDCALC
+
+
+
+ SUBROUTINE LOAD(FILNAM,ITYPE)
+C------------------------------------------------------
+C Reads airfoil file into buffer airfoil
+C and does various initial processesing on it.
+C------------------------------------------------------
+ INCLUDE 'XFOIL.INC'
+ CHARACTER*(*) FILNAM
+C
+ FNAME = FILNAM
+ IF(FNAME(1:1) .EQ. ' ') CALL ASKS('Enter filename^',FNAME)
+C
+ LU = 9
+ CALL AREAD(LU,FNAME,IBX,XB,YB,NB,NAME,ISPARS,ITYPE,1)
+ IF(ITYPE.EQ.0) RETURN
+C
+ IF(ITYPE.EQ.1) CALL ASKS('Enter airfoil name^',NAME)
+ CALL STRIP(NAME,NNAME)
+C
+C---- set default prefix for other filenames
+ KDOT = INDEX(FNAME,'.')
+ IF(KDOT.EQ.0) THEN
+ PREFIX = FNAME
+ ELSE
+ PREFIX = FNAME(1:KDOT-1)
+ ENDIF
+ CALL STRIP(PREFIX,NPREFIX)
+C
+C---- calculate airfoil area assuming counterclockwise ordering
+ AREA = 0.0
+ DO 50 I=1, NB
+ IP = I+1
+ IF(I.EQ.NB) IP = 1
+ AREA = AREA + 0.5*(YB(I)+YB(IP))*(XB(I)-XB(IP))
+ 50 CONTINUE
+C
+ IF(AREA.GE.0.0) THEN
+ LCLOCK = .FALSE.
+ WRITE(*,1010) NB
+ ELSE
+C----- if area is negative (clockwise order), reverse coordinate order
+ LCLOCK = .TRUE.
+ WRITE(*,1011) NB
+ DO 55 I=1, NB/2
+ XTMP = XB(NB-I+1)
+ YTMP = YB(NB-I+1)
+ XB(NB-I+1) = XB(I)
+ YB(NB-I+1) = YB(I)
+ XB(I) = XTMP
+ YB(I) = YTMP
+ 55 CONTINUE
+ ENDIF
+C
+ IF(LNORM) THEN
+ CALL NORM(XB,XBP,YB,YBP,SB,NB)
+ WRITE(*,1020)
+ ENDIF
+C
+ CALL SCALC(XB,YB,SB,NB)
+ CALL SEGSPL(XB,XBP,SB,NB)
+ CALL SEGSPL(YB,YBP,SB,NB)
+C
+ CALL GEOPAR(XB,XBP,YB,YBP,SB,NB, W1,
+ & SBLE,CHORDB,AREAB,RADBLE,ANGBTE,
+ & EI11BA,EI22BA,APX1BA,APX2BA,
+ & EI11BT,EI22BT,APX1BT,APX2BT,
+ & THICKB,CAMBRB )
+C
+ XBLE = SEVAL(SBLE,XB,XBP,SB,NB)
+ YBLE = SEVAL(SBLE,YB,YBP,SB,NB)
+ XBTE = 0.5*(XB(1) + XB(NB))
+ YBTE = 0.5*(YB(1) + YB(NB))
+C
+ WRITE(*,1050) XBLE,YBLE, CHORDB,
+ & XBTE,YBTE
+C
+C---- set reasonable MSES domain parameters for non-MSES coordinate file
+ IF(ITYPE.LE.2) THEN
+ XBLE = SEVAL(SBLE,XB,XBP,SB,NB)
+ YBLE = SEVAL(SBLE,YB,YBP,SB,NB)
+ XINL = XBLE - 2.0*CHORDB
+ XOUT = XBLE + 3.0*CHORDB
+ YBOT = YBLE - 2.5*CHORDB
+ YTOP = YBLE + 3.5*CHORDB
+ XINL = AINT(20.0*ABS(XINL/CHORDB)+0.5)/20.0 * SIGN(CHORDB,XINL)
+ XOUT = AINT(20.0*ABS(XOUT/CHORDB)+0.5)/20.0 * SIGN(CHORDB,XOUT)
+ YBOT = AINT(20.0*ABS(YBOT/CHORDB)+0.5)/20.0 * SIGN(CHORDB,YBOT)
+ YTOP = AINT(20.0*ABS(YTOP/CHORDB)+0.5)/20.0 * SIGN(CHORDB,YTOP)
+ WRITE(ISPARS,1005) XINL, XOUT, YBOT, YTOP
+ 1005 FORMAT(1X, 4F8.2 )
+ ENDIF
+C
+C---- wipe out old flap hinge location
+ XBF = 0.0
+ YBF = 0.0
+ LBFLAP = .FALSE.
+C
+C---- wipe out off-design alphas, CLs
+cc NALOFF = 0
+cc NCLOFF = 0
+C
+ RETURN
+C...............................................................
+ 1010 FORMAT(/' Number of input coordinate points:', I4
+ & /' Counterclockwise ordering')
+ 1011 FORMAT(/' Number of input coordinate points:', I4
+ & /' Clockwise ordering')
+ 1020 FORMAT(/' Airfoil has been normalized')
+ 1050 FORMAT(/' LE x,y =', 2F10.5,' | Chord =',F10.5
+ & /' TE x,y =', 2F10.5,' |' )
+ END ! LOAD
+
+
+
+ SUBROUTINE SAVE(IFTYP,FNAME1)
+C--------------------------------
+C Writes out current airfoil
+C--------------------------------
+ INCLUDE 'XFOIL.INC'
+ CHARACTER*(*) FNAME1
+C
+ CHARACTER*1 ANS, DELIM
+ CHARACTER*128 LINE
+C
+ IF (KDELIM.EQ.0) THEN
+ DELIM = ' '
+ ELSEIF(KDELIM.EQ.1) THEN
+ DELIM = ','
+ ELSEIF(KDELIM.EQ.2) THEN
+ DELIM = CHAR(9)
+ ELSE
+ WRITE(*,*) '? Illegal delimiter. Using blank.'
+ DELIM = ' '
+ ENDIF
+C
+C
+ LU = 2
+C
+C---- get output filename if it was not supplied
+ IF(FNAME1(1:1) .NE. ' ') THEN
+ FNAME = FNAME1
+ ELSE
+ CALL ASKS('Enter output filename^',FNAME)
+ ENDIF
+C
+ OPEN(LU,FILE=FNAME,STATUS='OLD',ERR=5)
+ WRITE(*,*)
+ WRITE(*,*) 'Output file exists. Overwrite? Y'
+ READ(*,1000) ANS
+ IF(INDEX('Nn',ANS).EQ.0) GO TO 6
+C
+ CLOSE(LU)
+ WRITE(*,*) 'Current airfoil not saved.'
+ RETURN
+C
+ 5 OPEN(LU,FILE=FNAME,STATUS='NEW',ERR=90)
+ 6 REWIND(LU)
+C
+ IF(IFTYP.GE.1) THEN
+C----- write name to first line
+ WRITE(LU,1000) NAME(1:NNAME)
+ ENDIF
+C
+ IF(IFTYP.GE.2) THEN
+C----- write MSES domain parameters to second line
+ DO K=80, 1, -1
+ IF(INDEX(ISPARS(K:K),' ') .NE. 1) GO TO 11
+ ENDDO
+ 11 CONTINUE
+C
+ WRITE(LU,1000) ISPARS(1:K)
+ ENDIF
+C
+ IF(LCLOCK) THEN
+C----- write out in clockwise order (reversed from internal XFOIL order)
+ IBEG = N
+ IEND = 1
+ INCR = -1
+ ELSE
+C----- write out in counterclockwise order (same as internal XFOIL order)
+ IBEG = 1
+ IEND = N
+ INCR = 1
+ ENDIF
+C
+ IF(IFTYP.EQ.-1) THEN
+ DO I = IBEG, IEND, INCR
+ WRITE(LU,1400) INT(X(I)+SIGN(0.5,X(I))),
+ & INT(Y(I)+SIGN(0.5,Y(I)))
+ ENDDO
+C
+ ELSE
+ DO I = IBEG, IEND, INCR
+ IF(KDELIM .EQ. 0) THEN
+ WRITE(LU,1100) X(I), Y(I)
+C
+ ELSE
+ WRITE(LINE,1200) X(I), DELIM, Y(I)
+ CALL BSTRIP(LINE,NLINE)
+ WRITE(LU,1000) LINE(1:NLINE)
+ ENDIF
+ ENDDO
+ ENDIF
+C
+ CLOSE(LU)
+ RETURN
+C
+ 90 WRITE(*,*) 'Bad filename.'
+ WRITE(*,*) 'Current airfoil not saved.'
+ RETURN
+C
+ 1000 FORMAT(A)
+ 1100 FORMAT(1X,G15.7, G15.7)
+ 1200 FORMAT(1X,F10.6, A, F10.6)
+ 1400 FORMAT(1X, I12, I12)
+ END ! SAVE
+
+
+
+ SUBROUTINE MSAVE(FNAME1)
+C------------------------------------------
+C Writes out current airfoil as one
+C element in a multielement MSES file.
+C------------------------------------------
+ INCLUDE 'XFOIL.INC'
+ CHARACTER*(*) FNAME1
+C
+ CHARACTER*80 NAME1, ISPARS1
+C
+ PARAMETER (NEX=5)
+ DIMENSION NTMP(NEX)
+ DIMENSION XTMP(2*IQX,NEX), YTMP(2*IQX,NEX)
+ EQUIVALENCE (Q(1,1),XTMP(1,1)), (Q(1,IQX/2),YTMP(1,1))
+C
+ LU = 2
+C
+C---- get output filename if it was not supplied
+ IF(FNAME1(1:1) .NE. ' ') THEN
+ FNAME = FNAME1
+ ELSE
+ CALL ASKS('Enter output filename for element replacement^',FNAME)
+ ENDIF
+C
+ OPEN(LU,FILE=FNAME,STATUS='OLD',ERR=9005)
+C
+ READ(LU,1000,ERR=9010) NAME1
+ READ(LU,1000,ERR=9010) ISPARS1
+C
+ DO NN1=80, 2, -1
+ IF(NAME1(NN1:NN1) .NE. ' ') GO TO 10
+ ENDDO
+ 10 CONTINUE
+C
+ DO NI1=80, 2, -1
+ IF(ISPARS1(NI1:NI1) .NE. ' ') GO TO 20
+ ENDDO
+ 20 CONTINUE
+C
+C---- read in existing airfoil coordinates
+ 40 DO 55 IEL=1, NEX
+ DO 50 I=1, 2*IQX+1
+ READ(LU,*,END=56) XTMP(I,IEL), YTMP(I,IEL)
+ IF(XTMP(I,IEL).EQ.999.0) THEN
+ NTMP(IEL) = I-1
+ GO TO 55
+ ENDIF
+ 50 CONTINUE
+ STOP 'LOAD: Array overflow'
+ 55 CONTINUE
+ NEL = NEX
+C
+ 56 IF(I.EQ.1) THEN
+C----- coordinate file has "999.0 999.0" at the end ...
+ NEL = IEL-1
+ ELSE
+C----- coordinate file has no ending line
+ NEL = IEL
+ NTMP(IEL) = I-1
+ ENDIF
+C
+C
+ WRITE(*,3010) NEL
+ CALL ASKI('Enter element to be replaced by current airfoil^',IEL)
+C
+ IF(IEL.LT.1 .OR. IEL.GT.NEL+1) THEN
+ WRITE(*,*) 'Element number inappropriate. Airfoil not written.'
+ CLOSE(LU)
+ RETURN
+ ELSE IF(IEL.EQ.NEL+1) THEN
+ NEL = NEL+1
+ ENDIF
+C
+C
+ NTMP(IEL) = N
+ DO 70 I = 1, NTMP(IEL)
+ IF(LCLOCK) THEN
+C------- write out in clockwise order (reversed from internal XFOIL order)
+ IDIR = NTMP(IEL) - I + 1
+ ELSE
+C------- write out in counterclockwise order (same as internal XFOIL order)
+ IDIR = I
+ ENDIF
+ XTMP(I,IEL) = X(IDIR)
+ YTMP(I,IEL) = Y(IDIR)
+ 70 CONTINUE
+C
+C
+ REWIND(LU)
+C
+C---- write first 2 lines of MSES format coordinate file
+ WRITE(LU,1000) NAME1(1:NN1)
+ WRITE(LU,1000) ISPARS1(1:NI1)
+C
+ DO 80 IEL=1, NEL
+ DO 805 I=1, NTMP(IEL)
+ WRITE(LU,1100) XTMP(I,IEL),YTMP(I,IEL)
+ 805 CONTINUE
+ IF(IEL.LT.NEL) WRITE(LU,*) ' 999.0 999.0'
+ 80 CONTINUE
+C
+ CLOSE(LU)
+ RETURN
+C
+ 9005 WRITE(*,*) 'Old file OPEN error. Airfoil not saved.'
+ RETURN
+C
+ 9010 WRITE(*,*) 'Old file READ error. Airfoil not saved.'
+ CLOSE(LU)
+ RETURN
+C
+ 1000 FORMAT(A)
+ 1100 FORMAT(1X,2G15.7)
+ 3010 FORMAT(/' Specified multielement airfoil has',I2,' elements.')
+ END ! MSAVE
+
+
+
+ SUBROUTINE ROTATE(X,Y,N,ALFA)
+ DIMENSION X(N), Y(N)
+C
+ SA = SIN(ALFA)
+ CA = COS(ALFA)
+CCC XOFF = 0.25*(1.0-CA)
+CCC YOFF = 0.25*SA
+ XOFF = 0.
+ YOFF = 0.
+ DO 8 I=1, N
+ XT = X(I)
+ YT = Y(I)
+ X(I) = CA*XT + SA*YT + XOFF
+ Y(I) = CA*YT - SA*XT + YOFF
+ 8 CONTINUE
+C
+ RETURN
+ END
+
+
+ SUBROUTINE NACA(IDES1)
+ INCLUDE 'XFOIL.INC'
+C
+C---- number of points per side
+ NSIDE = IQX/3
+C
+ IF(IDES1 .LE. 0) THEN
+ CALL ASKI('Enter NACA 4 or 5-digit airfoil designation^',IDES)
+ ELSE
+ IDES = IDES1
+ ENDIF
+C
+ ITYPE = 0
+ IF(IDES.LE.25099) ITYPE = 5
+ IF(IDES.LE.9999 ) ITYPE = 4
+C
+ IF(ITYPE.EQ.0) THEN
+ WRITE(*,*) 'This designation not implemented.'
+ RETURN
+ ENDIF
+C
+ IF(ITYPE.EQ.4) CALL NACA4(IDES,W1,W2,W3,NSIDE,XB,YB,NB,NAME)
+ IF(ITYPE.EQ.5) CALL NACA5(IDES,W1,W2,W3,NSIDE,XB,YB,NB,NAME)
+ CALL STRIP(NAME,NNAME)
+C
+C---- see if routines didn't recognize designator
+ IF(IDES.EQ.0) RETURN
+C
+ LCLOCK = .FALSE.
+C
+ XBF = 0.0
+ YBF = 0.0
+ LBFLAP = .FALSE.
+C
+ CALL SCALC(XB,YB,SB,NB)
+ CALL SEGSPL(XB,XBP,SB,NB)
+ CALL SEGSPL(YB,YBP,SB,NB)
+C
+ CALL GEOPAR(XB,XBP,YB,YBP,SB,NB, W1,
+ & SBLE,CHORDB,AREAB,RADBLE,ANGBTE,
+ & EI11BA,EI22BA,APX1BA,APX2BA,
+ & EI11BT,EI22BT,APX1BT,APX2BT,
+ & THICKB,CAMBRB )
+C
+ WRITE(*,1200) NB
+ 1200 FORMAT(/' Buffer airfoil set using', I4,' points')
+C
+C---- set paneling
+ CALL PANGEN(.TRUE.)
+ccc CALL PANPLT
+C
+ RETURN
+ END ! NACA
+
+
+ SUBROUTINE PANGEN(SHOPAR)
+C---------------------------------------------------
+C Set paneling distribution from buffer airfoil
+C geometry, thus creating current airfoil.
+C
+C If REFINE=True, bunch points at x=XSREF on
+C top side and at x=XPREF on bottom side
+C by setting a fictitious local curvature of
+C CTRRAT*(LE curvature) there.
+C---------------------------------------------------
+ INCLUDE 'XFOIL.INC'
+ LOGICAL SHOPAR
+C
+ IF(NB.LT.2) THEN
+ WRITE(*,*) 'PANGEN: Buffer airfoil not available.'
+ N = 0
+ RETURN
+ ENDIF
+C
+C---- Number of temporary nodes for panel distribution calculation
+C exceeds the specified panel number by factor of IPFAC.
+ IPFAC = 3
+ IPFAC = 5
+C
+C---- number of airfoil panel points
+ N = NPAN
+C
+cC---- number of wake points
+c NW = NPAN/8 + 2
+c IF(NW.GT.IWX) THEN
+c WRITE(*,*)
+c & 'Array size (IWX) too small. Last wake point index reduced.'
+c NW = IWX
+c ENDIF
+C
+C---- set arc length spline parameter
+ CALL SCALC(XB,YB,SB,NB)
+C
+C---- spline raw airfoil coordinates
+ CALL SEGSPL(XB,XBP,SB,NB)
+ CALL SEGSPL(YB,YBP,SB,NB)
+C
+C---- normalizing length (~ chord)
+ SBREF = 0.5*(SB(NB)-SB(1))
+C
+C---- set up curvature array
+ DO I = 1, NB
+ W5(I) = ABS( CURV(SB(I),XB,XBP,YB,YBP,SB,NB) ) * SBREF
+ ENDDO
+C
+C---- locate LE point arc length value and the normalized curvature there
+ CALL LEFIND(SBLE,XB,XBP,YB,YBP,SB,NB)
+ CVLE = ABS( CURV(SBLE,XB,XBP,YB,YBP,SB,NB) ) * SBREF
+C
+C---- check for doubled point (sharp corner) at LE
+ IBLE = 0
+ DO I = 1, NB-1
+ IF(SBLE.EQ.SB(I) .AND. SBLE.EQ.SB(I+1)) THEN
+ IBLE = I
+ WRITE(*,*)
+ WRITE(*,*) 'Sharp leading edge'
+ GO TO 21
+ ENDIF
+ ENDDO
+ 21 CONTINUE
+C
+C---- set LE, TE points
+ XBLE = SEVAL(SBLE,XB,XBP,SB,NB)
+ YBLE = SEVAL(SBLE,YB,YBP,SB,NB)
+ XBTE = 0.5*(XB(1)+XB(NB))
+ YBTE = 0.5*(YB(1)+YB(NB))
+ CHBSQ = (XBTE-XBLE)**2 + (YBTE-YBLE)**2
+C
+C---- set average curvature over 2*NK+1 points within Rcurv of LE point
+ NK = 3
+ CVSUM = 0.
+ DO K = -NK, NK
+ FRAC = FLOAT(K)/FLOAT(NK)
+ SBK = SBLE + FRAC*SBREF/MAX(CVLE,20.0)
+ CVK = ABS( CURV(SBK,XB,XBP,YB,YBP,SB,NB) ) * SBREF
+ CVSUM = CVSUM + CVK
+ ENDDO
+ CVAVG = CVSUM/FLOAT(2*NK+1)
+C
+C---- dummy curvature for sharp LE
+ IF(IBLE.NE.0) CVAVG = 10.0
+C
+C---- set curvature attraction coefficient actually used
+ CC = 6.0 * CVPAR
+C
+C---- set artificial curvature at TE to bunch panels there
+ CVTE = CVAVG * CTERAT
+ W5(1) = CVTE
+ W5(NB) = CVTE
+C
+C
+C**** smooth curvature array for smoother panel size distribution ****
+C
+CCC CALL ASKR('Enter curvature smoothing length/c^',SMOOL)
+CCC SMOOL = 0.010
+C
+C---- set smoothing length = 1 / averaged LE curvature, but
+C- no more than 5% of chord and no less than 1/4 average panel spacing
+ SMOOL = MAX( 1.0/MAX(CVAVG,20.0) , 0.25 /FLOAT(NPAN/2) )
+C
+ SMOOSQ = (SMOOL*SBREF) ** 2
+C
+C---- set up tri-diagonal system for smoothed curvatures
+ W2(1) = 1.0
+ W3(1) = 0.0
+ DO I=2, NB-1
+ DSM = SB(I) - SB(I-1)
+ DSP = SB(I+1) - SB(I)
+ DSO = 0.5*(SB(I+1) - SB(I-1))
+C
+ IF(DSM.EQ.0.0 .OR. DSP.EQ.0.0) THEN
+C------- leave curvature at corner point unchanged
+ W1(I) = 0.0
+ W2(I) = 1.0
+ W3(I) = 0.0
+ ELSE
+ W1(I) = SMOOSQ * ( - 1.0/DSM) / DSO
+ W2(I) = SMOOSQ * ( 1.0/DSP + 1.0/DSM) / DSO + 1.0
+ W3(I) = SMOOSQ * (-1.0/DSP ) / DSO
+ ENDIF
+ ENDDO
+C
+ W1(NB) = 0.0
+ W2(NB) = 1.0
+C
+C---- fix curvature at LE point by modifying equations adjacent to LE
+ DO I=2, NB-1
+ IF(SB(I).EQ.SBLE .OR. I.EQ.IBLE .OR. I.EQ.IBLE+1) THEN
+C------- if node falls right on LE point, fix curvature there
+ W1(I) = 0.
+ W2(I) = 1.0
+ W3(I) = 0.
+ W5(I) = CVLE
+ ELSE IF(SB(I-1).LT.SBLE .AND. SB(I).GT.SBLE) THEN
+C------- modify equation at node just before LE point
+ DSM = SB(I-1) - SB(I-2)
+ DSP = SBLE - SB(I-1)
+ DSO = 0.5*(SBLE - SB(I-2))
+C
+ W1(I-1) = SMOOSQ * ( - 1.0/DSM) / DSO
+ W2(I-1) = SMOOSQ * ( 1.0/DSP + 1.0/DSM) / DSO + 1.0
+ W3(I-1) = 0.
+ W5(I-1) = W5(I-1) + SMOOSQ*CVLE/(DSP*DSO)
+C
+C------- modify equation at node just after LE point
+ DSM = SB(I) - SBLE
+ DSP = SB(I+1) - SB(I)
+ DSO = 0.5*(SB(I+1) - SBLE)
+ W1(I) = 0.
+ W2(I) = SMOOSQ * ( 1.0/DSP + 1.0/DSM) / DSO + 1.0
+ W3(I) = SMOOSQ * (-1.0/DSP ) / DSO
+ W5(I) = W5(I) + SMOOSQ*CVLE/(DSM*DSO)
+C
+ GO TO 51
+ ENDIF
+ ENDDO
+ 51 CONTINUE
+C
+C---- set artificial curvature at bunching points and fix it there
+ DO I=2, NB-1
+C------ chord-based x/c coordinate
+ XOC = ( (XB(I)-XBLE)*(XBTE-XBLE)
+ & + (YB(I)-YBLE)*(YBTE-YBLE) ) / CHBSQ
+C
+ IF(SB(I).LT.SBLE) THEN
+C------- check if top side point is in refinement area
+ IF(XOC.GT.XSREF1 .AND. XOC.LT.XSREF2) THEN
+ W1(I) = 0.
+ W2(I) = 1.0
+ W3(I) = 0.
+ W5(I) = CVLE*CTRRAT
+ ENDIF
+ ELSE
+C------- check if bottom side point is in refinement area
+ IF(XOC.GT.XPREF1 .AND. XOC.LT.XPREF2) THEN
+ W1(I) = 0.
+ W2(I) = 1.0
+ W3(I) = 0.
+ W5(I) = CVLE*CTRRAT
+ ENDIF
+ ENDIF
+ ENDDO
+C
+C---- solve for smoothed curvature array W5
+ IF(IBLE.EQ.0) THEN
+ CALL TRISOL(W2,W1,W3,W5,NB)
+ ELSE
+ I = 1
+ CALL TRISOL(W2(I),W1(I),W3(I),W5(I),IBLE)
+ I = IBLE+1
+ CALL TRISOL(W2(I),W1(I),W3(I),W5(I),NB-IBLE)
+ ENDIF
+C
+C---- find max curvature
+ CVMAX = 0.
+ DO I=1, NB
+ CVMAX = MAX( CVMAX , ABS(W5(I)) )
+ ENDDO
+C
+C---- normalize curvature array
+ DO I=1, NB
+ W5(I) = W5(I) / CVMAX
+ ENDDO
+C
+C---- spline curvature array
+ CALL SEGSPL(W5,W6,SB,NB)
+C
+C---- Set initial guess for node positions uniform in s.
+C More nodes than specified (by factor of IPFAC) are
+C temporarily used for more reliable convergence.
+ NN = IPFAC*(N-1)+1
+C
+C---- ratio of lengths of panel at TE to one away from the TE
+ RDSTE = 0.667
+ RTF = (RDSTE-1.0)*2.0 + 1.0
+C
+ IF(IBLE.EQ.0) THEN
+C
+ DSAVG = (SB(NB)-SB(1))/(FLOAT(NN-3) + 2.0*RTF)
+ SNEW(1) = SB(1)
+ DO I=2, NN-1
+ SNEW(I) = SB(1) + DSAVG * (FLOAT(I-2) + RTF)
+ ENDDO
+ SNEW(NN) = SB(NB)
+C
+ ELSE
+C
+ NFRAC1 = (N * IBLE) / NB
+C
+ NN1 = IPFAC*(NFRAC1-1)+1
+ DSAVG1 = (SBLE-SB(1))/(FLOAT(NN1-2) + RTF)
+ SNEW(1) = SB(1)
+ DO I=2, NN1
+ SNEW(I) = SB(1) + DSAVG1 * (FLOAT(I-2) + RTF)
+ ENDDO
+C
+ NN2 = NN - NN1 + 1
+ DSAVG2 = (SB(NB)-SBLE)/(FLOAT(NN2-2) + RTF)
+ DO I=2, NN2-1
+ SNEW(I-1+NN1) = SBLE + DSAVG2 * (FLOAT(I-2) + RTF)
+ ENDDO
+ SNEW(NN) = SB(NB)
+C
+ ENDIF
+C
+C---- Newton iteration loop for new node positions
+ DO 10 ITER=1, 20
+C
+C------ set up tri-diagonal system for node position deltas
+ CV1 = SEVAL(SNEW(1),W5,W6,SB,NB)
+ CV2 = SEVAL(SNEW(2),W5,W6,SB,NB)
+ CVS1 = DEVAL(SNEW(1),W5,W6,SB,NB)
+ CVS2 = DEVAL(SNEW(2),W5,W6,SB,NB)
+C
+ CAVM = SQRT(CV1**2 + CV2**2)
+ IF(CAVM .EQ. 0.0) THEN
+ CAVM_S1 = 0.
+ CAVM_S2 = 0.
+ ELSE
+ CAVM_S1 = CVS1 * CV1/CAVM
+ CAVM_S2 = CVS2 * CV2/CAVM
+ ENDIF
+C
+ DO 110 I=2, NN-1
+ DSM = SNEW(I) - SNEW(I-1)
+ DSP = SNEW(I) - SNEW(I+1)
+ CV3 = SEVAL(SNEW(I+1),W5,W6,SB,NB)
+ CVS3 = DEVAL(SNEW(I+1),W5,W6,SB,NB)
+C
+ CAVP = SQRT(CV3**2 + CV2**2)
+ IF(CAVP .EQ. 0.0) THEN
+ CAVP_S2 = 0.
+ CAVP_S3 = 0.
+ ELSE
+ CAVP_S2 = CVS2 * CV2/CAVP
+ CAVP_S3 = CVS3 * CV3/CAVP
+ ENDIF
+C
+ FM = CC*CAVM + 1.0
+ FP = CC*CAVP + 1.0
+C
+ REZ = DSP*FP + DSM*FM
+C
+C-------- lower, main, and upper diagonals
+ W1(I) = -FM + CC* DSM*CAVM_S1
+ W2(I) = FP + FM + CC*(DSP*CAVP_S2 + DSM*CAVM_S2)
+ W3(I) = -FP + CC* DSP*CAVP_S3
+C
+C-------- residual, requiring that
+C (1 + C*curv)*deltaS is equal on both sides of node i
+ W4(I) = -REZ
+C
+ CV1 = CV2
+ CV2 = CV3
+ CVS1 = CVS2
+ CVS2 = CVS3
+ CAVM = CAVP
+ CAVM_S1 = CAVP_S2
+ CAVM_S2 = CAVP_S3
+ 110 CONTINUE
+C
+C------ fix endpoints (at TE)
+ W2(1) = 1.0
+ W3(1) = 0.0
+ W4(1) = 0.0
+ W1(NN) = 0.0
+ W2(NN) = 1.0
+ W4(NN) = 0.0
+C
+ IF(RTF .NE. 1.0) THEN
+C------- fudge equations adjacent to TE to get TE panel length ratio RTF
+C
+ I = 2
+ W4(I) = -((SNEW(I) - SNEW(I-1)) + RTF*(SNEW(I) - SNEW(I+1)))
+ W1(I) = -1.0
+ W2(I) = 1.0 + RTF
+ W3(I) = - RTF
+C
+ I = NN-1
+ W4(I) = -((SNEW(I) - SNEW(I+1)) + RTF*(SNEW(I) - SNEW(I-1)))
+ W3(I) = -1.0
+ W2(I) = 1.0 + RTF
+ W1(I) = - RTF
+ ENDIF
+C
+C
+C------ fix sharp LE point
+ IF(IBLE.NE.0) THEN
+ I = NN1
+ W1(I) = 0.0
+ W2(I) = 1.0
+ W3(I) = 0.0
+ W4(I) = SBLE - SNEW(I)
+ ENDIF
+C
+C------ solve for changes W4 in node position arc length values
+ CALL TRISOL(W2,W1,W3,W4,NN)
+C
+C------ find under-relaxation factor to keep nodes from changing order
+ RLX = 1.0
+ DMAX = 0.0
+ DO I=1, NN-1
+ DS = SNEW(I+1) - SNEW(I)
+ DDS = W4(I+1) - W4(I)
+ DSRAT = 1.0 + RLX*DDS/DS
+ IF(DSRAT.GT.4.0) RLX = (4.0-1.0)*DS/DDS
+ IF(DSRAT.LT.0.2) RLX = (0.2-1.0)*DS/DDS
+ DMAX = MAX(ABS(W4(I)),DMAX)
+ ENDDO
+C
+C------ update node position
+ DO I=2, NN-1
+ SNEW(I) = SNEW(I) + RLX*W4(I)
+ ENDDO
+C
+CCC IF(RLX.EQ.1.0) WRITE(*,*) DMAX
+CCC IF(RLX.NE.1.0) WRITE(*,*) DMAX,' RLX =',RLX
+ IF(ABS(DMAX).LT.1.E-3) GO TO 11
+ 10 CONTINUE
+ WRITE(*,*) 'Paneling convergence failed. Continuing anyway...'
+C
+ 11 CONTINUE
+C
+C---- set new panel node coordinates
+ DO I=1, N
+ IND = IPFAC*(I-1) + 1
+ S(I) = SNEW(IND)
+ X(I) = SEVAL(SNEW(IND),XB,XBP,SB,NB)
+ Y(I) = SEVAL(SNEW(IND),YB,YBP,SB,NB)
+ ENDDO
+C
+C
+C---- go over buffer airfoil again, checking for corners (double points)
+ NCORN = 0
+ DO 25 IB=1, NB-1
+ IF(SB(IB) .EQ. SB(IB+1)) THEN
+C------- found one !
+C
+ NCORN = NCORN+1
+ XBCORN = XB(IB)
+ YBCORN = YB(IB)
+ SBCORN = SB(IB)
+C
+C------- find current-airfoil panel which contains corner
+ DO 252 I=1, N
+C
+C--------- keep stepping until first node past corner
+ IF(S(I) .LE. SBCORN) GO TO 252
+C
+C---------- move remainder of panel nodes to make room for additional node
+ DO 2522 J=N, I, -1
+ X(J+1) = X(J)
+ Y(J+1) = Y(J)
+ S(J+1) = S(J)
+ 2522 CONTINUE
+ N = N+1
+C
+ IF(N .GT. IQX-1)
+ & STOP 'PANEL: Too many panels. Increase IQX in XFOIL.INC'
+C
+ X(I) = XBCORN
+ Y(I) = YBCORN
+ S(I) = SBCORN
+C
+C---------- shift nodes adjacent to corner to keep panel sizes comparable
+ IF(I-2 .GE. 1) THEN
+ S(I-1) = 0.5*(S(I) + S(I-2))
+ X(I-1) = SEVAL(S(I-1),XB,XBP,SB,NB)
+ Y(I-1) = SEVAL(S(I-1),YB,YBP,SB,NB)
+ ENDIF
+C
+ IF(I+2 .LE. N) THEN
+ S(I+1) = 0.5*(S(I) + S(I+2))
+ X(I+1) = SEVAL(S(I+1),XB,XBP,SB,NB)
+ Y(I+1) = SEVAL(S(I+1),YB,YBP,SB,NB)
+ ENDIF
+C
+C---------- go on to next input geometry point to check for corner
+ GO TO 25
+C
+ 252 CONTINUE
+ ENDIF
+ 25 CONTINUE
+C
+ CALL SCALC(X,Y,S,N)
+ CALL SEGSPL(X,XP,S,N)
+ CALL SEGSPL(Y,YP,S,N)
+ CALL LEFIND(SLE,X,XP,Y,YP,S,N)
+C
+ XLE = SEVAL(SLE,X,XP,S,N)
+ YLE = SEVAL(SLE,Y,YP,S,N)
+ XTE = 0.5*(X(1)+X(N))
+ YTE = 0.5*(Y(1)+Y(N))
+ CHORD = SQRT( (XTE-XLE)**2 + (YTE-YLE)**2 )
+C
+C---- calculate panel size ratios (user info)
+ DSMIN = 1000.0
+ DSMAX = -1000.0
+ DO 40 I=1, N-1
+ DS = S(I+1)-S(I)
+ IF(DS .EQ. 0.0) GO TO 40
+ DSMIN = MIN(DSMIN,DS)
+ DSMAX = MAX(DSMAX,DS)
+ 40 CONTINUE
+C
+ DSMIN = DSMIN*FLOAT(N-1)/S(N)
+ DSMAX = DSMAX*FLOAT(N-1)/S(N)
+ccc WRITE(*,*) 'DSmin/DSavg = ',DSMIN,' DSmax/DSavg = ',DSMAX
+C
+C---- set various flags for new airfoil
+ LGAMU = .FALSE.
+ LQINU = .FALSE.
+ LWAKE = .FALSE.
+ LQAIJ = .FALSE.
+ LADIJ = .FALSE.
+ LWDIJ = .FALSE.
+ LIPAN = .FALSE.
+ LBLINI = .FALSE.
+ LVCONV = .FALSE.
+ LSCINI = .FALSE.
+ LQSPEC = .FALSE.
+ LGSAME = .FALSE.
+C
+ IF(LBFLAP) THEN
+ XOF = XBF
+ YOF = YBF
+ LFLAP = .TRUE.
+ ENDIF
+C
+C---- determine if TE is blunt or sharp, calculate TE geometry parameters
+ CALL TECALC
+C
+C---- calculate normal vectors
+ CALL NCALC(X,Y,S,N,NX,NY)
+C
+C---- calculate panel angles for panel routines
+ CALL APCALC
+C
+ IF(SHARP) THEN
+ WRITE(*,1090) 'Sharp trailing edge'
+ ELSE
+ GAP = SQRT((X(1)-X(N))**2 + (Y(1)-Y(N))**2)
+ WRITE(*,1090) 'Blunt trailing edge. Gap =', GAP
+ ENDIF
+ 1090 FORMAT(/1X,A,F9.5)
+C
+ IF(SHOPAR) WRITE(*,1100) NPAN, CVPAR, CTERAT, CTRRAT,
+ & XSREF1, XSREF2, XPREF1, XPREF2
+ 1100 FORMAT(/' Paneling parameters used...'
+ & /' Number of panel nodes ' , I4
+ & /' Panel bunching parameter ' , F6.3
+ & /' TE/LE panel density ratio ' , F6.3
+ & /' Refined-area/LE panel density ratio ' , F6.3
+ & /' Top side refined area x/c limits ' , 2F6.3
+ & /' Bottom side refined area x/c limits ' , 2F6.3)
+C
+ RETURN
+ END ! PANGEN
+
+
+
+ SUBROUTINE GETPAN
+ INCLUDE 'XFOIL.INC'
+ LOGICAL LCHANGE
+ CHARACTER*4 VAR
+ CHARACTER*128 COMARG
+C
+ DIMENSION IINPUT(20)
+ DIMENSION RINPUT(20)
+ LOGICAL ERROR
+C
+ IF(NB.LE.1) THEN
+ WRITE(*,*) 'GETPAN: Buffer airfoil not available.'
+ RETURN
+ ENDIF
+C
+ 5 CONTINUE
+ IF(N.LE.1) THEN
+ WRITE(*,*) 'No current airfoil to plot'
+ ELSE
+ CALL PANPLT
+ ENDIF
+ LCHANGE = .FALSE.
+C
+ 10 WRITE(*,1000) NPAN, CVPAR, CTERAT, CTRRAT,
+ & XSREF1, XSREF2, XPREF1, XPREF2
+ 1000 FORMAT(
+ & /' Present paneling parameters...'
+ & /' N i Number of panel nodes ' , I4
+ & /' P r Panel bunching parameter ' , F6.3
+ & /' T r TE/LE panel density ratio ' , F6.3
+ & /' R r Refined area/LE panel density ratio ' , F6.3
+ & /' XT rr Top side refined area x/c limits ' , 2F6.3
+ & /' XB rr Bottom side refined area x/c limits ' , 2F6.3
+ & /' Z oom'
+ & /' U nzoom' )
+C
+ 12 CALL ASKC('Change what ? (<cr> if nothing else)^',VAR,COMARG)
+C
+ IF(VAR.EQ.'Z ') THEN
+ CALL USETZOOM(.TRUE.,.TRUE.)
+ CALL REPLOT(IDEV)
+ GO TO 12
+ ENDIF
+C
+ IF(VAR.EQ.'U ') THEN
+ CALL CLRZOOM
+ CALL REPLOT(IDEV)
+ GO TO 12
+ ENDIF
+C
+C
+ DO I=1, 20
+ IINPUT(I) = 0
+ RINPUT(I) = 0.0
+ ENDDO
+ NINPUT = 0
+ CALL GETINT(COMARG,IINPUT,NINPUT,ERROR)
+ NINPUT = 0
+ CALL GETFLT(COMARG,RINPUT,NINPUT,ERROR)
+C
+ IF (VAR.EQ.' ') THEN
+C
+ IF(LCHANGE) THEN
+C
+C-------- set new panel distribution, and display max panel corner angle
+ CALL PANGEN(.FALSE.)
+ IF(N.GT.0) CALL CANG(X,Y,N,1,IMAX,AMAX)
+C
+C-------- go back to paneling menu
+ GO TO 5
+ ENDIF
+C
+ CALL CLRZOOM
+ RETURN
+C
+ ELSE IF(VAR.EQ.'N ' .OR. VAR.EQ.'n ') THEN
+C
+ IF(NINPUT.GE.1) THEN
+ NPAN = IINPUT(1)
+ ELSE
+ CALL ASKI('Enter number of panel nodes^',NPAN)
+ ENDIF
+ IF(NPAN .GT. IQX-6) THEN
+ NPAN = IQX - 6
+ WRITE(*,1200) NPAN
+ 1200 FORMAT(1X,' Number of panel nodes reduced to array limit:',I4)
+ ENDIF
+ LCHANGE = .TRUE.
+C
+ ELSE IF(VAR.EQ.'P ' .OR. VAR.EQ.'p ') THEN
+C
+ IF(NINPUT.GE.1) THEN
+ CVPAR = RINPUT(1)
+ ELSE
+ CALL ASKR('Enter panel bunching parameter (0 to ~1)^',CVPAR)
+ ENDIF
+ LCHANGE = .TRUE.
+C
+ ELSE IF(VAR.EQ.'T ' .OR. VAR.EQ.'t ') THEN
+C
+ IF(NINPUT.GE.1) THEN
+ CTERAT = RINPUT(1)
+ ELSE
+ CALL ASKR('Enter TE/LE panel density ratio^',CTERAT)
+ ENDIF
+ LCHANGE = .TRUE.
+C
+ ELSE IF(VAR.EQ.'R ' .OR. VAR.EQ.'r ') THEN
+C
+ IF(NINPUT.GE.1) THEN
+ CTRRAT = RINPUT(1)
+ ELSE
+ CALL ASKR('Enter refined-area panel density ratio^',CTRRAT)
+ ENDIF
+ LCHANGE = .TRUE.
+C
+ ELSE IF(VAR.EQ.'XT ' .OR. VAR.EQ.'xt ') THEN
+C
+ IF(NINPUT.GE.2) THEN
+ XSREF1 = RINPUT(1)
+ XSREF2 = RINPUT(2)
+ ELSE
+ CALL ASKR('Enter left top side refinement limit^',XSREF1)
+ CALL ASKR('Enter right top side refinement limit^',XSREF2)
+ ENDIF
+ LCHANGE = .TRUE.
+C
+ ELSE IF(VAR.EQ.'XB ' .OR. VAR.EQ.'xb ') THEN
+C
+ IF(NINPUT.GE.2) THEN
+ XPREF1 = RINPUT(1)
+ XPREF2 = RINPUT(2)
+ ELSE
+ CALL ASKR('Enter left bottom side refinement limit^',XPREF1)
+ CALL ASKR('Enter right bottom side refinement limit^',XPREF2)
+ ENDIF
+ LCHANGE = .TRUE.
+C
+ ELSE
+C
+ WRITE(*,*)
+ WRITE(*,*) '*** Input not recognized ***'
+ GO TO 10
+C
+ ENDIF
+C
+ GO TO 12
+C
+ END ! GETPAN
+
+
+ SUBROUTINE TECALC
+C-------------------------------------------
+C Calculates total and projected TE gap
+C areas and TE panel strengths.
+C-------------------------------------------
+ INCLUDE 'XFOIL.INC'
+C
+C---- set TE base vector and TE bisector components
+ DXTE = X(1) - X(N)
+ DYTE = Y(1) - Y(N)
+ DXS = 0.5*(-XP(1) + XP(N))
+ DYS = 0.5*(-YP(1) + YP(N))
+C
+C---- normal and streamwise projected TE gap areas
+ ANTE = DXS*DYTE - DYS*DXTE
+ ASTE = DXS*DXTE + DYS*DYTE
+C
+C---- total TE gap area
+ DSTE = SQRT(DXTE**2 + DYTE**2)
+C
+ SHARP = DSTE .LT. 0.0001*CHORD
+C
+ IF(SHARP) THEN
+ SCS = 1.0
+ SDS = 0.0
+ ELSE
+ SCS = ANTE/DSTE
+ SDS = ASTE/DSTE
+ ENDIF
+C
+C---- TE panel source and vorticity strengths
+ SIGTE = 0.5*(GAM(1) - GAM(N))*SCS
+ GAMTE = -.5*(GAM(1) - GAM(N))*SDS
+C
+ SIGTE_A = 0.5*(GAM_A(1) - GAM_A(N))*SCS
+ GAMTE_A = -.5*(GAM_A(1) - GAM_A(N))*SDS
+C
+ RETURN
+ END ! TECALC
+
+
+
+ SUBROUTINE INTE
+C-----------------------------------------------------------
+C Interpolates two airfoils into an intermediate shape.
+C Extrapolation is also possible to a reasonable extent.
+C-----------------------------------------------------------
+ INCLUDE 'XFOIL.INC'
+ CHARACTER*2 CAIR
+ INTEGER NINT(2)
+ REAL SINT(IBX,2),
+ & XINT(IBX,2), XPINT(IBX,2),
+ & YINT(IBX,2), YPINT(IBX,2),
+ & SLEINT(2)
+ CHARACTER*20 PROMPTN
+ CHARACTER*48 NAMEINT(2)
+ CHARACTER*80 ISPARST
+C
+ LU = 21
+C
+ 1000 FORMAT(A)
+C
+ WRITE(*,1100) NAME
+ DO IP=1, NPOL
+ IF(NXYPOL(IP).GT.0) THEN
+ WRITE(*,1200) IP, NAMEPOL(IP)
+ ENDIF
+ ENDDO
+ IF (NPOL.EQ.0) THEN
+ PROMPTN = '" ( F C ): '
+ NPR = 12
+ ELSEIF(NPOL.EQ.1) THEN
+ PROMPTN = '" ( F C 1 ): '
+ NPR = 14
+ ELSEIF(NPOL.EQ.2) THEN
+ PROMPTN = '" ( F C 1 2 ): '
+ NPR = 16
+ ELSE
+ PROMPTN = '" ( F C 1 2.. ): '
+ NPR = 18
+ ENDIF
+C
+ 1100 FORMAT(/ ' F disk file'
+ & / ' C current airfoil ', A)
+ 1200 FORMAT( 1X,I2,' polar airfoil ', A)
+C
+ 2100 FORMAT(/' Select source of airfoil "',I1, A, $)
+C
+ DO 40 K = 1, 2
+ IAIR = K - 1
+ 20 WRITE(*,2100) IAIR, PROMPTN(1:NPR)
+ READ(*,1000) CAIR
+C
+ IF (INDEX('Ff',CAIR(1:1)).NE.0) THEN
+ CALL ASKS('Enter filename^',FNAME)
+ CALL AREAD(LU,FNAME,IBX,
+ & XINT(1,K),YINT(1,K),NINT(K),
+ & NAMEINT(K),ISPARST,ITYPE,0)
+ IF(ITYPE.EQ.0) RETURN
+C
+ ELSEIF(INDEX('Cc',CAIR(1:1)).NE.0) THEN
+ IF(N.LE.1) THEN
+ WRITE(*,*) 'No current airfoil available'
+ GO TO 20
+ ENDIF
+C
+ NINT(K) = N
+ DO I = 1, N
+ XINT(I,K) = X(I)
+ YINT(I,K) = Y(I)
+ ENDDO
+ NAMEINT(K) = NAME
+C
+ ELSE
+ READ(CAIR,*,ERR=90) IP
+ IF(IP.LT.1 .OR. IP.GT.NPOL) THEN
+ GO TO 90
+ ELSEIF(NXYPOL(IP).LE.0) THEN
+ GO TO 90
+ ELSE
+ NINT(K) = NXYPOL(IP)
+ DO I = 1, N
+ XINT(I,K) = CPOLXY(I,1,IP)
+ YINT(I,K) = CPOLXY(I,2,IP)
+ ENDDO
+ ENDIF
+ NAMEINT(K) = NAMEPOL(IP)
+C
+ ENDIF
+C
+ CALL SCALC(XINT(1,K),YINT(1,K),SINT(1,K),NINT(K))
+ CALL SEGSPLD(XINT(1,K),XPINT(1,K),SINT(1,K),NINT(K),-999.,-999.)
+ CALL SEGSPLD(YINT(1,K),YPINT(1,K),SINT(1,K),NINT(K),-999.,-999.)
+ CALL LEFIND(SLEINT(K),
+ & XINT(1,K),XPINT(1,K),
+ & YINT(1,K),YPINT(1,K),SINT(1,K),NINT(K))
+ 40 CONTINUE
+C
+ WRITE(*,*)
+ WRITE(*,*) 'airfoil "0": ', NAMEINT(1)
+ WRITE(*,*) 'airfoil "1": ', NAMEINT(2)
+ FRAC = 0.5
+ CALL ASKR('Specify interpolating fraction 0...1^',FRAC)
+C
+ CALL INTER(XINT(1,1),XPINT(1,1),
+ & YINT(1,1),YPINT(1,1),SINT(1,1),NINT(1),SLEINT(1),
+ & XINT(1,2),XPINT(1,2),
+ & YINT(1,2),YPINT(1,2),SINT(1,2),NINT(2),SLEINT(2),
+ & XB,YB,NB,FRAC)
+C
+ CALL SCALC(XB,YB,SB,NB)
+ CALL SEGSPL(XB,XBP,SB,NB)
+ CALL SEGSPL(YB,YBP,SB,NB)
+C
+ CALL GEOPAR(XB,XBP,YB,YBP,SB,NB, W1,
+ & SBLE,CHORDB,AREAB,RADBLE,ANGBTE,
+ & EI11BA,EI22BA,APX1BA,APX2BA,
+ & EI11BT,EI22BT,APX1BT,APX2BT,
+ & THICKB,CAMBRB )
+C
+ CALL ASKS('Enter new airfoil name^',NAME)
+ CALL STRIP(NAME,NNAME)
+ WRITE(*,*)
+ WRITE(*,*) 'Result has been placed in buffer airfoil'
+ WRITE(*,*) 'Execute PCOP or PANE to set new current airfoil'
+ RETURN
+C
+ 90 CONTINUE
+ WRITE(*,*)
+ WRITE(*,*) 'Invalid response'
+ RETURN
+ END ! INTE
+
+
+ SUBROUTINE INTX
+C-----------------------------------------------------------
+C Interpolates two airfoils into an intermediate shape.
+C Extrapolation is also possible to a reasonable extent.
+C-----------------------------------------------------------
+ INCLUDE 'XFOIL.INC'
+ CHARACTER*2 CAIR
+ INTEGER NINT(2)
+ REAL SINT(IBX,2),
+ & XINT(IBX,2), XPINT(IBX,2),
+ & YINT(IBX,2), YPINT(IBX,2),
+ & SLEINT(2)
+ CHARACTER*20 PROMPTN
+ CHARACTER*48 NAMEINT(2)
+ CHARACTER*80 ISPARST
+C
+ LU = 21
+C
+ 1000 FORMAT(A)
+C
+ WRITE(*,1100) NAME
+ DO IP=1, NPOL
+ IF(NXYPOL(IP).GT.0) THEN
+ WRITE(*,1200) IP, NAMEPOL(IP)
+ ENDIF
+ ENDDO
+ IF (NPOL.EQ.0) THEN
+ PROMPTN = '" ( F C ): '
+ NPR = 12
+ ELSEIF(NPOL.EQ.1) THEN
+ PROMPTN = '" ( F C 1 ): '
+ NPR = 14
+ ELSEIF(NPOL.EQ.2) THEN
+ PROMPTN = '" ( F C 1 2 ): '
+ NPR = 16
+ ELSE
+ PROMPTN = '" ( F C 1 2.. ): '
+ NPR = 18
+ ENDIF
+C
+ 1100 FORMAT(/ ' F disk file'
+ & / ' C current airfoil ', A)
+ 1200 FORMAT( 1X,I2,' polar airfoil ', A)
+C
+ 2100 FORMAT(/' Select source of airfoil "',I1, A, $)
+C
+ DO 40 K = 1, 2
+ IAIR = K - 1
+ 20 WRITE(*,2100) IAIR, PROMPTN(1:NPR)
+ READ(*,1000,ERR=90,END=90) CAIR
+C
+ IF(CAIR .EQ. ' ') THEN
+ GO TO 90
+C
+ ELSEIF(INDEX('Ff',CAIR(1:1)).NE.0) THEN
+ CALL ASKS('Enter filename^',FNAME)
+ CALL AREAD(LU,FNAME,IBX,
+ & XINT(1,K),YINT(1,K),NINT(K),
+ & NAMEINT(K),ISPARST,ITYPE,0)
+ IF(ITYPE.EQ.0) RETURN
+C
+ ELSEIF(INDEX('Cc',CAIR(1:1)).NE.0) THEN
+ IF(N.LE.1) THEN
+ WRITE(*,*) 'No current airfoil available'
+ GO TO 20
+ ENDIF
+C
+ NINT(K) = N
+ DO I = 1, N
+ XINT(I,K) = X(I)
+ YINT(I,K) = Y(I)
+ ENDDO
+ NAMEINT(K) = NAME
+C
+ ELSE
+ READ(CAIR,*,ERR=90) IP
+ IF(IP.LT.1 .OR. IP.GT.NPOL) THEN
+ GO TO 90
+ ELSEIF(NXYPOL(IP).LE.0) THEN
+ GO TO 90
+ ELSE
+ NINT(K) = NXYPOL(IP)
+ DO I = 1, N
+ XINT(I,K) = CPOLXY(I,1,IP)
+ YINT(I,K) = CPOLXY(I,2,IP)
+ ENDDO
+ ENDIF
+ NAMEINT(K) = NAMEPOL(IP)
+C
+ ENDIF
+C
+ CALL SCALC(XINT(1,K),YINT(1,K),SINT(1,K),NINT(K))
+ CALL SEGSPLD(XINT(1,K),XPINT(1,K),SINT(1,K),NINT(K),-999.,-999.)
+ CALL SEGSPLD(YINT(1,K),YPINT(1,K),SINT(1,K),NINT(K),-999.,-999.)
+ CALL LEFIND(SLEINT(K),
+ & XINT(1,K),XPINT(1,K),
+ & YINT(1,K),YPINT(1,K),SINT(1,K),NINT(K))
+ 40 CONTINUE
+C
+ WRITE(*,*)
+ WRITE(*,*) 'airfoil "0": ', NAMEINT(1)
+ WRITE(*,*) 'airfoil "1": ', NAMEINT(2)
+ FRAC = 0.5
+ CALL ASKR('Specify interpolating fraction 0...1^',FRAC)
+C
+ CALL INTERX(XINT(1,1),XPINT(1,1),
+ & YINT(1,1),YPINT(1,1),SINT(1,1),NINT(1),SLEINT(1),
+ & XINT(1,2),XPINT(1,2),
+ & YINT(1,2),YPINT(1,2),SINT(1,2),NINT(2),SLEINT(2),
+ & XB,YB,NB,FRAC)
+C
+ CALL SCALC(XB,YB,SB,NB)
+ CALL SEGSPL(XB,XBP,SB,NB)
+ CALL SEGSPL(YB,YBP,SB,NB)
+C
+ CALL GEOPAR(XB,XBP,YB,YBP,SB,NB, W1,
+ & SBLE,CHORDB,AREAB,RADBLE,ANGBTE,
+ & EI11BA,EI22BA,APX1BA,APX2BA,
+ & EI11BT,EI22BT,APX1BT,APX2BT,
+ & THICKB,CAMBRB )
+C
+ CALL ASKS('Enter new airfoil name^',NAME)
+ CALL STRIP(NAME,NNAME)
+ WRITE(*,*)
+ WRITE(*,*) 'Result has been placed in buffer airfoil'
+ WRITE(*,*) 'Execute PCOP or PANE to set new current airfoil'
+ RETURN
+C
+ 90 CONTINUE
+ WRITE(*,*)
+ WRITE(*,*) 'Invalid response. No action taken.'
+ RETURN
+ END ! INTX
+
+
+
+