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 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( & /' 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