aboutsummaryrefslogtreecommitdiff
path: root/src/xgdes.f
diff options
context:
space:
mode:
Diffstat (limited to 'src/xgdes.f')
-rw-r--r--src/xgdes.f2314
1 files changed, 2314 insertions, 0 deletions
diff --git a/src/xgdes.f b/src/xgdes.f
new file mode 100644
index 0000000..d6c9d86
--- /dev/null
+++ b/src/xgdes.f
@@ -0,0 +1,2314 @@
+C***********************************************************************
+C Module: xgdes.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
+ SUBROUTINE GDES
+ INCLUDE 'XFOIL.INC'
+ CHARACTER*4 COMAND, COMOLD
+ LOGICAL LRECALC, LMODPL, LPLNEW
+ DIMENSION XBOX(2), YBOX(2), XRF(2)
+C
+ CHARACTER*128 COMARG, ARGOLD
+ CHARACTER*1 CHKEY
+C
+ DIMENSION IINPUT(20)
+ DIMENSION RINPUT(20)
+ LOGICAL ERROR
+C
+ EXTERNAL NEWPLOTG
+C
+ SAVE COMOLD, ARGOLD
+C
+ COMAND = '****'
+ COMARG = ' '
+ LRECALC = .FALSE.
+C
+ IF(NB.EQ.0) THEN
+ WRITE(*,*)
+ WRITE(*,*) '*** No airfoil available ***'
+ RETURN
+ ENDIF
+C
+ LPLCAM = .FALSE.
+ LSYM = .TRUE.
+C
+ WRITE(*,*)
+ WRITE(*,*) 'You are working with the buffer airfoil'
+C
+ CALL PLTINI
+ CALL GOFINI
+ CALL PLOTG
+C
+C====================================================
+C---- start of menu loop
+ 500 CONTINUE
+ COMOLD = COMAND
+ ARGOLD = COMARG
+C
+ 501 IF(LGSYM) THEN
+ CALL ASKC('.GDESs^',COMAND,COMARG)
+ ELSE
+ CALL ASKC('.GDES^',COMAND,COMARG)
+ ENDIF
+C
+C--------------------------------------------------------
+C---- process previous command ?
+ IF(COMAND(1:1).EQ.'!') THEN
+ IF(COMOLD.EQ.'****') THEN
+ WRITE(*,*) 'Previous .GDES command not valid'
+ GO TO 501
+ ELSE
+ COMAND = COMOLD
+ COMARG = ARGOLD
+ LRECALC = .TRUE.
+ ENDIF
+ ELSE
+ LRECALC = .FALSE.
+ ENDIF
+C
+ IF(COMAND.EQ.' ') THEN
+C----- just <return> was typed... clean up plotting and exit OPER
+ IF(LPLOT) CALL PLEND
+ LPLOT = .FALSE.
+ LGSYM = .FALSE.
+ LGEOPL = .FALSE.
+ IF(.NOT.LGSAME) THEN
+ WRITE(*,*)
+ WRITE(*,*) 'Buffer airfoil is not identical to current airfoil'
+ ENDIF
+ CALL CLRZOOM
+ RETURN
+ ENDIF
+C
+C---- extract command line numeric arguments
+ DO I=1, 20
+ IINPUT(I) = 0
+ RINPUT(I) = 0.0
+ ENDDO
+ NINPUT = 20
+ CALL GETINT(COMARG,IINPUT,NINPUT,ERROR)
+ NINPUT = 20
+ CALL GETFLT(COMARG,RINPUT,NINPUT,ERROR)
+C
+C--------------------------------------------------------
+ IF(COMAND.EQ.'? ') THEN
+ WRITE(*,1050)
+ 1050 FORMAT(
+ & /' <cr> Return to Top Level'
+ & /' ! Redo previous command'
+ &//' GSET Set buffer airfoil <== current airfoil'
+ & /' eXec Set current airfoil <== buffer airfoil'
+ & /' SYMM Toggle y-symmetry flag'
+ &//' ADEG r Rotate about origin (degrees)'
+ & /' ARAD r Rotate about origin (radians)'
+ & /' Tran rr Translate'
+ & /' Scal r Scale about origin'
+ & /' LINS rr. Linearly-varying y scale'
+ & /' DERO Derotate (set chord line level)'
+ &//' TGAP rr Change trailing edge gap'
+ & /' LERA rr Change leading edge radius'
+ &//' TCPL Toggle thickness and camber plotting'
+ & /' TFAC rr Scale existing thickness and camber'
+ & /' TSET rr Set new thickness and camber'
+ & /' HIGH rr Move camber and thickness highpoints'
+ & /' .CAMB Modify camber shape directly or via loading'
+ &//' BEND Display structural properties of buffer airfoil'
+ &//' Flap rrr Deflect trailing edge flap'
+ &//' Modi Modify contour via cursor'
+ & /' SLOP Toggle modified-contour slope matching flag'
+ &//' CORN Double point with cursor (set sharp corner)'
+ & /' ADDP Add point with cursor or keyboard x,y'
+ & /' MOVP Move point with cursor or keyboard x,y'
+ & /' DELP Delete point with cursor'
+ &//' UNIT Normalize buffer airfoil to unit chord'
+ & /' Dist Determine distance between 2 cursor points'
+ & /' CLIS List curvatures'
+ & /' CPLO Plot curvatures'
+ & /' CANG List panel corner angles'
+ & /' CADD ri. Add points at corners exceeding angle threshold'
+ &//' Plot Replot buffer airfoil'
+ & /' INPL Replot buffer airfoil without scaling (in inches)'
+ & /' Blow Blowup plot region'
+ & /' Rese Reset plot scale and origin'
+ & /' Wind Plot window adjust via cursor and keys'
+ &//' TSIZ r Change tick-mark size'
+ & /' TICK Toggle node tick-mark plotting'
+ & /' GRID Toggle grid plotting'
+ & /' GPAR Toggle geometric parameter plotting'
+ & /' Over f Overlay disk file airfoil'
+ &//' SIZE r Change absolute plot-object size'
+ & /' .ANNO Annotate plot'
+ & /' HARD Hardcopy current plot'
+ &//' NAME s Specify new airfoil name'
+ & /' NINC Increment name version number')
+
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'Z ') THEN
+ CALL USETZOOM(.TRUE.,.TRUE.)
+ CALL REPLOT(IDEV)
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'U ') THEN
+ CALL CLRZOOM
+ CALL REPLOT(IDEV)
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'GSET') THEN
+ NB = N
+ DO I=1, NB
+ XB(I) = X(I)
+ YB(I) = Y(I)
+ ENDDO
+ LGSAME = .TRUE.
+ 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 PLTINI
+ CALL PLOTG
+ IF(LGSYM) CALL ZERCAM
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'EXEC' .OR.
+ & COMAND.EQ.'X ' ) THEN
+ CALL ABCOPY(.TRUE.)
+cc CALL NAMMOD(NAME,1,1)
+cc CALL STRIP(NAME,NNAME)
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'SYMM') THEN
+ LGSYM = .NOT.LGSYM
+ IF(LGSYM) THEN
+ WRITE(*,*) 'y-symmetry forcing enabled.'
+ CALL ZERCAM
+ ELSE
+ WRITE(*,*) 'y-symmetry forcing disabled.'
+ ENDIF
+C
+C=================================================
+C---- rotate airfoil by degrees
+ ELSEIF(COMAND.EQ.'ADEG' .OR.
+ & COMAND.EQ.'ARAD' ) THEN
+ IF(COMAND.EQ.'ADEG') THEN
+ IF(NINPUT.GE.1) THEN
+ ADEG = RINPUT(1)
+ ELSE
+ ADEG = 0.0
+ CALL ASKR('Enter angle change (deg)^',ADEG)
+ ENDIF
+ ARAD = ADEG*PI/180.0
+ ELSE
+ IF(NINPUT.GE.1) THEN
+ ARAD = RINPUT(1)
+ ELSE
+ ARAD = 0.0
+ CALL ASKR('Enter angle change (rad)^',ARAD)
+ ENDIF
+ ENDIF
+C
+ CALL ROTATE(XB,YB,NB,ARAD)
+CCC CALL SCALC(XB,YB,SB,NB)
+ CALL SEGSPL(XB,XBP,SB,NB)
+ CALL SEGSPL(YB,YBP,SB,NB)
+C
+ APX1BA = APX1BA - ARAD
+ APX2BA = APX2BA - ARAD
+ APX1BT = APX1BT - ARAD
+ APX2BT = APX2BT - ARAD
+C
+ CALL NEWPEN(2)
+ CALL PLTAIR(XB,XBP,YB,YBP,SB,NB, XOFF,XSF,YOFF,YSF,'magenta')
+ CALL PLNEWP('magenta')
+ LGEOPL = .FALSE.
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'TRAN' .OR.
+ & COMAND.EQ.'T ' ) THEN
+ IF (NINPUT.GE.2) THEN
+ DELX = RINPUT(1)
+ DELY = RINPUT(2)
+ ELSEIF(NINPUT.GE.1) THEN
+ DELX = RINPUT(1)
+ DELY = 0.0
+ CALL ASKR('Enter delta(y)^',DELY)
+ ELSE
+ DELX = 0.0
+ CALL ASKR('Enter delta(x)^',DELX)
+ DELY = 0.0
+ CALL ASKR('Enter delta(y)^',DELY)
+ ENDIF
+ DO I=1, NB
+ XB(I) = XB(I) + DELX
+ YB(I) = YB(I) + DELY
+ ENDDO
+C
+ CALL NEWPEN(2)
+ CALL PLTAIR(XB,XBP,YB,YBP,SB,NB, XOFF,XSF,YOFF,YSF,'magenta')
+ CALL PLNEWP('magenta')
+ LGEOPL = .FALSE.
+ LGSAME = .FALSE.
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'SCAL' .OR.
+ & COMAND.EQ.'S ' ) THEN
+ IF(NINPUT.GE.1) THEN
+ FAC = RINPUT(1)
+ XXFAC = FAC
+ YYFAC = FAC
+ ELSE
+ FAC = 1.0
+ CALL ASKR('Enter scale factor (0 for separate x,y scales)^',FAC)
+ XXFAC = FAC
+ YYFAC = FAC
+ ENDIF
+C
+ IF(FAC .EQ. 0.0) THEN
+ IF(NINPUT.GE.3) THEN
+ XXFAC = RINPUT(2)
+ YYFAC = RINPUT(3)
+ ELSE
+ XXFAC = 1.0
+ CALL ASKR('Enter x scale factor^',XXFAC)
+ YYFAC = 1.0
+ CALL ASKR('Enter y scale factor^',YYFAC)
+ ENDIF
+ ENDIF
+C
+ DO I=1, NB
+ XB(I) = XB(I)*XXFAC
+ YB(I) = YB(I)*YYFAC
+ ENDDO
+C
+C----- re-order if necessary to maintain counterclockwise ordering
+ IF(XXFAC*YYFAC .LT. 0.0) THEN
+ DO I=1, NB/2
+ XTMP = XB(I)
+ YTMP = YB(I)
+ XB(I) = XB(NB-I+1)
+ YB(I) = YB(NB-I+1)
+ XB(NB-I+1) = XTMP
+ YB(NB-I+1) = YTMP
+ ENDDO
+ ENDIF
+C
+C----- re-spline new geometry
+ 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 NEWPEN(2)
+ CALL PLTAIR(XB,XBP,YB,YBP,SB,NB, XOFF,XSF,YOFF,YSF,'magenta')
+ CALL PLNEWP('magenta')
+ LGEOPL = .FALSE.
+ LGSAME = .FALSE.
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'LINS') THEN
+ 40 CONTINUE
+ IF(NINPUT.GE.4) THEN
+ XOC1 = RINPUT(1)
+ YFAC1 = RINPUT(2)
+ XOC2 = RINPUT(3)
+ YFAC2 = RINPUT(4)
+ ELSE
+ 1001 FORMAT(/1X,A,$)
+ 41 WRITE(*,1001) 'Location 1... enter x/c, y-scale : '
+ READ(*,*,ERR=41) XOC1, YFAC1
+ 42 WRITE(*,1001) 'Location 2... enter x/c, y-scale : '
+ READ(*,*,ERR=42) XOC2, YFAC2
+ ENDIF
+C
+ IF(ABS(XOC1-XOC2) .LT. 1.0E-5) THEN
+ WRITE(*,*) 'x/c locations 1 and 2 must be different'
+ NINPUT = 0
+ GO TO 40
+ ENDIF
+C
+ CALL LEFIND(SBLE,XB,XBP,YB,YBP,SB,NB)
+ XLE = SEVAL(SBLE,XB,XBP,SB,NB)
+ YLE = SEVAL(SBLE,YB,YBP,SB,NB)
+ XTE = 0.5*(XB(1) + XB(NB))
+ YTE = 0.5*(YB(1) + YB(NB))
+ DO I=1, NB
+ XOC = (XB(I)-XLE) / (XTE-XLE)
+ FR1 = (XOC2-XOC )/(XOC2-XOC1)
+ FR2 = (XOC -XOC1)/(XOC2-XOC1)
+ YYFAC = FR1*YFAC1 + FR2*YFAC2
+ YB(I) = YB(I)*YYFAC
+ ENDDO
+C
+C----- re-spline new geometry
+ 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 NEWPEN(2)
+ CALL PLTAIR(XB,XBP,YB,YBP,SB,NB, XOFF,XSF,YOFF,YSF,'magenta')
+ CALL PLNEWP('magenta')
+ LGEOPL = .FALSE.
+ LGSAME = .FALSE.
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'DERO') THEN
+ CALL LEFIND(SBLE,XB,XBP,YB,YBP,SB,NB)
+ XLE = SEVAL(SBLE,XB,XBP,SB,NB)
+ YLE = SEVAL(SBLE,YB,YBP,SB,NB)
+ XTE = 0.5*(XB(1) + XB(NB))
+ YTE = 0.5*(YB(1) + YB(NB))
+C
+ ARAD = ATAN2(YTE-YLE,XTE-XLE)
+ CALL ROTATE(XB,YB,NB,ARAD)
+ WRITE(*,1080) ARAD / DTOR
+ 1080 FORMAT(/'Rotating buffer airfoil by ',F8.3,' deg.')
+C
+ CALL SCALC(XB,YB,SB,NB)
+ CALL SEGSPL(XB,XBP,SB,NB)
+ CALL SEGSPL(YB,YBP,SB,NB)
+ 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 NEWPEN(2)
+ CALL PLTAIR(XB,XBP,YB,YBP,SB,NB, XOFF,XSF,YOFF,YSF,'magenta')
+ CALL PLNEWP('magenta')
+ LGEOPL = .FALSE.
+ LGSAME = .FALSE.
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'TGAP') THEN
+ CALL TGAP(RINPUT,NINPUT)
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'LERA') THEN
+ CALL LERAD(RINPUT,NINPUT)
+C
+C--------------------------------------------------------
+cc ELSEIF(COMAND.EQ.'TC ') THEN
+cc CALL TCBUF
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'TCPL') THEN
+ LPLCAM = .NOT.LPLCAM
+ CALL PLTINI
+ CALL GOFINI
+ CALL PLOTG
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'TFAC') THEN
+ IF(.NOT.LPLCAM) THEN
+ WRITE(*,*) 'Enabling camber,thickness plotting'
+ LPLCAM = .TRUE.
+ CALL PLTINI
+ CALL GOFINI
+ CALL PLOTG
+ ENDIF
+ CALL TCSCAL(RINPUT,NINPUT)
+ CALL NEWPEN(2)
+ CALL PLTAIR(XB,XBP,YB,YBP,SB,NB, XOFF,XSF,YOFF,YSF,'magenta')
+ CALL PLNEWP('magenta')
+ CALL PLTCAM('magenta')
+ LGEOPL = .FALSE.
+ LGSAME = .FALSE.
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'TSET') THEN
+ IF(.NOT.LPLCAM) THEN
+ WRITE(*,*) 'Enabling camber,thickness plotting'
+ LPLCAM = .TRUE.
+ CALL PLTINI
+ CALL GOFINI
+ CALL PLOTG
+ ENDIF
+ CALL TCSET(RINPUT,NINPUT)
+ CALL NEWPEN(2)
+ CALL PLTAIR(XB,XBP,YB,YBP,SB,NB, XOFF,XSF,YOFF,YSF,'magenta')
+ CALL PLNEWP('magenta')
+ CALL PLTCAM('magenta')
+ LGEOPL = .FALSE.
+ LGSAME = .FALSE.
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'HIGH') THEN
+ IF(.NOT.LPLCAM) THEN
+ WRITE(*,*) 'Enabling camber,thickness plotting'
+ LPLCAM = .TRUE.
+ CALL PLTINI
+ CALL GOFINI
+ CALL PLOTG
+ ENDIF
+ CALL HIPNT(RINPUT,NINPUT)
+ CALL NEWPEN(2)
+ CALL PLTAIR(XB,XBP,YB,YBP,SB,NB, XOFF,XSF,YOFF,YSF,'magenta')
+ CALL PLNEWP('magenta')
+ CALL PLTCAM('magenta')
+ LGEOPL = .FALSE.
+ LGSAME = .FALSE.
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'CAMB') THEN
+ IF(LGSYM) THEN
+ WRITE(*,*) 'Disabling symmetry enforcement.'
+ LGSYM = .FALSE.
+ ENDIF
+ CALL CAMB
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'BEND') THEN
+ CALL BENDUMP(NB,XB,YB)
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'CANG') THEN
+ CALL CANG(XB,YB,NB,2, IMAX,AMAX)
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'CADD') THEN
+ CALL CANG(XB,YB,NB,2, IMAX,AMAX)
+ WRITE(*,*)
+C
+ XBMIN = XB(1)
+ XBMAX = XB(1)
+ DO I=1, NB
+ XBMIN = MIN(XBMIN,XB(I))
+ XBMAX = MAX(XBMAX,XB(I))
+ ENDDO
+C
+C----- default inputs
+ ATOL = 0.5*AMAX
+ ISPL = 1
+ XRF(1) = XBMIN - 0.1*(XBMAX-XBMIN)
+ XRF(2) = XBMAX + 0.1*(XBMAX-XBMIN)
+C
+ IF (NINPUT.LE.0) THEN
+ GO TO 70
+ ELSEIF(NINPUT.LE.1) THEN
+ ATOL = RINPUT(1)
+ GO TO 71
+ ELSEIF(NINPUT.LE.2) THEN
+ ATOL = RINPUT(1)
+ ISPL = IINPUT(2)
+ GO TO 72
+ ELSEIF(NINPUT.LE.4) THEN
+ ATOL = RINPUT(1)
+ ISPL = IINPUT(2)
+ XRF(1) = RINPUT(3)
+ XRF(2) = RINPUT(4)
+ GO TO 74
+ ENDIF
+C
+ 70 WRITE(*,1090) ATOL
+ 1090 FORMAT(1X,
+ & 'Enter corner angle criterion for refinement (deg):', F8.3)
+ CALL READR(1,ATOL,ERROR)
+ IF(ERROR) GO TO 70
+C
+ 71 WRITE(*,1091) ISPL
+ 1091 FORMAT(1X,
+ & 'Enter type of spline parameter (1=uniform, 2=arclength):', I4)
+ CALL READI(1,ISPL,ERROR)
+ IF(ERROR) GO TO 71
+ IF(ISPL.LE.0) GO TO 500
+ IF(ISPL.GT.2) GO TO 71
+C
+ 72 WRITE(*,1092) XRF(1), XRF(2)
+ 1092 FORMAT(1X,
+ & 'Enter refinement x limits:', 2F10.5)
+ CALL READR(2,XRF,ERROR)
+ IF(ERROR) GO TO 72
+C
+ 74 CONTINUE
+ IF(ISPL.EQ.1) THEN
+ SB(1) = 0.0
+ DO I = 2, NB
+ IF(XB(I).EQ.XB(I-1) .AND. YB(I).EQ.YB(I-1)) THEN
+ SB(I) = SB(I-1)
+ ELSE
+ SB(I) = SB(I-1) + 1.0
+ ENDIF
+ ENDDO
+ CALL SEGSPL(XB,XBP,SB,NB)
+ CALL SEGSPL(YB,YBP,SB,NB)
+ ENDIF
+C
+ CALL AREFINE(XB,YB,SB,XBP,YBP,NB, ATOL,
+ & IBX,NNEW,W1,W2,XRF(1),XRF(2))
+C
+ NBADD = NNEW - NB
+ WRITE(*,*) 'Number of points added: ', NBADD
+C
+ NB = NNEW
+ DO I = 1, NB
+ XB(I) = W1(I)
+ YB(I) = W2(I)
+ ENDDO
+ LGSAME = .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
+ CALL NEWPEN(2)
+ CALL PLTAIR(XB,XBP,YB,YBP,SB,NB, XOFF,XSF,YOFF,YSF,'magenta')
+ CALL PLNEWP('magenta')
+C
+ CALL CANG(XB,YB,NB,1, IMAX,AMAX)
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'CLIS') THEN
+ CALL CLIS(XB,XBP,YB,YBP,SB,NB)
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'CPLO') THEN
+ CALL PLTCRV(SBLE,XB,XBP,YB,YBP,SB,NB,W1)
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'FLAP' .OR.
+ & COMAND.EQ.'F ' ) THEN
+ CALL FLAP(RINPUT,NINPUT)
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'MODI' .OR.
+ & COMAND.EQ.'M ' ) THEN
+C----- plot current geometry if it's not on the screen
+ IF(.NOT.LGEOPL) THEN
+ CALL PLTINI
+ CALL PLOTG
+ ENDIF
+C
+ IF(LGSYM) THEN
+ DO I = 1, NB
+ W1(I) = XB(I)
+ W2(I) = YB(I)
+ ENDDO
+ ENDIF
+C
+ IBFRST = 1
+ IBLAST = NB
+ NSIDE = 1
+ XBOX(1) = XMARG
+ XBOX(2) = XPAGE-XMARG
+ YBOX(1) = YMARG
+ YBOX(2) = YPAGE-YMARG
+ LMODPL = .FALSE.
+ CALL MODIXY(IBX,IBFRST,IBLAST,NSIDE,
+ & XB,YB,XBP,YBP,SB, LGSLOP,
+ & IGMOD1,IGMOD2,ISMOD,
+ & XBOX,YBOX, XBOX,YBOX,SIZE,
+ & XOFF,YOFF,XSF,YSF, LMODPL,
+ & NEWPLOTG)
+C
+ IF(LGSYM) THEN
+ DO I = 1, NB
+ XBDEL = XB(I) - W1(I)
+ YBDEL = YB(I) - W2(I)
+ XB(I) = XB(I) + XBDEL
+ YB(I) = YB(I) + YBDEL
+ ENDDO
+ CALL ZERCAM
+ ENDIF
+ LGSAME = .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
+ CALL NEWPEN(2)
+ CALL PLTAIR(XB,XBP,YB,YBP,SB,NB, XOFF,XSF,YOFF,YSF,'magenta')
+ CALL PLNEWP('magenta')
+ LGEOPL = .FALSE.
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'SLOP') THEN
+ LGSLOP = .NOT.LGSLOP
+ IF(LGSLOP) THEN
+ WRITE(*,*) 'Modified segment will be',
+ & ' made tangent at endpoints'
+ ELSE
+ WRITE(*,*) 'Modified segment will not be',
+ & ' made tangent at endpoints'
+ ENDIF
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'CORN') THEN
+ IF(NB.EQ.2*IQX) THEN
+ WRITE(*,*)
+ & 'Buffer airfoil arrays will overflow. No action taken.'
+ GO TO 500
+ ENDIF
+C
+ XWS = XWIND/SIZE
+ YWS = YWIND/SIZE
+ CALL POINTF(XB,XBP,YB,YBP,SB,NB, XWS,YWS, XOFF,YOFF,XSF,YSF,
+ & IPNT,XC,YC)
+ IF(IPNT.EQ.0) GO TO 500
+ IF(IPNT.EQ.1 .OR. IPNT.EQ.NB) THEN
+ WRITE(*,*) 'Cannot double trailing edge point. No action taken.'
+ GO TO 500
+ ENDIF
+C
+C----- add doubled point
+ DO I=NB, IPNT, -1
+ XB(I+1) = XB(I)
+ YB(I+1) = YB(I)
+ ENDDO
+ NB = NB+1
+ LGSAME = .FALSE.
+C
+C----- spline new geometry
+ CALL SCALC(XB,YB,SB,NB)
+ CALL SEGSPL(XB,XBP,SB,NB)
+ CALL SEGSPL(YB,YBP,SB,NB)
+ CALL PLTINI
+ CALL PLOTG
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'ADDP') THEN
+ CALL ADDP
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'DELP') THEN
+ CALL DELP
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'MOVP') THEN
+ CALL MOVP(NEWPLOTG)
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'UNIT') THEN
+ CALL NORM(XB,XBP,YB,YBP,SB,NB)
+ LGSAME = .FALSE.
+C
+C----- re-spline new geometry
+ 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 NEWPEN(2)
+ CALL PLTAIR(XB,XBP,YB,YBP,SB,NB, XOFF,XSF,YOFF,YSF,'magenta')
+ CALL PLNEWP('magenta')
+ LGEOPL = .FALSE.
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'DIST' .OR.
+ & COMAND.EQ.'D ' ) THEN
+ CALL DIST
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'HARD') THEN
+ IF(LPLOT) CALL PLEND
+ LPLOT = .FALSE.
+ CALL REPLOT(IDEVRP)
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'PLOT' .OR.
+ & COMAND.EQ.'P ' ) THEN
+ CALL PLTINI
+ CALL PLOTG
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'INPL') THEN
+ CALL PLTINI
+ XOFF0 = XOFF
+ YOFF0 = YOFF
+ XSF0 = XSF
+ YSF0 = YSF
+C
+ XSF = 1.0/SIZE
+ YSF = 1.0/SIZE
+c write(*,*) 'Enter Xoff, Yoff'
+c read (*,*) xoff, yoff
+c xoff = -xoff
+c yoff = -yoff
+c
+ CALL PLOTG
+ XOFF = XOFF0
+ YOFF = YOFF0
+ XSF = XSF0
+ YSF = YSF0
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'BLOW' .OR.
+ & COMAND.EQ.'B ' ) THEN
+ XWS = XWIND/SIZE
+ YWS = YWIND/SIZE
+ CALL OFFGET(XOFF,YOFF,XSF,YSF,XWS,YWS, .TRUE. , .TRUE. )
+ CALL GOFSET
+ CALL PLTINI
+ CALL PLOTG
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'RESE' .OR.
+ & COMAND.EQ.'R ' ) THEN
+ CALL PLTINI
+ CALL GOFINI
+ CALL PLOTG
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'WIND' .OR.
+ & COMAND.EQ.'W ' ) THEN
+ XWS = XWIND/SIZE
+ YWS = YWIND/SIZE
+C
+ WRITE(*,*) ' '
+ WRITE(*,*) 'Type I,O,P to In,Out,Pan with cursor...'
+C
+ 80 CALL PLTINI
+ CALL PLOTG
+C
+ CALL GETCURSORXY(XCRS,YCRS,CHKEY)
+C
+C----- do possible pan,zoom operations based on CHKEY
+ CALL KEYOFF(XCRS,YCRS,CHKEY, XWS,YWS, XOFF,YOFF,XSF,YSF, LPLNEW)
+C
+ IF(LPLNEW) THEN
+ CALL GOFSET
+ GO TO 80
+ ENDIF
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'TSIZ') THEN
+ IF(NINPUT.GE.1) THEN
+ GTICK = RINPUT(1)
+ ELSE
+ WRITE(*,*)
+ & 'Current tick-mark size (as fraction of perimeter) =', GTICK
+ CALL ASKR('Enter new tick-mark size^',GTICK)
+ ENDIF
+ CALL PLTINI
+ CALL PLOTG
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'TICK') THEN
+ LGTICK = .NOT.LGTICK
+ CALL PLTINI
+ CALL PLOTG
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'GRID') THEN
+ LGGRID = .NOT.LGGRID
+ CALL PLTINI
+ CALL PLOTG
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'GPAR') THEN
+ LGPARM = .NOT.LGPARM
+ CALL PLTINI
+ CALL PLOTG
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'SIZE') THEN
+ IF(NINPUT.GE.1) THEN
+ SIZE = RINPUT(1)
+ ELSE
+ WRITE(*,*) 'Current plot-object size =', SIZE
+ CALL ASKR('Enter new plot-object size^',SIZE)
+ ENDIF
+ CALL PLTINI
+ CALL PLOTG
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'OVER' .OR.
+ & COMAND.EQ.'O ' ) THEN
+ CALL OVER(COMARG)
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'ANNO') THEN
+ IF(LPLOT) THEN
+ CALL ANNOT(CH)
+ ELSE
+ WRITE(*,*) 'No active plot to annotate'
+ ENDIF
+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.'NDEC') THEN
+ CALL NAMMOD(NAME,-1,1)
+ CALL STRIP(NAME,NNAME)
+C
+C--------------------------------------------------------
+ ELSEIF(COMAND.EQ.'SINT') THEN
+ CALL SPLNXY(XB,XBP,YB,YBP,SB,NB)
+ CALL PLTAIR(XB,XBP,YB,YBP,SB,NB, XOFF,XSF,YOFF,YSF,'cyan')
+ CALL NEWPEN(2)
+ CALL PLTAIR(XB,XBP,YB,YBP,SB,NB, XOFF,XSF,YOFF,YSF,'magenta')
+ CALL PLNEWP('magenta')
+ LGEOPL = .FALSE.
+C
+C--------------------------------------------------------
+ ELSE
+ WRITE(*,1100) COMAND
+ 1100 FORMAT(' Command ',A4,' not recognized. Type a " ? " for list.')
+ COMAND = '****'
+ ENDIF
+C
+ GO TO 500
+ END ! GDES
+
+
+ SUBROUTINE NEWPLOTG
+ CALL GOFSET
+ CALL PLTINI
+ CALL PLOTG
+ RETURN
+ END
+
+
+ SUBROUTINE ABCOPY(LCONF)
+ INCLUDE 'XFOIL.INC'
+ LOGICAL LCONF
+C
+ IF(NB.LE.1) THEN
+ WRITE(*,*) 'ABCOPY: Buffer airfoil not available.'
+ RETURN
+ ELSEIF(NB.GT.IQX-5) THEN
+ WRITE(*,*) 'Maximum number of panel nodes : ',IQX-5
+ WRITE(*,*) 'Number of buffer airfoil points: ',NB
+ WRITE(*,*) 'Current airfoil cannot be set.'
+ WRITE(*,*) 'Try executing PANE at Top Level instead.'
+ RETURN
+ ENDIF
+ IF(N.NE.NB) LBLINI = .FALSE.
+C
+ N = NB
+ DO 101 I=1, N
+ X(I) = XB(I)
+ Y(I) = YB(I)
+ 101 CONTINUE
+ LGSAME = .TRUE.
+C
+ IF(LBFLAP) THEN
+ XOF = XBF
+ YOF = YBF
+ LFLAP = .TRUE.
+ ENDIF
+C
+C---- strip out doubled points
+ I = 1
+ 102 CONTINUE
+ I = I+1
+ IF(X(I-1).EQ.X(I) .AND. Y(I-1).EQ.Y(I)) THEN
+ DO 104 J=I, N-1
+ X(J) = X(J+1)
+ Y(J) = Y(J+1)
+ 104 CONTINUE
+ N = N-1
+ ENDIF
+ IF(I.LT.N) GO TO 102
+C
+ CALL SCALC(X,Y,S,N)
+ CALL SEGSPL(X,XP,S,N)
+ CALL SEGSPL(Y,YP,S,N)
+ CALL NCALC(X,Y,S,N,NX,NY)
+ CALL LEFIND(SLE,X,XP,Y,YP,S,N)
+ 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 )
+ CALL TECALC
+ CALL APCALC
+C
+ LGAMU = .FALSE.
+ LQINU = .FALSE.
+ LWAKE = .FALSE.
+ LQAIJ = .FALSE.
+ LADIJ = .FALSE.
+ LWDIJ = .FALSE.
+ LIPAN = .FALSE.
+ LVCONV = .FALSE.
+ LSCINI = .FALSE.
+CCC LBLINI = .FALSE.
+C
+ IF(LCONF) WRITE(*,1200) N
+ 1200 FORMAT(/' Current airfoil nodes set from buffer airfoil nodes (',
+ & I4,' )')
+C
+ RETURN
+ END ! ABCOPY
+
+
+ SUBROUTINE GOFINI
+C----------------------------------------------------------
+C Sets initial airfoil scaling and offset parameters
+C----------------------------------------------------------
+ INCLUDE 'XFOIL.INC'
+C
+C---- get airfoil bounding box
+ XBMIN = XB(1)
+ YBMIN = YB(1)
+ XBMAX = XB(1)
+ YBMAX = YB(1)
+ DO I=1, NB
+ XBMIN = MIN(XBMIN,XB(I))
+ YBMIN = MIN(YBMIN,YB(I))
+ XBMAX = MAX(XBMAX,XB(I))
+ YBMAX = MAX(YBMAX,YB(I))
+ ENDDO
+C
+C---- set camber and thickness distributions
+ CALL GETCAM(XCM,YCM,NCM,XTK,YTK,NTK,
+ & XB,XBP,YB,YBP,SB,NB )
+C
+C---- get camber,thickness y bounds
+ CMMIN = 0.
+ CMMAX = 0.
+ DO I=1, NCM
+ CMMIN = MIN(CMMIN,YCM(I))
+ CMMAX = MAX(CMMAX,YCM(I))
+ ENDDO
+ TKMIN = 0.
+ TKMAX = 0.
+ DO I=1, NTK
+ TKMIN = MIN(TKMIN,YTK(I))
+ TKMAX = MAX(TKMAX,YTK(I))
+ ENDDO
+C
+ XRANGE = XBMAX - XBMIN
+ YRANGE = YBMAX - YBMIN
+C
+C---- set x,y scaling factors needed for O(1) size plot with "nice" limits
+ CALL SCALIT(1,0.95*XRANGE,0.0,XSF)
+ CALL SCALIT(1,0.95*YRANGE,0.0,YSF)
+C
+C---- grid increment as a fraction of a nice upper bound on delta x
+cc DXYG = 0.1 / XSF
+ DXYG = 0.1 / MIN(XSF,YSF)
+C
+C---- set "nice" grid limits as integer multiples of DXYG
+c XGMAX = DXYG*(INT(XBMAX/DXYG+1000.05) - 999)
+c XGMIN = DXYG*(INT(XBMIN/DXYG-1000.05) + 999)
+c YGMAX = DXYG*(INT(YBMAX/DXYG+1000.25) - 999)
+c YGMIN = DXYG*(INT(YBMIN/DXYG-1000.25) + 999)
+C
+C---- set "nice" grid limits as integer multiples of DXYG
+ XGMAX = DXYG*(INT(XBMAX/DXYG+1001.01) - 1000)
+ XGMIN = DXYG*(INT(XBMIN/DXYG-1001.01) + 1000)
+ YGMAX = DXYG*(INT(YBMAX/DXYG+1001.01) - 1000)
+ YGMIN = DXYG*(INT(YBMIN/DXYG-1001.01) + 1000)
+C
+C---- set bounding box for thickness/camber plot
+ DXYC = DXYG
+ XCMIN = XGMIN
+ XCMAX = XGMAX
+ YCMIN = MIN(CMMIN,-TKMAX)
+ YCMAX = MAX(CMMAX, TKMAX)
+ YCMAX = DXYC*(INT(YCMAX/DXYC+1000.25) - 999)
+ YCMIN = DXYC*(INT(YCMIN/DXYC-1000.25) + 999)
+ YCMAX = MAX(YCMAX,YCMIN+DXYC)
+C
+C---- set minimum scaling factor to fit airfoil or grid
+ IF(LGGRID) THEN
+ XRANGE = XGMAX - XGMIN
+ YRANGE = YGMAX - YGMIN
+ ELSE
+ XRANGE = XBMAX - XBMIN
+ YRANGE = YBMAX - YBMIN
+ ENDIF
+C
+C---- include y range from thickness/camber plot if present
+ IF(LPLCAM) THEN
+ YRANGE = YRANGE + (YCMAX - YCMIN)
+ ENDIF
+C
+ RANGE = MAX(XRANGE,YRANGE)
+C
+ SF = MIN( 1.0/XRANGE , PLOTAR/YRANGE )
+ XSF = SF
+ YSF = SF
+ CHG = 0.75*CH * RANGE*SF
+C--- HHY 4/24/01 keep the character size from getting too low
+
+ CHG = MAX(CHG,0.0075)
+C
+ IF(LGGRID) THEN
+C------ set offsets to position grid, with space for numerical axis annotations
+ XOFF = XGMIN - 0.05*RANGE - 3.0*CHG/SF
+ YOFF = YGMIN - 0.05*RANGE - 2.0*CHG/SF
+ ELSE
+C------ set offsets to position airfoil
+ XOFF = XBMIN - 0.05*RANGE
+ YOFF = YBMIN - 0.05*RANGE
+ ENDIF
+C
+C---- set plot limits for DCp plot (y-axis limit defaults set in INIT)
+ XPMIN = XGMIN
+ XPMAX = XGMAX
+ccc DXYP = DXYG
+ CALL AXISADJ(YPMIN,YPMAX,PSPAN,DXYP,NTICS)
+C
+C---- set Yoffset for camber plot in scale factor YSF for geom plots
+ DYOFFC = - YGMAX + YCMIN - 2.2*CHG/YSF
+C
+C---- set the Cp scale factor for DCp plots
+ PAR = (YPAGE-2.0*YMARG)/(XPAGE-2.0*XMARG)
+ DPRANGE = YPMAX-YPMIN
+ DYPLT = MAX(0.1,PAR-PLOTAR)
+ YSFP = 0.8*DYPLT/DPRANGE
+ YSFP = YSFP/YSF
+C
+C---- set shifts to YOFF for DCp plots in scale factor YSF for geom plots
+ DYOFFP = -YCMAX+DYOFFC + YPMIN*YSFP - 2.2*CHG/YSF
+C
+ RETURN
+ END ! GOFINI
+
+
+
+ SUBROUTINE GOFSET
+C----------------------------------------------------------
+C Sets grid-overlay parameters
+C----------------------------------------------------------
+ INCLUDE 'XFOIL.INC'
+C
+C---- airfoil extent
+ XBMIN = XB(1)
+ YBMIN = YB(1)
+ XBMAX = XB(1)
+ YBMAX = YB(1)
+ DO I=1, NB
+ XBMIN = MIN(XBMIN,XB(I))
+ YBMIN = MIN(YBMIN,YB(I))
+ XBMAX = MAX(XBMAX,XB(I))
+ YBMAX = MAX(YBMAX,YB(I))
+ ENDDO
+C
+ RANGE = MAX( (XWIND/SIZE)/XSF , (YWIND/SIZE)/YSF )
+C
+C---- set bounding-box corner locations in user coordinates
+ XG1 = XOFF + 0.1*RANGE + 4.0*CHG/XSF
+ YG1 = YOFF + 0.1*RANGE + 2.0*CHG/YSF
+ XG2 = XOFF - 0.1*RANGE + (XWIND/SIZE)/XSF
+ YG2 = YOFF - 0.1*RANGE + (YWIND/SIZE)/YSF
+C
+C---- crunch down onto airfoil limits
+ XG1 = MAX(XG1,XBMIN)
+ XG2 = MIN(XG2,XBMAX)
+ YG1 = MAX(YG1,YBMIN)
+ YG2 = MIN(YG2,YBMAX)
+C
+C---- set x,y scaling factors needed for O(1) size plot with "nice" limits
+ CALL SCALIT(1,0.95*(XG2-XG1),0.0,GXSF)
+ CALL SCALIT(1,0.95*(YG2-YG1),0.0,GYSF)
+C
+ GSF = GXSF
+ccc GSF = MIN(GXSF,GYSF)
+C
+C---- grid increment as a fraction of a nice upper bound on delta x
+ DXYG = 0.1 / GSF
+C
+C---- set "nice" grid limits as integer multiples of DXYG
+ XGMAX = DXYG*(INT(XG2/DXYG+1001.01) - 1000)
+ XGMIN = DXYG*(INT(XG1/DXYG-1001.01) + 1000)
+ YGMAX = DXYG*(INT(YG2/DXYG+1001.01) - 1000)
+ YGMIN = DXYG*(INT(YG1/DXYG-1001.01) + 1000)
+C
+ RETURN
+ END ! GOFSET
+
+
+
+ SUBROUTINE TGAP(RINPUT,NINPUT)
+C----------------------------------
+C Used to set buffer airfoil
+C trailing edge gap
+C----------------------------------
+ INCLUDE 'XFOIL.INC'
+ DIMENSION RINPUT(*)
+C
+ CALL LEFIND(SBLE,XB,XBP,YB,YBP,SB,NB)
+ 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
+ DXN = XB(1) - XB(NB)
+ DYN = YB(1) - YB(NB)
+ GAP = SQRT(DXN**2 + DYN**2)
+C
+C---- components of unit vector parallel to TE gap
+ IF(GAP.GT.0.0) THEN
+ DXU = DXN / GAP
+ DYU = DYN / GAP
+ ELSE
+ DXU = -.5*(YBP(NB) - YBP(1))
+ DYU = 0.5*(XBP(NB) - XBP(1))
+ ENDIF
+C
+ IF (NINPUT .GE. 2) THEN
+ GAPNEW = RINPUT(1)
+ DOC = RINPUT(2)
+ ELSEIF(NINPUT .GE. 1) THEN
+ GAPNEW = RINPUT(1)
+ DOC = 1.0
+ CALL ASKR('Enter blending distance/c (0..1)^',DOC)
+ ELSE
+ WRITE(*,1000) GAP
+ 1000 FORMAT(/' Current gap =',F9.5)
+ GAPNEW = 0.0
+ CALL ASKR('Enter new gap^',GAPNEW)
+ DOC = 1.0
+ CALL ASKR('Enter blending distance/c (0..1)^',DOC)
+ ENDIF
+C
+ DOC = MIN( MAX( DOC , 0.0 ) , 1.0 )
+C
+ DGAP = GAPNEW - GAP
+C
+C---- go over each point, changing the y-thickness appropriately
+ DO 30 I=1, NB
+C
+C------ chord-based x/c
+ XOC = ( (XB(I)-XBLE)*(XBTE-XBLE)
+ & + (YB(I)-YBLE)*(YBTE-YBLE) ) / CHBSQ
+C
+C------ thickness factor tails off exponentially away from trailing edge
+ IF(DOC .EQ. 0.0) THEN
+ TFAC = 0.0
+ IF(I.EQ.1 .OR. I.EQ.NB) TFAC = 1.0
+ ELSE
+ ARG = MIN( (1.0-XOC)*(1.0/DOC-1.0) , 15.0 )
+ TFAC = EXP(-ARG)
+ ENDIF
+C
+ IF(SB(I).LE.SBLE) THEN
+ XB(I) = XB(I) + 0.5*DGAP*XOC*TFAC*DXU
+ YB(I) = YB(I) + 0.5*DGAP*XOC*TFAC*DYU
+ ELSE
+ XB(I) = XB(I) - 0.5*DGAP*XOC*TFAC*DXU
+ YB(I) = YB(I) - 0.5*DGAP*XOC*TFAC*DYU
+ ENDIF
+ 30 CONTINUE
+ LGSAME = .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
+ CALL PLTAIR(XB,XBP,YB,YBP,SB,NB, XOFF,XSF,YOFF,YSF,'magenta')
+ CALL PLNEWP('magenta')
+C
+ LGEOPL = .FALSE.
+C
+ RETURN
+ END ! TGAP
+
+
+
+ SUBROUTINE LERAD(RINPUT,NINPUT)
+C----------------------------
+C Changes buffer airfoil
+C leading edge radius.
+C----------------------------
+ INCLUDE 'XFOIL.INC'
+ DIMENSION RINPUT(*)
+C
+ IF (NINPUT .GE. 2) THEN
+ RFAC = RINPUT(1)
+ DOC = RINPUT(2)
+ ELSEIF(NINPUT .GE. 1) THEN
+ RFAC = RINPUT(1)
+ DOC = 1.0
+ CALL ASKR('Enter blending distance/c from LE^',DOC)
+ ELSE
+ RFAC = 1.0
+ CALL ASKR('Enter approx. new/old LE radius scaling ratio^',RFAC)
+ DOC = 1.0
+ CALL ASKR('Enter blending distance/c from LE^',DOC)
+ ENDIF
+C
+ DOC = MAX( DOC , 0.001 )
+C
+ CALL LERSCL(XB,XBP,YB,YBP,SB,NB, DOC,RFAC, W1,W2)
+C
+ DO 40 I=1, NB
+ XB(I) = W1(I)
+ YB(I) = W2(I)
+ 40 CONTINUE
+ LGSAME = .FALSE.
+C
+C---- spline new coordinates
+ 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---- find max curvature
+ CVMAX = 0.
+ DO 6 I=NB/4, (3*NB)/4
+ CV = CURV(SB(I),XB,XBP,YB,YBP,SB,NB)
+ CVMAX = MAX( ABS(CV) , CVMAX )
+ 6 CONTINUE
+C
+ RADIUS = 1.0/CVMAX
+C
+ WRITE(*,1000) RADIUS
+ 1000 FORMAT(/' New LE radius = ',F7.5)
+C
+ CALL PLTAIR(XB,XBP,YB,YBP,SB,NB, XOFF,XSF,YOFF,YSF,'magenta')
+ CALL PLNEWP('magenta')
+C
+ LGEOPL = .FALSE.
+C
+ RETURN
+ END ! LERAD
+
+
+
+ SUBROUTINE SCLXY
+C---------------------------------------------------
+C Scale airfoil about LE, TE, or selected point
+C---------------------------------------------------
+ INCLUDE 'XFOIL.INC'
+ CHARACTER*1 VAR
+C
+ CALL LEFIND(SBLE,XB,XBP,YB,YBP,SB,NB)
+ XLE = SEVAL(SBLE,XB,XBP,SB,NB)
+ YLE = SEVAL(SBLE,YB,YBP,SB,NB)
+ XTE = 0.5*(XB(1) + XB(NB))
+ YTE = 0.5*(YB(1) + YB(NB))
+C
+ WRITE(*,*) 'Enter origin for airfoil scaling:'
+ WRITE(*,*) ' L scales about LE'
+ WRITE(*,*) ' T scales about TE'
+ WRITE(*,*) ' P scales about input point'
+C
+ CALL ASKS('Select origin for scaling^',VAR)
+ IF (VAR.EQ.'L') THEN
+ XORG = XLE
+ YORG = YLE
+ ELSE IF (VAR.EQ.'T') THEN
+ XORG = XTE
+ YORG = YTE
+ ELSE
+ XORG = 0.25
+ YORG = 0.0
+ CALL ASKR('Enter X origin for scaling^',XORG)
+ CALL ASKR('Enter Y origin for scaling^',YORG)
+ ENDIF
+C
+ SCL = 1.0
+ CALL ASKR('Enter scaling factor about selected point^',SCL)
+C
+ DO 10 I=1, NB
+ XB(I) = SCL*(XB(I) - XORG) + XORG
+ YB(I) = SCL*(YB(I) - YORG) + YORG
+ 10 CONTINUE
+ LGSAME = .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
+ RETURN
+ END ! SCLXY
+
+
+
+ SUBROUTINE FLAP(RINPUT,NINPUT)
+C----------------------------------------------------
+C Modifies buffer airfoil for a deflected flap.
+C Points may be added/subtracted in the flap
+C break vicinity to clean things up.
+C----------------------------------------------------
+ INCLUDE 'XFOIL.INC'
+ LOGICAL LCHANGE
+ DIMENSION RINPUT(*)
+C
+ LOGICAL INSID
+ LOGICAL INSIDE
+ LOGICAL LT1NEW,LT2NEW,LB1NEW,LB2NEW
+C
+ SHT = CH * MAX(XSF,YSF)
+C
+ IF(NINPUT.GE.2) THEN
+ XBF = RINPUT(1)
+ YBF = RINPUT(2)
+ ELSE
+ XBF = -999.0
+ YBF = -999.0
+ ENDIF
+C
+ CALL GETXYF(XB,XBP,YB,YBP,SB,NB, TOPS,BOTS,XBF,YBF)
+ INSID = INSIDE(XB,YB,NB,XBF,YBF)
+C
+ WRITE(*,1050) XBF, YBF
+ 1050 FORMAT(/' Flap hinge: x,y =', 2F9.5 )
+C
+ IF(NINPUT.GE.3) THEN
+ DDEF = RINPUT(3)
+ ELSE
+ DDEF = 0.0
+ CALL ASKR('Enter flap deflection in degrees (+ down)^',DDEF)
+ ENDIF
+ RDEF = DDEF*PI/180.0
+ IF(RDEF .EQ. 0.0) RETURN
+C
+C
+ IF(INSID) THEN
+ ATOP = MAX( 0.0 , -RDEF )
+ ABOT = MAX( 0.0 , RDEF )
+ ELSE
+ CHX = DEVAL(BOTS,XB,XBP,SB,NB) - DEVAL(TOPS,XB,XBP,SB,NB)
+ CHY = DEVAL(BOTS,YB,YBP,SB,NB) - DEVAL(TOPS,YB,YBP,SB,NB)
+ FVX = SEVAL(BOTS,XB,XBP,SB,NB) + SEVAL(TOPS,XB,XBP,SB,NB)
+ FVY = SEVAL(BOTS,YB,YBP,SB,NB) + SEVAL(TOPS,YB,YBP,SB,NB)
+ CRSP = CHX*(YBF-0.5*FVY) - CHY*(XBF-0.5*FVX)
+ IF(CRSP .GT. 0.0) THEN
+C-------- flap hinge is above airfoil
+ ATOP = MAX( 0.0 , RDEF )
+ ABOT = MAX( 0.0 , RDEF )
+ ELSE
+C-------- flap hinge is below airfoil
+ ATOP = MAX( 0.0 , -RDEF )
+ ABOT = MAX( 0.0 , -RDEF )
+ ENDIF
+ ENDIF
+C
+C---- find upper and lower surface break arc length values...
+ CALL SSS(TOPS,ST1,ST2,ATOP,XBF,YBF,XB,XBP,YB,YBP,SB,NB,1)
+ CALL SSS(BOTS,SB1,SB2,ABOT,XBF,YBF,XB,XBP,YB,YBP,SB,NB,2)
+C
+C---- ... and x,y coordinates
+ XT1 = SEVAL(ST1,XB,XBP,SB,NB)
+ YT1 = SEVAL(ST1,YB,YBP,SB,NB)
+ XT2 = SEVAL(ST2,XB,XBP,SB,NB)
+ YT2 = SEVAL(ST2,YB,YBP,SB,NB)
+ XB1 = SEVAL(SB1,XB,XBP,SB,NB)
+ YB1 = SEVAL(SB1,YB,YBP,SB,NB)
+ XB2 = SEVAL(SB2,XB,XBP,SB,NB)
+ YB2 = SEVAL(SB2,YB,YBP,SB,NB)
+C
+C
+ WRITE(*,1100) XT1, YT1, XT2, YT2,
+ & XB1, YB1, XB2, YB2
+ 1100 FORMAT(/' Top breaks: x,y = ', 2F9.5, 4X, 2F9.5
+ & /' Bot breaks: x,y = ', 2F9.5, 4X, 2F9.5)
+C
+C---- find points adjacent to breaks
+ DO 5 I=1, NB-1
+ IF(SB(I).LE.ST1 .AND. SB(I+1).GT.ST1) IT1 = I+1
+ IF(SB(I).LT.ST2 .AND. SB(I+1).GE.ST2) IT2 = I
+ IF(SB(I).LE.SB1 .AND. SB(I+1).GT.SB1) IB1 = I
+ IF(SB(I).LT.SB2 .AND. SB(I+1).GE.SB2) IB2 = I+1
+ 5 CONTINUE
+C
+ DSAVG = (SB(NB)-SB(1))/FLOAT(NB-1)
+C
+C---- smallest fraction of s increments i+1 and i+2 away from break point
+ SFRAC = 0.33333
+C
+ IF(ATOP .NE. 0.0) THEN
+ ST1P = ST1 + SFRAC*(SB(IT1 )-ST1)
+ ST1Q = ST1 + SFRAC*(SB(IT1+1)-ST1)
+ IF(SB(IT1) .LT. ST1Q) THEN
+C-------- simply move adjacent point to ideal SFRAC location
+ XT1NEW = SEVAL(ST1Q,XB,XBP,SB,NB)
+ YT1NEW = SEVAL(ST1Q,YB,YBP,SB,NB)
+ LT1NEW = .FALSE.
+ ELSE
+C-------- make new point at SFRAC location
+ XT1NEW = SEVAL(ST1P,XB,XBP,SB,NB)
+ YT1NEW = SEVAL(ST1P,YB,YBP,SB,NB)
+ LT1NEW = .TRUE.
+ ENDIF
+C
+ ST2P = ST2 + SFRAC*(SB(IT2 )-ST2)
+ IT2Q = MAX(IT2-1,1)
+ ST2Q = ST2 + SFRAC*(SB(IT2Q)-ST2)
+ IF(SB(IT2) .GT. ST2Q) THEN
+C-------- simply move adjacent point
+ XT2NEW = SEVAL(ST2Q,XB,XBP,SB,NB)
+ YT2NEW = SEVAL(ST2Q,YB,YBP,SB,NB)
+ LT2NEW = .FALSE.
+ ELSE
+C-------- make new point
+ XT2NEW = SEVAL(ST2P,XB,XBP,SB,NB)
+ YT2NEW = SEVAL(ST2P,YB,YBP,SB,NB)
+ LT2NEW = .TRUE.
+ ENDIF
+ ENDIF
+C
+ IF(ABOT .NE. 0.0) THEN
+ SB1P = SB1 + SFRAC*(SB(IB1 )-SB1)
+ SB1Q = SB1 + SFRAC*(SB(IB1-1)-SB1)
+ IF(SB(IB1) .GT. SB1Q) THEN
+C-------- simply move adjacent point
+ XB1NEW = SEVAL(SB1Q,XB,XBP,SB,NB)
+ YB1NEW = SEVAL(SB1Q,YB,YBP,SB,NB)
+ LB1NEW = .FALSE.
+ ELSE
+C-------- make new point
+ XB1NEW = SEVAL(SB1P,XB,XBP,SB,NB)
+ YB1NEW = SEVAL(SB1P,YB,YBP,SB,NB)
+ LB1NEW = .TRUE.
+ ENDIF
+C
+ SB2P = SB2 + SFRAC*(SB(IB2 )-SB2)
+ IB2Q = MIN(IB2+1,NB)
+ SB2Q = SB2 + SFRAC*(SB(IB2Q)-SB2)
+ IF(SB(IB2) .LT. SB2Q) THEN
+C-------- simply move adjacent point
+ XB2NEW = SEVAL(SB2Q,XB,XBP,SB,NB)
+ YB2NEW = SEVAL(SB2Q,YB,YBP,SB,NB)
+ LB2NEW = .FALSE.
+ ELSE
+C-------- make new point
+ XB2NEW = SEVAL(SB2P,XB,XBP,SB,NB)
+ YB2NEW = SEVAL(SB2P,YB,YBP,SB,NB)
+ LB2NEW = .TRUE.
+ ENDIF
+ ENDIF
+C
+cc DSTOP = ABS(SB(IT2)-SB(IT1))
+cc DSBOT = ABS(SB(IB2)-SB(IB1))
+C
+ SIND = SIN(RDEF)
+ COSD = COS(RDEF)
+C
+C---- rotate flap points about the hinge point (XBF,YBF)
+ DO 10 I=1, NB
+ IF(I.GE.IT1 .AND. I.LE.IB1) GO TO 10
+C
+ XBAR = XB(I) - XBF
+ YBAR = YB(I) - YBF
+C
+ XB(I) = XBF + XBAR*COSD + YBAR*SIND
+ YB(I) = YBF - XBAR*SIND + YBAR*COSD
+ 10 CONTINUE
+C
+ IDIF = IT1-IT2-1
+ IF(IDIF.GT.0) THEN
+C----- delete points on upper airfoil surface which "disappeared".
+ NB = NB -IDIF
+ IT1 = IT1-IDIF
+ IB1 = IB1-IDIF
+ IB2 = IB2-IDIF
+ DO 21 I=IT2+1, NB
+ SB(I) = SB(I+IDIF)
+ XB(I) = XB(I+IDIF)
+ YB(I) = YB(I+IDIF)
+ 21 CONTINUE
+ ENDIF
+C
+ IDIF = IB2-IB1-1
+ IF(IDIF.GT.0) THEN
+C----- delete points on lower airfoil surface which "disappeared".
+ NB = NB -IDIF
+ IB2 = IB2-IDIF
+ DO 22 I=IB1+1, NB
+ SB(I) = SB(I+IDIF)
+ XB(I) = XB(I+IDIF)
+ YB(I) = YB(I+IDIF)
+ 22 CONTINUE
+ ENDIF
+C
+C
+ IF(ATOP .EQ. 0.0) THEN
+C
+C------ arc length of newly created surface on top of airfoil
+ DSNEW = ABS(RDEF)*SQRT((XT1-XBF)**2 + (YT1-YBF)**2)
+C
+C------ number of points to be added to define newly created surface
+ NPADD = INT(1.5*DSNEW/DSAVG + 1.0)
+ccc NPADD = INT(1.5*DSNEW/DSTOP + 1.0)
+C
+C------ skip everything if no points are to be added
+ IF(NPADD.EQ.0) GO TO 35
+C
+C------ increase coordinate array length to make room for the new point(s)
+ NB = NB +NPADD
+ IT1 = IT1+NPADD
+ IB1 = IB1+NPADD
+ IB2 = IB2+NPADD
+ DO 30 I=NB, IT1, -1
+ XB(I) = XB(I-NPADD)
+ YB(I) = YB(I-NPADD)
+ 30 CONTINUE
+C
+C------ add new points along the new surface circular arc segment
+ DANG = RDEF / FLOAT(NPADD)
+ XBAR = XT1 - XBF
+ YBAR = YT1 - YBF
+ DO 31 IP=1, NPADD
+ ANG = DANG*(FLOAT(IP) - 0.5)
+ CA = COS(ANG)
+ SA = SIN(ANG)
+C
+ XB(IT1-IP) = XBF + XBAR*CA + YBAR*SA
+ YB(IT1-IP) = YBF - XBAR*SA + YBAR*CA
+ 31 CONTINUE
+C
+ ELSE
+C
+C------ set point in the corner and possibly two adjacent points
+ NPADD = 1
+ IF(LT2NEW) NPADD = NPADD+1
+ IF(LT1NEW) NPADD = NPADD+1
+C
+ NB = NB +NPADD
+ IT1 = IT1+NPADD
+ IB1 = IB1+NPADD
+ IB2 = IB2+NPADD
+ DO 33 I=NB, IT1, -1
+ XB(I) = XB(I-NPADD)
+ YB(I) = YB(I-NPADD)
+ 33 CONTINUE
+C
+ IF(LT1NEW) THEN
+ XB(IT1-1) = XT1NEW
+ YB(IT1-1) = YT1NEW
+ XB(IT1-2) = XT1
+ YB(IT1-2) = YT1
+ ELSE
+ XB(IT1 ) = XT1NEW
+ YB(IT1 ) = YT1NEW
+ XB(IT1-1) = XT1
+ YB(IT1-1) = YT1
+ ENDIF
+C
+ XBAR = XT2NEW - XBF
+ YBAR = YT2NEW - YBF
+ IF(LT2NEW) THEN
+ XB(IT2+1) = XBF + XBAR*COSD + YBAR*SIND
+ YB(IT2+1) = YBF - XBAR*SIND + YBAR*COSD
+ ELSE
+ XB(IT2 ) = XBF + XBAR*COSD + YBAR*SIND
+ YB(IT2 ) = YBF - XBAR*SIND + YBAR*COSD
+ ENDIF
+C
+ ENDIF
+ 35 CONTINUE
+C
+C
+ IF(ABOT .EQ. 0.0) THEN
+C
+C------ arc length of newly created surface on top of airfoil
+ DSNEW = ABS(RDEF)*SQRT((XB1-XBF)**2 + (YB1-YBF)**2)
+C
+C------ number of points to be added to define newly created surface
+ NPADD = INT(1.5*DSNEW/DSAVG + 1.0)
+ccc NPADD = INT(1.5*DSNEW/DSBOT + 1.0)
+C
+C------ skip everything if no points are to be added
+ IF(NPADD.EQ.0) GO TO 45
+C
+C------ increase coordinate array length to make room for the new point(s)
+ NB = NB +NPADD
+ IB2 = IB2+NPADD
+ DO 40 I=NB, IB2, -1
+ XB(I) = XB(I-NPADD)
+ YB(I) = YB(I-NPADD)
+ 40 CONTINUE
+C
+C------ add new points along the new surface circular arc segment
+ DANG = RDEF / FLOAT(NPADD)
+ XBAR = XB1 - XBF
+ YBAR = YB1 - YBF
+ DO 41 IP=1, NPADD
+ ANG = DANG*(FLOAT(IP) - 0.5)
+ CA = COS(ANG)
+ SA = SIN(ANG)
+C
+ XB(IB1+IP) = XBF + XBAR*CA + YBAR*SA
+ YB(IB1+IP) = YBF - XBAR*SA + YBAR*CA
+ 41 CONTINUE
+C
+ ELSE
+
+C------ set point in the corner and possibly two adjacent points
+ NPADD = 1
+ IF(LB2NEW) NPADD = NPADD+1
+ IF(LB1NEW) NPADD = NPADD+1
+C
+ NB = NB +NPADD
+ IB2 = IB2+NPADD
+ DO 43 I=NB, IB2, -1
+ XB(I) = XB(I-NPADD)
+ YB(I) = YB(I-NPADD)
+ 43 CONTINUE
+C
+ IF(LB1NEW) THEN
+ XB(IB1+1) = XB1NEW
+ YB(IB1+1) = YB1NEW
+ XB(IB1+2) = XB1
+ YB(IB1+2) = YB1
+ ELSE
+ XB(IB1 ) = XB1NEW
+ YB(IB1 ) = YB1NEW
+ XB(IB1+1) = XB1
+ YB(IB1+1) = YB1
+ ENDIF
+C
+ XBAR = XB2NEW - XBF
+ YBAR = YB2NEW - YBF
+ IF(LB2NEW) THEN
+ XB(IB2-1) = XBF + XBAR*COSD + YBAR*SIND
+ YB(IB2-1) = YBF - XBAR*SIND + YBAR*COSD
+ ELSE
+ XB(IB2 ) = XBF + XBAR*COSD + YBAR*SIND
+ YB(IB2 ) = YBF - XBAR*SIND + YBAR*COSD
+ ENDIF
+C
+ ENDIF
+ 45 CONTINUE
+C
+ LGSAME = .FALSE.
+C
+C
+C---- check new geometry for splinter segments
+ STOL = 0.2
+ CALL SCHECK(XB,YB,NB, STOL, LCHANGE)
+C
+C---- spline new geometry
+ 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
+ LBFLAP = .TRUE.
+C
+ IF(LGSYM) THEN
+ WRITE(*,*)
+ WRITE(*,*) 'Disabling symmetry enforcement'
+ LGSYM = .FALSE.
+ ENDIF
+C
+C
+ IF(.NOT.LPLOT) THEN
+ CALL PLTINI
+ ENDIF
+C
+C---- save current color and set new color
+ CALL GETCOLOR(ICOL0)
+C
+ CALL NEWCOLORNAME('green')
+ CALL PLOT((XBF-XOFF)*XSF,(YBF-YOFF)*YSF,3)
+ CALL PLOT((XT1-XOFF)*XSF,(YT1-YOFF)*YSF,2)
+ CALL PLOT((XBF-XOFF)*XSF,(YBF-YOFF)*YSF,3)
+ CALL PLOT((XB1-XOFF)*XSF,(YB1-YOFF)*YSF,2)
+C
+ IF(ATOP .EQ. 0.0) THEN
+ XBAR = XT1 - XBF
+ YBAR = YT1 - YBF
+ XT1C = XBF + XBAR*COSD + YBAR*SIND
+ YT1C = YBF - XBAR*SIND + YBAR*COSD
+ CALL PLOT((XBF -XOFF)*XSF,(YBF -YOFF)*YSF,3)
+ CALL PLOT((XT1C-XOFF)*XSF,(YT1C-YOFF)*YSF,2)
+ ENDIF
+C
+ IF(ABOT .EQ. 0.0) THEN
+ XBAR = XB1 - XBF
+ YBAR = YB1 - YBF
+ XB1C = XBF + XBAR*COSD + YBAR*SIND
+ YB1C = YBF - XBAR*SIND + YBAR*COSD
+ CALL PLOT((XBF -XOFF)*XSF,(YBF -YOFF)*YSF,3)
+ CALL PLOT((XB1C-XOFF)*XSF,(YB1C-YOFF)*YSF,2)
+ ENDIF
+C
+ CALL NEWCOLORNAME('red')
+ CALL PLSYMB((XBF-XOFF)*XSF,(YBF-YOFF)*YSF,0.5*SHT,1,0.0,0)
+C
+ CALL PLTAIR(XB,XBP,YB,YBP,SB,NB, XOFF,XSF,YOFF,YSF,'magenta')
+ CALL PLNEWP('magenta')
+C
+ LGEOPL = .FALSE.
+C
+ CALL NEWCOLOR(ICOL0)
+ RETURN
+ END ! FLAP
+
+
+ LOGICAL FUNCTION INSIDE(X,Y,N, XF,YF)
+ DIMENSION X(N),Y(N)
+C-------------------------------------
+C Returns .TRUE. if point XF,YF
+C is inside contour X(i),Y(i).
+C-------------------------------------
+C
+C---- integrate subtended angle around airfoil perimeter
+ ANGLE = 0.0
+ DO 10 I=1, N
+ IP = I+1
+ IF(I.EQ.N) IP = 1
+ XB1 = X(I) - XF
+ YB1 = Y(I) - YF
+ XB2 = X(IP) - XF
+ YB2 = Y(IP) - YF
+ ANGLE = ANGLE + (XB1*YB2 - YB1*XB2)
+ & / SQRT((XB1**2 + YB1**2)*(XB2**2 + YB2**2))
+ 10 CONTINUE
+C
+C---- angle = 0 if XF,YF is outside, angle = +/- 2 pi if XF,YF is inside
+ INSIDE = ABS(ANGLE) .GT. 1.0
+C
+ RETURN
+ END ! INSIDE
+
+
+
+ SUBROUTINE GETXYF(X,XP,Y,YP,S,N, TOPS,BOTS,XF,YF)
+ DIMENSION X(N),XP(N),Y(N),YP(N),S(N)
+C
+ IF(XF .EQ. -999.0)
+ & CALL ASKR('Enter flap hinge x location^',XF)
+C
+C---- find top and bottom y at hinge x location
+ TOPS = S(1) + (X(1) - XF)
+ BOTS = S(N) - (X(N) - XF)
+ CALL SINVRT(TOPS,XF,X,XP,S,N)
+ CALL SINVRT(BOTS,XF,X,XP,S,N)
+ TOPY = SEVAL(TOPS,Y,YP,S,N)
+ BOTY = SEVAL(BOTS,Y,YP,S,N)
+C
+ WRITE(*,1000) TOPY, BOTY
+ 1000 FORMAT(/' Top surface: y =', F8.4,' y/t = 1.0'
+ & /' Bottom surface: y =', F8.4,' y/t = 0.0')
+C
+ IF(YF .EQ. -999.0)
+ & CALL ASKR(
+ & 'Enter flap hinge y location (or 999 to specify y/t)^',YF)
+C
+ IF(YF .EQ. 999.0) THEN
+ CALL ASKR('Enter flap hinge relative y/t location^',YREL)
+ YF = TOPY*YREL + BOTY*(1.0-YREL)
+ ENDIF
+C
+ RETURN
+ END ! GETXYF
+
+
+
+ SUBROUTINE PLOTG
+C--------------------------------------------------------------
+C Plots buffer airfoil with ticked chord line or grid
+C--------------------------------------------------------------
+ INCLUDE 'XFOIL.INC'
+C
+ DATA LMASK1, LMASK2, LMASK3 / -32640, -30584, -21846 /
+ INCLUDE 'XDES.INC'
+C
+C---- node tick mark size and corner symbol size
+ DTICK = GTICK*(SB(NB)-SB(1))
+ SSH = DTICK * 3.0
+C
+ CALL NCALC(XB,YB,SB,NB,W1,W2)
+C
+ IF(LGGRID) THEN
+ CALL GRDAIR(XGMIN,XGMAX,YGMIN,YGMAX,DXYG,DXYG,CHG,.TRUE.,.TRUE.,
+ & XOFF,XSF,YOFF,YSF, LMASK2)
+ XL0 = XMOD(XGMIN)
+ YL0 = YMOD(YGMAX) + 2.0*CH
+ ELSE
+C------ plot chord line and tick marks every 10% chord
+ CALL NEWPEN(1)
+ CALL PLOT(XMOD(0.0),YMOD(0.0),3)
+ CALL PLOT(XMOD(1.0),YMOD(0.0),2)
+ DO ITICK=1, 10
+ XPLT = FLOAT(ITICK)/10.0
+ CALL PLOT(XMOD(XPLT),YMOD(0.003),3)
+ CALL PLOT(XMOD(XPLT),YMOD(-.003),2)
+ ENDDO
+C
+ XL0 = XMOD(XBMIN)
+ YL0 = YMOD(YBMAX) + 2.0*CH
+ ENDIF
+ IF(LPLCAM) YL0 = YSF*(YCMAX-DYOFFC-YOFF) + 2.0*CH
+C
+ CALL PLFLUSH
+C
+ CALL NEWPEN(2)
+ CALL PLTAIR(XB,XBP,YB,YBP,SB,NB, XOFF,XSF,YOFF,YSF,'black')
+C
+ IF(LGTICK) THEN
+C----- draw tiny tick mark normal to airfoil surface at each panel node
+ DO I=2, NB-1
+ CALL PLOT(XMOD(XB(I) ),YMOD(YB(I) ),3)
+ CALL PLOT(XMOD(XB(I)-DTICK*W1(I)),YMOD(YB(I)-DTICK*W2(I)),2)
+ ENDDO
+ ENDIF
+c
+cC---- plot symbol at nose
+c CALL NSFIND(STLE,XB,XBP,YB,YBP,SB,NB)
+c XT = SEVAL(STLE,XB,XBP,SB,NB)
+c YT = SEVAL(STLE,YB,YBP,SB,NB)
+c CALL PLSYMB(XMOD(XT),YMOD(YT),0.005*XSF,5,0.0,0)
+c
+C---- put symbol at any doubled point
+ DO I=1, NB-1
+ IF(SB(I) .EQ. SB(I+1))
+ & CALL PLSYMB(XMOD(XB(I)),YMOD(YB(I)),SSH,5,0.0,0)
+ ENDDO
+C
+ IF(LPLCAM) THEN
+ CALL PLTCAM(' ')
+ ENDIF
+C
+ IF(LGPARM) THEN
+ CALL NEWPEN(3)
+ CALL GPARPL(XL0,YL0,0.7*CH,.TRUE.,NAME,
+ & CHORDB,AREAB,RADBLE,ANGBTE,
+ & EI11BA,EI22BA,APX1BA,APX2BA,
+ & EI11BT,EI22BT,APX1BT,APX2BT,
+ & THICKB,CAMBRB)
+ ENDIF
+C
+ CALL PLFLUSH
+C
+ LGEOPL = .TRUE.
+ NOVER = 0
+C
+ RETURN
+ END ! PLOTG
+
+
+
+ SUBROUTINE PLTCAM(COLIN)
+C--------------------------------------------
+C Plots camber & thickness distributions
+C--------------------------------------------
+ INCLUDE 'XFOIL.INC'
+ CHARACTER*(*) COLIN
+ CHARACTER*32 COLC, COLT
+ DATA LMASK1, LMASK2, LMASK3 / -32640, -30584, -21846 /
+C
+C---- plot camber/thickness only if camber/tickness plot is being shown
+ IF(.NOT.LPLCAM) RETURN
+C
+ CALL NEWPEN(1)
+ CALL GRDAIR(XGMIN,XGMAX,YCMIN,YCMAX,DXYG,DXYG,CHG,.FALSE.,.TRUE.,
+ & XOFF,XSF,DYOFFC+YOFF,YSF, LMASK2)
+C
+ CALL GETCAM(XCM,YCM,NCM,XTK,YTK,NTK,
+ & XB,XBP,YB,YBP,SB,NB )
+ CALL SCALC(XCM,YCM,SCM,NCM)
+ CALL SEGSPL(XCM,XCMP,SCM,NCM)
+ CALL SEGSPL(YCM,YCMP,SCM,NCM)
+ CALL SCALC(XTK,YTK,STK,NTK)
+ CALL SEGSPL(XTK,XTKP,STK,NTK)
+ CALL SEGSPL(YTK,YTKP,STK,NTK)
+C
+ IF(COLIN(1:1) .EQ. ' ') THEN
+ COLC = 'green'
+ COLT = 'cyan'
+ ELSE
+ COLC = COLIN
+ COLT = COLIN
+ ENDIF
+C
+ CALL NEWPEN(2)
+ YOF = YOFF + DYOFFC
+ CALL PLTAIR(XTK,XTKP,YTK,YTKP,STK,NTK,XOFF,XSF, YOF, YSF,COLT)
+ CALL PLTAIR(XTK,XTKP,YTK,YTKP,STK,NTK,XOFF,XSF,-YOF,-YSF,COLT)
+C--- Offset for camber includes offset for LE camber point
+ YOFF1C = YOFF + DYOFFC + YCM(1)
+ CALL PLTAIR(XCM,XCMP,YCM,YCMP,SCM,NCM,XOFF,XSF, YOFF1C,YSF,COLC)
+C
+ RETURN
+ END ! PLTCAM
+
+
+ SUBROUTINE PLNEWP(COLOR)
+ INCLUDE 'XFOIL.INC'
+ CHARACTER*(*) COLOR
+C
+ INCLUDE 'XDES.INC'
+C
+C---- don't plot geometric parameters if camber/tickness plot is being shown
+ IF(LPLCAM) RETURN
+C
+ CALL GETCOLOR(ICOL0)
+C
+ CALL NEWCOLORNAME(COLOR)
+ CALL NEWPEN(3)
+C
+ NOVER = NOVER + 1
+ IF(LGGRID) THEN
+ XL0 = XMOD(XGMIN) + 2.0*CH + 9.0*CH*FLOAT(NOVER)
+ YL0 = YMOD(YGMAX) + 2.0*CH
+ ELSE
+ XL0 = XMOD(XBMIN) + 2.0*CH + 9.0*CH*FLOAT(NOVER)
+ YL0 = YMOD(YBMAX) + 2.0*CH
+ ENDIF
+
+ IF(LPLCAM) YL0 = YSF*(YCMAX-YOFF-DYOFFC) + 2.0*CH
+C
+ IF(LGPARM) THEN
+ CALL GPARPL(XL0,YL0,0.7*CH,.FALSE.,NAME,
+ & CHORDB,AREAB,RADBLE,ANGBTE,
+ & EI11BA,EI22BA,APX1BA,APX2BA,
+ & EI11BT,EI22BT,APX1BT,APX2BT,
+ & THICKB,CAMBRB)
+ ENDIF
+C
+ CALL NEWCOLOR(ICOL0)
+ CALL PLFLUSH
+C
+ RETURN
+ END ! PLNEWP
+
+
+
+ SUBROUTINE GPARPL(X0,Y0,CH, LABEL, NAME,
+ & CHORD,AREA,RADLE,ANGTE,
+ & EI11A,EI22A,APX1A,APX2A,
+ & EI11T,EI22T,APX1T,APX2T,
+ & THICK,CAMBR)
+ LOGICAL LABEL
+ EXTERNAL PLCHAR
+ CHARACTER NAME*(*)
+C
+ RTD = 45.0/ATAN(1.0)
+C
+ XSPACE = 30.0*CH
+ YSPACE = 2.0*CH
+C
+ X = X0
+ Y = Y0
+C
+ IF(LABEL) THEN
+ CALL PLCHAR(X,Y,CH,' = ',0.0, 9)
+ CALL PLMATH(X,Y,CH,' Oq ',0.0, 9)
+ CALL PLSUBS(X+3.0*CH,Y,CH,'TE',0.0, 2, PLCHAR)
+ ENDIF
+ CALL PLNUMB(X+9.0*CH,Y,CH,ANGTE*RTD ,0.0, 2)
+ CALL PLMATH(999.,Y,CH,'"' ,0.0, 1)
+ Y = Y + YSPACE
+C
+ IF(LABEL) THEN
+ CALL PLCHAR(X,Y,CH,' r = ',0.0, 9)
+ CALL PLSUBS(X+3.0*CH,Y,CH,'LE',0.0, 2, PLCHAR)
+ ENDIF
+ CALL PLNUMB(X+9.0*CH,Y,CH,RADLE,0.0, 5)
+ Y = Y + YSPACE
+C
+ IF(LABEL) THEN
+ CALL PLCHAR(X,Y,CH,'camber = ',0.0, 9)
+ ENDIF
+ CALL PLNUMB(X+9.0*CH,Y,CH,CAMBR,0.0, 5)
+ Y = Y + YSPACE
+C
+ IF(LABEL) THEN
+ CALL PLCHAR(X,Y,CH,'thick. = ',0.0, 9)
+ ENDIF
+ CALL PLNUMB(X+9.0*CH,Y,CH,THICK,0.0, 5)
+ Y = Y + YSPACE
+C
+ IF(LABEL) THEN
+ CALL PLCHAR(X,Y,CH,' area = ',0.0, 9)
+ ENDIF
+ CALL PLNUMB(X+9.0*CH,Y,CH, AREA,0.0, 5)
+ Y = Y + YSPACE
+C
+C
+c X = X0 + XSPACE
+c Y = Y0
+cC
+c Y = Y + YSPACE
+cC
+c CALL PLMATH(X,Y,1.4*CH,'I',0.0,1)
+c CALL PLMATH(X,Y,CH,' 2 ',0.0,-1)
+c CALL PLCHAR(X,Y,CH,' (y-y ) ds = ',0.0,-1)
+c CALL PLNUMB(999.,Y,CH, 1000.0*EI11T,0.0,4)
+c CALL PLMATH(999.,Y,CH,'#' ,0.0,1)
+c CALL PLCHAR(999.,Y,CH, '10' ,0.0,2)
+c CALL PLMATH(999.,Y,CH, '3',0.0,1)
+c CALL PLSUBS(X+4.0*CH,Y,CH,'o',0.0,1,PLCHAR)
+c Y = Y + YSPACE
+cC
+c CALL PLMATH(X,Y,1.4*CH,'I',0.0,1)
+c CALL PLMATH(X,Y,CH,' 2 ',0.0,-1)
+c CALL PLCHAR(X,Y,CH,' (y-y ) dA = ',0.0,-1)
+c CALL PLNUMB(999.,Y,CH, 1000.0*EI11A,0.0,4)
+c CALL PLMATH(999.,Y,CH,'#' ,0.0,1)
+c CALL PLCHAR(999.,Y,CH, '10' ,0.0,2)
+c CALL PLMATH(999.,Y,CH, '3',0.0,1)
+c CALL PLSUBS(X+4.0*CH,Y,CH,'o',0.0,1,PLCHAR)
+c Y = Y + YSPACE
+cC
+c CALL PLMATH(X,Y,CH,' ',0.0,-1)
+c CALL PLCHAR(X,Y,CH,' area = ',0.0,-1)
+c CALL PLNUMB(999.,Y,CH, AREA,0.0, 5)
+c Y = Y + YSPACE
+C
+C--- Plot airfoil name over data list
+ CALL PLCHAR(X+9.0*CH,Y,CH,NAME,0.0, 12)
+C
+ RETURN
+ END ! GPARPL
+
+
+
+ SUBROUTINE GRDAIR(XGMIN,XGMAX, YGMIN,YGMAX,DXGN,DYGN,CHG,
+ & LXAXIS,LYAXIS,
+ & XOFF,XSF,YOFF,YSF, LMASK)
+ LOGICAL LXAXIS,LYAXIS
+C----------------------------------------
+C Plots grid with axes.
+C Intended for airfoil plot.
+C----------------------------------------
+ INCLUDE 'XDES.INC'
+C
+ CALL NEWPEN(1)
+C
+C---- plot outline
+ CALL PLOT(XMOD(XGMIN),YMOD(YGMIN),3)
+ CALL PLOT(XMOD(XGMAX),YMOD(YGMIN),2)
+ CALL PLOT(XMOD(XGMAX),YMOD(YGMAX),2)
+ CALL PLOT(XMOD(XGMIN),YMOD(YGMAX),2)
+ CALL PLOT(XMOD(XGMIN),YMOD(YGMIN),2)
+C
+ IF(LXAXIS)
+ & CALL XAXIS(XMOD(XGMIN),YMOD(YGMIN),(XGMAX-XGMIN)*XSF,
+ & DXGN*XSF, XGMIN,DXGN,CHG,-2)
+ IF(LYAXIS)
+ & CALL YAXIS(XMOD(XGMIN),YMOD(YGMIN),(YGMAX-YGMIN)*YSF,
+ & DYGN*YSF, YGMIN,DYGN,CHG,-2)
+C
+C---- fine grid
+ NXG = INT((XGMAX-XGMIN)/DXGN + 0.1)
+ NYG = INT((YGMAX-YGMIN)/DYGN + 0.1)
+ NXG = MAX(1,NXG)
+ NYG = MAX(1,NYG)
+C
+ X0 = XMOD(XGMIN)
+ Y0 = YMOD(YGMIN)
+ DXG = (XMOD(XGMAX)-X0)/NXG
+ DYG = (YMOD(YGMAX)-Y0)/NYG
+ CALL PLGRID(X0,Y0,NXG,DXG,NYG,DYG, LMASK)
+C
+ RETURN
+ END ! GRDAIR
+
+
+
+ SUBROUTINE PLTAIR(XX,XXP,YY,YYP,SS,NN, XOFF,XSF,YOFF,YSF,COLOR)
+ DIMENSION XX(NN), XXP(NN), YY(NN), YYP(NN), SS(NN)
+ CHARACTER*(*) COLOR
+C-----------------------------
+C Plots passed-in airfoil
+C-----------------------------
+ LOGICAL LCOLOR
+ XMOD(XTMP) = XSF * (XTMP - XOFF)
+ YMOD(YTMP) = YSF * (YTMP - YOFF)
+C
+ NT = 20
+ccc NT = 50
+C
+ LCOLOR = COLOR(1:1) .NE. ' '
+C
+ IF(LCOLOR) THEN
+ CALL GETCOLOR(ICOL0)
+ CALL NEWCOLORNAME(COLOR)
+ ENDIF
+C
+ DO 60 I=2, NN
+ DS = SS(I) - SS(I-1)
+ CALL PLOT(XMOD(XX(I-1)),YMOD(YY(I-1)),3)
+C
+C------ subdivide current panel into NT segments for smoother airfoil plot
+ DO 610 IT=1, NT
+ ST = SS(I-1) + DS*FLOAT(IT)/FLOAT(NT)
+ XT = SEVAL(ST,XX,XXP,SS,NN)
+ YT = SEVAL(ST,YY,YYP,SS,NN)
+ CALL PLOT(XMOD(XT),YMOD(YT),2)
+ 610 CONTINUE
+ 60 CONTINUE
+C
+ IF(LCOLOR) CALL NEWCOLOR(ICOL0)
+C
+ CALL PLFLUSH
+C
+ RETURN
+ END ! PLTAIR
+
+
+
+ SUBROUTINE OVER(FNAME1)
+C----------------------------------------------------
+C Overlays plot of airfoil from coordinate file.
+C----------------------------------------------------
+ INCLUDE 'XFOIL.INC'
+ CHARACTER*(*) FNAME1
+C
+ CHARACTER*32 NAME0, NAMEW
+ CHARACTER*80 ISPARS0
+C
+ IF(FNAME1(1:1).NE.' ') THEN
+ FNAME = FNAME1
+ ELSE
+C----- no argument... get it somehow
+ IF(ONAME(1:1).NE.' ') THEN
+C------ offer existing default
+ WRITE(*,1100) ONAME
+ 1100 FORMAT(/' Enter filename: ', A)
+ READ(*,1000) FNAME
+ 1000 FORMAT(A)
+ CALL STRIP(FNAME,NFN)
+ IF(NFN.EQ.0) FNAME = ONAME
+ ELSE
+C------ just ask for filename
+ CALL ASKS('Enter filename^',FNAME)
+ ENDIF
+ ENDIF
+C
+ LU = 9
+ CALL AREAD(LU,FNAME,2*IQX,W1,W2,NN,NAME0,ISPARS0,ITYPE,1)
+ IF(ITYPE.EQ.0) RETURN
+C
+C---- set new default filename
+ ONAME = FNAME
+C
+ IF(LNORM) THEN
+C----- normalize to unit chord
+ CALL NORM(W1,W3,W2,W4,W5,NN)
+ ELSE
+ CALL SCALC(W1,W2,W5,NN)
+ CALL SEGSPL(W1,W3,W5,NN)
+ CALL SEGSPL(W2,W4,W5,NN)
+ ENDIF
+C
+ NAMEW = NAME
+ SWLE = SBLE
+ CHORDW = CHORDB
+ AREAW = AREAB
+ RADWLE = RADBLE
+ ANGWTE = ANGBTE
+ EI11WA = EI11BA
+ EI22WA = EI22BA
+ APX1WA = APX1BA
+ APX2WA = APX2BA
+ EI11WT = EI11BT
+ EI22WT = EI22BT
+ APX1WT = APX1BT
+ APX2WT = APX2BT
+ THICKW = THICKB
+ CAMBRW = CAMBRB
+C
+ NAME = NAME0
+ CALL GEOPAR(W1,W3,W2,W4,W5,NN,W6,
+ & SBLE,CHORDB,AREAB,RADBLE,ANGBTE,
+ & EI11BA,EI22BA,APX1BA,APX2BA,
+ & EI11BT,EI22BT,APX1BT,APX2BT,
+ & THICKB,CAMBRB )
+C
+ IF(.NOT.LPLOT) THEN
+ CALL PLTINI
+ccc CALL PLOT(0.05,0.30,-3)
+ ENDIF
+C
+ CALL NEWPEN(2)
+ CALL PLTAIR(W1,W3,W2,W4,W5,NN, XOFF,XSF, YOFF,YSF,'cyan')
+ CALL PLNEWP('cyan')
+C
+C---- restore parameters
+ NAME = NAMEW
+ SBLE = SWLE
+ CHORDB = CHORDW
+ AREAB = AREAW
+ RADBLE = RADWLE
+ ANGBTE = ANGWTE
+ EI11BA = EI11WA
+ EI22BA = EI22WA
+ APX1BA = APX1WA
+ APX2BA = APX2WA
+ EI11BT = EI11WT
+ EI22BT = EI22WT
+ APX1BT = APX1WT
+ APX2BT = APX2WT
+ THICKB = THICKW
+ CAMBRB = CAMBRW
+C
+ RETURN
+ END ! OVER
+
+