*********************************************************************** C Module: xtcam.f C C Copyright (C) 2000 Harold Youngren, 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*********************************************************************** SUBROUTINE CAMB C------------------------------------------- C Camber modification routine. C------------------------------------------- INCLUDE 'XFOIL.INC' C CHARACTER*72 LINE CHARACTER*4 COMAND, COMOLD CHARACTER*128 COMARG, ARGOLD CHARACTER*1 ANS C REAL XBOX(2), YBOX(2) DIMENSION IINPUT(20) DIMENSION RINPUT(20) LOGICAL ERROR, LRECALC, LCLEAR, LGPARSAVE C EXTERNAL NEWPLOTC C DATA LMASK0, LMASK1, LMASK2, LMASK3 / -1, -32640, -30584, -21846 / C 1000 FORMAT(A) C LGPARSAVE = LGPARM COMAND = '****' COMARG = ' ' LRECALC = .FALSE. LCLEAR = .TRUE. LU = 8 C COMOLD = COMAND ARGOLD = COMARG C IF(.NOT.LPLCAM) THEN WRITE(*,*) 'Enabling camber,thickness plotting' LPLCAM = .TRUE. CALL GOFINI ENDIF C C--- Check chordline direction (should be unrotated for camber routines C to function correctly 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)) AROT = ATAN2(YLE-YTE,XTE-XLE) / DTOR IF(ABS(AROT).GT.1.0) THEN WRITE(*,*) ' ' WRITE(*,*) 'Warning: CAMB does not work well on rotated foils' WRITE(*,*) 'Current chordline angle: ',AROT WRITE(*,*) 'Proceeding anyway...' ENDIF C CHS = 0.5*CHG LDCPLOT = .FALSE. LGPARM = .NOT.LDCPLOT C WRITE(*,1200) C C-------------------------------------------------------------- C---- pick up here to initialize camber and loading 100 CONTINUE C C---- find leftmost point cc CALL LEFIND(SBL,XB,XBP,YB,YBP,SB,NB) CALL XLFIND(SBL,XB,XBP,YB,YBP,SB,NB) C XBL = SEVAL(SBL, XB,XBP,SB,NB) YBL = SEVAL(SBL, YB,YBP,SB,NB) XBR = 0.5*(XB(1)+XB(NB)) YBR = 0.5*(YB(1)+YB(NB)) C C---- set "chordline" axis vector for camber,thickness definitions XBCH = XBR - XBL cc YBCH = YBR - YBL YBCH = 0. SBCH = SQRT(XBCH**2 + YBCH**2) C C---- find the current buffer airfoil camber and thickness CALL GETCAM(XCM,YCM,NCM,XTK,YTK,NTK, & XB,XBP,YB,YBP,SB,NB ) C write(*,*) 'xc0 xc1', xcm(1) , xcm(ncm) NCAM = MIN( 201 , NTX ) DO K=1, NCAM XCAM(K) = XCM(1) + (XCM(NCM)-XCM(1))*FLOAT(K-1)/FLOAT(NCAM-1) ENDDO C IF(LCLEAR) THEN C---- initialize added camber to zero NCADD = 2 XCADD(1) = XCAM(1) XCADD(2) = XCAM(NCAM) YCADD(1) = 0.0 YCADD(2) = 0.0 C---- initialize added loading to zero NPADD = 2 XPADD(1) = XCAM(1) XPADD(2) = XCAM(NCAM) YPADD(1) = 0.0 YPADD(2) = 0.0 C---- spline added camber line y(x) and added loading dCp(x) CALL SEGSPL(YCADD,YCADDP,XCADD,NCADD) CALL SEGSPL(YPADD,YPADDP,XPADD,NPADD) C----- interpolate to dense plotting array DO K=1, NCAM YCAM(K) = SEVAL(XCAM(K),YCADD,YCADDP,XCADD,NCADD) YCAMP(K) = DEVAL(XCAM(K),YCADD,YCADDP,XCADD,NCADD) PCAM(K) = SEVAL(XCAM(K),YPADD,YPADDP,XPADD,NPADD) PCAMP(K) = DEVAL(XCAM(K),YPADD,YPADDP,XPADD,NPADD) ENDDO LCLEAR = .FALSE. ENDIF C C-------------------------------------------------------------- C---- pick up here to find and display current camber and added camber line properties 200 CONTINUE C WRITE(*,*) WRITE(*,*) 'Buffer airfoil thickness and camber:' CALL TCBUF C XMX = 0.0 YMX = 0.0 DO K=1, NCAM IF(ABS(YCAM(K)) .GT. ABS(YMX)) THEN XMX = XCAM(K) YMX = YCAM(K) ENDIF END DO CHRD = XCAM(NCAM) - XCAM(1) ALE = ATAN( DEVAL(XCAM(1) ,YCAM,YCAMP,XCAM,NCAM) ) / DTOR ATE = ATAN( DEVAL(XCAM(NCAM),YCAM,YCAMP,XCAM,NCAM) ) / DTOR WRITE(*,1100) ALE, ATE, YMX/CHRD, XMX/CHRD 1100 FORMAT(/' Added camber line incidence at LE = ', F6.2, ' deg.', & /' Added camber line incidence at TE = ', F6.2, ' deg.', & /' Max added camber y/c = ', F8.4, ' at x/c = ', F7.3 ) C C-------------------------------------------------------------- C---- pick up here to replot everything 300 CONTINUE LGPARM = .NOT.LDCPLOT CALL PLTINI CALL PLOTG CALL PLOTC C C================================================== C---- top of menu loop 500 CALL ASKC('..CAMB^',COMAND,COMARG) C C---- process previous command ? IF(COMAND(1:1).EQ.'!') THEN IF(COMOLD.EQ.'****') THEN WRITE(*,*) 'Previous ..CAMB command not valid' GO TO 500 ELSE COMAND = COMOLD COMARG = ARGOLD ENDIF ENDIF C IF(COMAND.EQ.' ') THEN C----- just was typed... clean up plotting and exit CAMP IF(LPLOT) CALL PLEND LPLOT = .FALSE. CALL CLRZOOM LGPARM = LGPARSAVE 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 IF(COMAND.EQ.' ') THEN IF(LPLOT) CALL PLEND RETURN C ELSEIF(COMAND.EQ.'? ') THEN WRITE(*,1200) 1200 FORMAT( & /' Return to GDES' & /' TFAC rr Scale existing thickness and camber' & /' TSET rr Set new thickness and camber' & /' HIGH rr Move camber and thickness highpoints' & /' WRTC Write airfoil camber x/c,y/c to file' & //' RDAC Read added camber x/c,y/c from file' & /' SETC Set added camber x/c,y/c from camberline' & /' INPC Input added camber x/c,y/c from keyboard' & /' MODC Modify added camber x/c,y/c with cursor' & /' INPP Input added loading x/c,DCp from keyboard' & /' MODP Modify added loading x/c,DCp with cursor' & /' SLOP Toggle modified-camber,dCp slope matching flag' & /' SCAL r Scale the added camber' & /' CLR Clear the added camber' & /' ADD Add added camber to the existing camberline' & //' DCPL Toggle DCp plot' & /' CPLI rr Change DCp axis plot limits' & //' Blow Blowup plot region' & /' Rese Reset plot scale and origin' & //' SIZE r Change absolute plot-object size' & /' .ANNO Annotate plot' & /' HARD Hardcopy current plot') 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.'TFAC') THEN CALL TCSCAL(RINPUT,NINPUT) GO TO 100 C C-------------------------------------------------------- ELSEIF(COMAND.EQ.'TSET') THEN CALL TCSET(RINPUT,NINPUT) GO TO 100 C C-------------------------------------------------------- ELSEIF(COMAND.EQ.'HIGH') THEN CALL HIPNT(RINPUT,NINPUT) GO TO 100 C C-------------------------------------------------------- ELSEIF(COMAND.EQ.'WRTC') THEN CALL ASKS('Enter output camber filename^',FNAME) C OPEN(LU,FILE=FNAME,STATUS='OLD',ERR=12) WRITE(*,*) WRITE(*,*) 'Output file exists. Overwrite? Y' READ(*,1000) ANS IF(INDEX('Nn',ANS).EQ.0) GO TO 13 C CLOSE(LU) WRITE(*,*) 'Current camber not saved.' GO TO 500 C 12 OPEN(LU,FILE=FNAME,STATUS='NEW',ERR=15) 13 REWIND(LU) C C--- Write out normalized camber coordinates (x/c in range 0->1, y/c) WRITE(LU,1000) 'Camber: '//NAME DO K = 1, NCM WRITE(LU,14) (XCM(K)-XCM(1))/XBCH,(YCM(K)-YCM(1))/XBCH END DO CLOSE(LU) GO TO 500 C 14 FORMAT(2(1X,F12.6)) C 15 WRITE(*,*) 'Error opening camber save file' GO TO 500 C C-------------------------------------------------------- ELSEIF(COMAND.EQ.'RDAC ') THEN CALL ASKS('Enter added camber filename^',FNAME) OPEN(LU,FILE=FNAME,STATUS='OLD',ERR=19) READ(LU,1000,ERR=18,END=18) LINE NCADD = 0 DO K = 1, NTX READ(LU,*,ERR=18,END=18) XX,YY NCADD = NCADD + 1 XCADD(NCADD) = XX YCADD(NCADD) = YY END DO 18 CLOSE(LU) IF(NCADD.LE.1 .OR. (XCADD(NCADD)-XCADD(1)).EQ.0.0) THEN NCADD = 2 XCADD(1) = XCAM(1) XCADD(2) = XCAM(NCAM) YCADD(1) = 0.0 YCADD(2) = 0.0 WRITE(*,*) 'No added camber points found' GO TO 100 ENDIF C----- normalize input camber to x/c range 0->1 XCORG = XCADD(1) XCSCL = XCADD(NCADD) - XCORG DO K=1, NCADD XCADD(K) = (XCADD(K)-XCORG) / XCSCL YCADD(K) = YCADD(K) / XCSCL ENDDO C----- reorigin and scale added camber to camber line coordinates DO K=1, NCADD XCADD(K) = XCAM(1) + XCADD(K)*XBCH - YCADD(K)*YBCH YCADD(K) = XCADD(K)*YBCH + YCADD(K)*XBCH ENDDO C----- spline camber line y(x) CALL SEGSPL(YCADD,YCADDP,XCADD,NCADD) C----- interpolate to dense plotting array DO K=1, NCAM YCAM(K) = SEVAL(XCAM(K),YCADD,YCADDP,XCADD,NCADD) YCAMP(K) = DEVAL(XCAM(K),YCADD,YCADDP,XCADD,NCADD) ENDDO LDCPLOT = .FALSE. GO TO 200 19 WRITE(*,*) WRITE(*,*) 'Error opening added camber file' GO TO 500 C C-------------------------------------------------------- ELSEIF(COMAND.EQ.'SETC') THEN C----- Set added camber from camberline NCADD = NCM DO K=1, NCM XCADD(K) = XCM(K) YCADD(K) = YCM(K) END DO C----- spline added camber line y(x) CALL SEGSPL(YCADD,YCADDP,XCADD,NCADD) C C----- interpolate to dense plotting array DO K=1, NCAM YCAM(K) = SEVAL(XCAM(K),YCADD,YCADDP,XCADD,NCADD) YCAMP(K) = DEVAL(XCAM(K),YCADD,YCADDP,XCADD,NCADD) ENDDO LDCPLOT = .FALSE. GO TO 200 C C-------------------------------------------------------- ELSEIF(COMAND.EQ.'INPC') THEN C----- Manual input of camber points 20 WRITE(*,2000) 2000 FORMAT(/' Manual input of camber x/c,y/c:', & //' Input x/c, y/c pairs from x/c = 0 to x/c = 1', cc & /' Identical successive points enable a slope break', & /' ends input') C C--- Points of x/c, y/c are added to existing definition of added camber line CALL GETCOLOR(ICOL0) CALL NEWCOLORNAME('RED') NCADD = 0 DO 25 I=1, 2*IQX 23 READ(*,1000,ERR=24) LINE IF(LINE.EQ.' ') GO TO 26 READ(LINE,*,ERR=24,END=24) XX,YY IF(XX.LE.0.0) THEN XX = 0.0 ELSEIF(XX.GE.1.0) THEN XX = 1.0 ENDIF NCADD = NCADD + 1 XCADD(NCADD) = XCAM(1) + XX*XBCH - YY*YBCH YCADD(NCADD) = XX*YBCH + YY*XBCH C XPL = XSF*(XCADD(NCADD)-XOFF) YPL = YSF*(YCADD(NCADD)-YOFF-DYOFFC) CALL PLSYMB(XPL,YPL,CHS*XSF,1,0.0,I-1) CALL PLFLUSH GO TO 25 24 WRITE(*,*) 'try again' GO TO 23 25 CONTINUE C----- Sort points allowing duplicates for slope breaks 26 CALL SORTDUP(NCADD,XCADD,YCADD) CALL FIXDUP (NCADD,XCADD,YCADD) CALL NEWCOLOR(ICOL0) C----- spline camber line y(x) CALL SEGSPL(YCADD,YCADDP,XCADD,NCADD) C C----- interpolate to dense plotting array DO K=1, NCAM YCAM(K) = SEVAL(XCAM(K),YCADD,YCADDP,XCADD,NCADD) YCAMP(K) = DEVAL(XCAM(K),YCADD,YCADDP,XCADD,NCADD) ENDDO LDCPLOT = .FALSE. GO TO 200 C C-------------------------------------------------------- ELSEIF(COMAND.EQ.'INPP') THEN C----- Manual input of loading points 30 WRITE(*,3000) 3000 FORMAT(/' Manual input of loading x/c, DCp:', & //' Input x/c, DCp pairs from x/c = 0 to x/c = 1', cc & /' Identical successive points enable a slope break', & /' ends input') C CALL GETPEN(IPN) CALL GETCOLOR(ICOL0) C CALL NEWPEN(1) CHL = 1.5*CHG YOFFP = (DYOFFP+YOFF)/YSFP CALL GRDAIR(XPMIN,XPMAX,YPMIN,YPMAX,DXYG,DXYP,CHG,.FALSE.,.TRUE., & XOFF,XSF,YOFFP,YSF*YSFP, LMASK2) CALL NEWCOLORNAME('RED') CALL NEWPEN(2) XLAB = (XPMIN -XOFF )*XSF - 4.0*CHL YLAB = (YPMAX-0.5*DXYP-YOFFP)*YSFP*YSF - 0.6*CHL CALL PLCHAR(XLAB,YLAB,CHL,' Cp',0.0,3) CALL PLMATH(XLAB,YLAB,CHL,'O ',0.0,3) C C--- Points of x/c, dCp are added to existing definition of loading line DO 35 I=1, 2*IQX 33 READ(*,1000,ERR=34) LINE IF(LINE.EQ.' ') GO TO 36 READ(LINE,*,ERR=34) XX,YY IF(XX.LE.0.0) THEN XX = 0.0 ELSEIF(XX.GE.1.0) THEN XX = 1.0 ENDIF NPADD = NPADD + 1 XPADD(NPADD) = XCAM(1) + XX*XBCH YPADD(NPADD) = YY C YOFFP = (DYOFFP*YOFF)/YSFP XPL = (XPADD(NPADD)-XOFF )*XSF YPL = (YPADD(NPADD)-YOFFP)*YSFP*YSF CALL PLSYMB(XPL,YPL,CHS,1,0.0,I-1) CALL PLFLUSH GO TO 35 34 WRITE(*,*) 'try again' GO TO 33 35 CONTINUE C----- Sort points allowing duplicates for slope breaks 36 CONTINUE CALL SORTDUP(NPADD,XPADD,YPADD) CALL FIXDUP (NPADD,XPADD,YPADD) C CALL NEWCOLOR(ICOL0) CALL NEWPEN(IPN) C C----- spline loading DCp(x) CALL SEGSPL(YPADD,YPADDP,XPADD,NPADD) C C----- interpolate to dense plotting array DO K=1, NCAM PCAM(K) = SEVAL(XCAM(K),YPADD,YPADDP,XPADD,NPADD) PCAMP(K) = DEVAL(XCAM(K),YPADD,YPADDP,XPADD,NPADD) ENDDO C C----- calculate camber line corresponding to specified loading CALL CPCAM(NCAM,XCAM,YCAM,YCAMP,PCAM,PCAMP) C C----- calculate added lift and moment from added loading CLX = 0.0 CMX = 0.0 DO K=1, NCAM-1 DX = XCAM(K+1) - XCAM(K) XA = 0.5*(XCAM(K+1) + XCAM(K)) PA = 0.5*(PCAM(K+1) + PCAM(K)) CLX = CLX + PA*DX CMX = CMX + PA*DX*(XCMREF-XA) END DO WRITE(*,1110) CLX, CMX C LDCPLOT = .TRUE. GO TO 200 C C-------------------------------------------------------- ELSEIF(COMAND.EQ.'MODC') THEN C----- Interactively modify camber XBOX(1) = XMARG XBOX(2) = XPAGE-XMARG YBOX(1) = YMARG YBOX(2) = YPAGE-YMARG XOFF1 = XOFF YOFF1 = YOFF+DYOFFC XSF1 = XSF YSF1 = YSF CALL MODIFY(NTX,1,NCAM,1,1, & XCAM,YCAM,YCAMP, LCSLOP, & K1,K2,ISMOD,IFMOD, & XBOX,YBOX, XBOX,YBOX,SIZE, & XOFF1,YOFF1,XSF1,YSF1, 'RED',' ', & NEWPLOTC) LDCPLOT = .FALSE. GO TO 200 C C-------------------------------------------------------- ELSEIF(COMAND.EQ.'MODP') THEN C----- Interactively modify loading IF(.NOT.LDCPLOT) THEN LDCPLOT = .TRUE. LGPARM = .NOT.LDCPLOT CALL PLTINI CALL PLOTG CALL PLOTC ENDIF XBOX(1) = XMARG XBOX(2) = XPAGE-XMARG YBOX(1) = YMARG YBOX(2) = YPAGE-YMARG XOFF1 = XOFF YOFF1 = (DYOFFP+YOFF)/YSFP XSF1 = XSF YSF1 = YSF*YSFP CALL MODIFY(NTX,1,NCAM,1,1, & XCAM,PCAM,PCAMP, LCSLOP, & K1,K2,ISMOD,IFMOD, & XBOX,YBOX, XBOX,YBOX,SIZE, & XOFF1,YOFF1,XSF1,YSF1, 'RED',' ', & NEWPLOTC) C C----- calculate camber line corresponding to specified loading CALL CPCAM(NCAM,XCAM,YCAM,YCAMP,PCAM,PCAMP) C C----- calculate added lift and moment from added loading CLX = 0.0 CMX = 0.0 DO K=1, NCAM-1 DX = XCAM(K+1) - XCAM(K) XA = 0.5*(XCAM(K+1) + XCAM(K)) PA = 0.5*(PCAM(K+1) + PCAM(K)) CLX = CLX + PA*DX CMX = CMX + PA*DX*(XCMREF-XA) END DO WRITE(*,1110) CLX, CMX C GO TO 200 C C-------------------------------------------------------- ELSEIF(COMAND.EQ.'CLR ') THEN C----- Clear the added camber LCLEAR = .TRUE. LDCPLOT = .FALSE. GO TO 100 C C-------------------------------------------------------- ELSEIF(COMAND.EQ.'SCAL') THEN C----- Scale camber IF(NINPUT.GE.1) THEN SCAL = RINPUT(1) ELSE SCAL = 1.0 CALL ASKR('Enter camber scaling factor^',SCAL) ENDIF C C--- Scale added camber user arrays DO I = 1, NCADD YCADD(I) = YCADD(I) *SCAL YPADD(I) = YPADD(I) *SCAL YCADDP(I) = YCADDP(I)*SCAL YPADDP(I) = YPADDP(I)*SCAL END DO C C--- Scale added camber arrays DO I = 1, NCAM YCAM(I) = YCAM(I) *SCAL YCAMP(I) = YCAMP(I)*SCAL PCAM(I) = PCAM(I) *SCAL PCAMP(I) = PCAMP(I)*SCAL END DO C C----- calculate added lift and moment from added loading CLX = 0.0 CMX = 0.0 DO K=1, NCAM-1 DX = XCAM(K+1) - XCAM(K) XA = 0.5*(XCAM(K+1) + XCAM(K)) PA = 0.5*(PCAM(K+1) + PCAM(K)) CLX = CLX + PA*DX CMX = CMX + PA*DX*(XCMREF-XA) END DO IF(CLX.NE.0.0 .AND. CMX.NE.0.0) WRITE(*,1110) CLX, CMX C COMOLD = COMAND ARGOLD = COMARG C C C----- go replot new shape and resume menu loop GO TO 200 C C-------------------------------------------------------- ELSEIF(COMAND.EQ.'ADD ') THEN C----- Add camber to camberline CALL SEGSPL(YCAM,YCAMP,XCAM,NCAM) C C----- go over each point, changing the camber line appropriately DO I=1, NB C------- coordinates of point on the opposite side with the same chord x value CALL SOPPS(SBOPP, SB(I),XB,XBP,YB,YBP,SB,NB,SBL) XBOPP = SEVAL(SBOPP,XB,XBP,SB,NB) YBOPP = SEVAL(SBOPP,YB,YBP,SB,NB) C C------- set present camber height OLDCAM = 0.5*(YB(I)+YBOPP)*XBCH/SBCH & - 0.5*(XB(I)+XBOPP)*YBCH/SBCH C C------- add on new camber CAM = OLDCAM & + SEVAL(XB(I),YCAM,YCAMP,XCAM,NCAM) C C------- set new y coordinate by changing camber & thickness appropriately W1(I) = CAM + 0.5*(YB(I)-YBOPP) END DO C DO I=1, NB YB(I) = W1(I) END DO 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 LDCPLOT = .FALSE. C---- reinitialize added camber to zero LCLEAR = .TRUE. GO TO 100 C C-------------------------------------------------------- ELSEIF(COMAND.EQ.'SLOP') THEN LCSLOP = .NOT.LCSLOP IF(LCSLOP) 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.'BLOW' .OR. & COMAND.EQ.'B ') THEN XWS = XWIND/SIZE YWS = YWIND/SIZE CALL OFFGET(XOFF,YOFF,XSF,YSF,XWS,YWS, .TRUE. , .TRUE. ) SF = MIN(XSF,YSF) XSF = SF YSF = SF GO TO 300 C C-------------------------------------------------------- ELSEIF(COMAND.EQ.'RESE' .OR. & COMAND.EQ.'R ') THEN CALL PLTINI CALL GOFINI CALL PLOTG cc CALL RESETSCL GO TO 300 C C-------------------------------------------------------- ELSEIF(COMAND.EQ.'DCPL') THEN C----- Toggle DCp plot flag LDCPLOT = .NOT.LDCPLOT GO TO 200 C C-------------------------------------------------------- ELSEIF(COMAND.EQ.'CPLI') THEN IF (NINPUT.GE.2) THEN YPMN = RINPUT(1) YPMX = RINPUT(2) ELSEIF(NINPUT.GE.1) THEN YPMIN = RINPUT(1) CALL ASKR('Enter max DCp^',YPMX) ELSE CALL ASKR('Enter min DCp^',YPMN) CALL ASKR('Enter max DCp^',YPMX) ENDIF IF(YPMX-YPMN.GT.0.0) THEN CALL AXISADJ(YPMN,YPMX,PSPAN,DXYP,NTICS) YPMIN = YPMN YPMAX = YPMX CALL GOFINI ENDIF GO TO 300 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 GOFINI GO TO 300 C C-------------------------------------------------------- ELSEIF(COMAND.EQ.'ANNO') THEN IF(LPLOT) THEN CALL ANNOT(CH) ELSE WRITE(*,*) 'No active plot to annotate' ENDIF GO TO 300 C C-------------------------------------------------------- ELSEIF(COMAND.EQ.'HARD') THEN IF(LPLOT) CALL PLEND LPLOT = .FALSE. CALL REPLOT(IDEVRP) GO TO 300 C C------------------------------------------------------- ELSE WRITE(*,8000) COMAND 8000 FORMAT(1X,A4,' command not recognized. Type a "?" for list') C ENDIF GO TO 500 C 1110 FORMAT(/' Delta Cp loading gives delta CL = ',F7.3, & /' delta CM = ',F7.3) C END ! CAMB SUBROUTINE NEWPLOTC CALL PLTINI CALL PLOTG CALL PLOTC RETURN END SUBROUTINE PLOTC C------------------------------------------------------ C Plots camber, thickness on its own axis above airfoil plot C Also plots deltaCP distribution above the camber,thickness C on its own axis above airfoil plot if LDCPLOT is set C------------------------------------------------------ INCLUDE 'XFOIL.INC' C DATA LMASK0, LMASK1, LMASK2, LMASK3 / -1, -32640, -30584, -21846 / C CALL GETCOLOR(ICOL0) CALL GETPAT(IPAT0) C CHS = 0.6*CHG C IF(LDCPLOT) THEN C----- current DCp loading is valid... plot it CALL NEWPEN(1) CHL = 1.5*CHG YOFFP = (DYOFFP+YOFF)/YSFP CALL GRDAIR(XGMIN,XGMAX,YPMIN,YPMAX,DXYG,DXYP,CHG,.FALSE.,.TRUE., & XOFF,XSF,YOFFP,YSFP*YSF, LMASK2) C CALL NEWCOLORNAME('RED') CALL NEWPEN(2) XLAB = (XPMIN -XOFF )*XSF - 4.0*CHL YLAB = (YPMAX-0.5*DXYP-YOFFP)*YSFP*YSF - 0.6*CHL CALL PLCHAR(XLAB,YLAB,CHL,' Cp',0.0,3) CALL PLMATH(XLAB,YLAB,CHL,'O ',0.0,3) C CALL XYLINE(NCAM,XCAM,PCAM,XOFF,XSF,YOFFP,YSFP*YSF,1) ccc CALL XYSYMB(NCAM,XCAM,PCAM,XOFF,XSF,YOFFP,YSFP*YSF,CHS,1) CALL NEWCOLOR(ICOL0) C cC----- plot derived camber line in dotted line c CALL NEWPAT(LMASK3) c ELSE cC----- plot specified camber line in solid line c CALL NEWPAT(LMASK0) ENDIF C CALL NEWPEN(1) CALL NEWCOLORNAME('RED') CALL XYLINE(NCAM,XCAM,YCAM,XOFF,XSF,YOFF+DYOFFC,YSF,1) ccc CALL XYSYMB(NCAM,XCAM,YCAM,XOFF,XSF,YOFF+DYOFFC,YSF,CHS,1) C CALL NEWCOLOR(ICOL0) CALL NEWPAT(IPAT0) CALL PLFLUSH C RETURN END ! PLOTC SUBROUTINE ZERCAM C----------------------------------------- C Zeros out camber of buffer airfoil C----------------------------------------- INCLUDE 'XFOIL.INC' C WRITE(*,*) 'Setting current camber to zero.' TFAC = 1.0 CFAC = 0.0 CALL THKCAM(TFAC,CFAC) C C---- make points exact mirror images CALL YSYM(XB,XBP,YB,YBP,SB,2*IQX,NB,1,W1,W2) 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 ! ZERCAM SUBROUTINE TCBUF C------------------------------------------------------ C Reports buffer airfoil thickness and camber C------------------------------------------------------ INCLUDE 'XFOIL.INC' C C--- find the current buffer airfoil camber and thickness CALL GETCAM(XCM,YCM,NCM,XTK,YTK,NTK, & XB,XBP,YB,YBP,SB,NB ) CALL GETMAX(XCM,YCM,YCMP,NCM,CXMAX,CYMAX) CALL GETMAX(XTK,YTK,YTKP,NTK,TXMAX,TYMAX) C WRITE(*,1000) 2.0*TYMAX,TXMAX, CYMAX,CXMAX 1000 FORMAT( ' Max thickness = ',F8.4,' at x = ',F7.3, & /' Max camber = ',F8.4,' at x = ',F7.3) C RETURN END ! TCBUF SUBROUTINE TCSCAL(RINPUT,NINPUT) DIMENSION RINPUT(*) C------------------------------------------------------ C Finds buffer airfoil thickness and/or camber, C plots thickness, camber and airfoil, C and scales t and/or c by user input factors C------------------------------------------------------ INCLUDE 'XFOIL.INC' C C--- find the current buffer airfoil camber and thickness CALL GETCAM(XCM,YCM,NCM,XTK,YTK,NTK, & XB,XBP,YB,YBP,SB,NB ) CALL GETMAX(XCM,YCM,YCMP,NCM,CXMAX,CYMAX) CALL GETMAX(XTK,YTK,YTKP,NTK,TXMAX,TYMAX) C WRITE(*,1000) 2.0*TYMAX,TXMAX, CYMAX,CXMAX C IF (NINPUT .GE. 2) THEN TFAC = RINPUT(1) CFAC = RINPUT(2) ELSEIF(NINPUT .GE. 1) THEN TFAC = RINPUT(1) IF(LGSYM) THEN WRITE(*,*) 'Symmetry enforced: Maintaining zero camber.' ELSE CFAC = 1.0 CALL ASKR('Enter new/old camber scale factor^',CFAC) ENDIF ELSE TFAC = 1.0 CALL ASKR( 'Enter new/old thickness scale factor^',TFAC) IF(LGSYM) THEN WRITE(*,*) 'Symmetry enforced: Maintaining zero camber.' ELSE CFAC = 1.0 CALL ASKR('Enter new/old camber scale factor^',CFAC) ENDIF ENDIF C ccc IF (TFAC.LT.0.0) TFAC = 0.0 CALL THKCAM(TFAC,CFAC) C CALL GETCAM(XCM,YCM,NCM,XTK,YTK,NTK, & XB,XBP,YB,YBP,SB,NB ) cc IPLT = 1 cc CALL PLOTC C 1000 FORMAT(/' Max thickness = ',F8.4,' at x = ',F7.3, & /' Max camber = ',F8.4,' at x = ',F7.3/) C RETURN END ! TCSCAL SUBROUTINE TCSET(RINPUT,NINPUT) DIMENSION RINPUT(*) C------------------------------------------------------ C Finds buffer airfoil thickness and/or camber, C plots thickness, camber and airfoil, C and scales t and/or c by user input factors C------------------------------------------------------ INCLUDE 'XFOIL.INC' C C--- find the current buffer airfoil camber and thickness CALL GETCAM(XCM,YCM,NCM,XTK,YTK,NTK, & XB,XBP,YB,YBP,SB,NB ) CALL GETMAX(XCM,YCM,YCMP,NCM,CXMAX,CYMAX) CALL GETMAX(XTK,YTK,YTKP,NTK,TXMAX,TYMAX) C WRITE(*,1000) 2.0*TYMAX,TXMAX, CYMAX,CXMAX 1000 FORMAT(/' Max thickness = ',F8.4,' at x = ',F7.3, & /' Max camber = ',F8.4,' at x = ',F7.3/) C cc IPLT = 0 cc CALL PLOTC C IF (NINPUT .GE. 2) THEN TNEW = RINPUT(1) CNEW = RINPUT(2) ELSEIF(NINPUT .GE. 1) THEN TNEW = RINPUT(1) IF(LGSYM) THEN WRITE(*,*) 'Symmetry enforced: Maintaining zero camber.' ELSE CNEW = 999 CALL ASKR('Enter new max camber to skip^',CNEW) ENDIF ELSE TNEW = 999 CALL ASKR('Enter new max thickness to skip^',TNEW) IF(LGSYM) THEN WRITE(*,*) 'Symmetry enforced: Maintaining zero camber.' ELSE CNEW = 999 CALL ASKR('Enter new max camber to skip^',CNEW) ENDIF ENDIF C CFAC = 1.0 TFAC = 1.0 IF(CYMAX.NE.0.0 .AND. CNEW.NE.999.0) CFAC = CNEW / ( CYMAX) IF(TYMAX.NE.0.0 .AND. TNEW.NE.999.0) TFAC = TNEW / (2.0*TYMAX) C C---- sanity checks on scaling factors IF(ABS(TFAC) .GT. 100.0 .OR. ABS(CFAC) .GT. 100.0) THEN WRITE(*,1100) TFAC, CFAC 1100 FORMAT(/' Questionable input...' & /' Implied scaling factors are:', F13.2,' x thickness' & /' ', F13.2,' x camber ') CALL ASKL('Apply scaling factors?^',OK) IF(.NOT.OK) THEN WRITE(*,*) 'No action taken' RETURN ENDIF ENDIF C ccc IF (TFAC.LT.0.0) TFAC = 0.0 CALL THKCAM(TFAC,CFAC) C CALL GETCAM(XCM,YCM,NCM,XTK,YTK,NTK, & XB,XBP,YB,YBP,SB,NB ) cc IPLT = 1 cc CALL PLOTC C RETURN END ! TCSET SUBROUTINE THKCAM(TFAC,CFAC) C--------------------------------------------------- C Changes buffer airfoil thickness and camber C--------------------------------------------------- INCLUDE 'XFOIL.INC' C CALL LEFIND(SBLE,XB,XBP,YB,YBP,SB,NB) C C---This fails miserably with sharp LE foils, tsk,tsk,tsk HHY 4/24/01 C---- set baseline vector normal to surface at LE point c DXC = -DEVAL(SBLE,YB,YBP,SB,NB) c DYC = DEVAL(SBLE,XB,XBP,SB,NB) c DSC = SQRT(DXC**2 + DYC**2) c DXC = DXC/DSC c DYC = DYC/DSC C C---Rational alternative 4/24/01 HHY 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)) CHORD = SQRT((XTE-XLE)**2 + (YTE-YLE)**2) C---- set unit chord-line vector DXC = (XTE-XLE) / CHORD DYC = (YTE-YLE) / CHORD C C---- go over each point, changing the y-thickness appropriately DO I=1, NB C------ coordinates of point on the opposite side with the same x value CALL SOPPS(SBOPP, SB(I),XB,XBP,YB,YBP,SB,NB,SBLE) XBOPP = SEVAL(SBOPP,XB,XBP,SB,NB) YBOPP = SEVAL(SBOPP,YB,YBP,SB,NB) C C------ set new y coordinate by changing camber & thickness appropriately XCAVG = ( 0.5*(XB(I)+XBOPP)*DXC + 0.5*(YB(I)+YBOPP)*DYC ) YCAVG = CFAC * ( 0.5*(YB(I)+YBOPP)*DXC - 0.5*(XB(I)+XBOPP)*DYC ) XCDEL = ( 0.5*(XB(I)-XBOPP)*DXC + 0.5*(YB(I)-YBOPP)*DYC ) YCDEL = TFAC * ( 0.5*(YB(I)-YBOPP)*DXC - 0.5*(XB(I)-XBOPP)*DYC ) C W1(I) = (XCAVG+XCDEL)*DXC - (YCAVG+YCDEL)*DYC W2(I) = (YCAVG+YCDEL)*DXC + (XCAVG+XCDEL)*DYC ENDDO C 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 RETURN END ! THKCAM SUBROUTINE HIPNT(RINPUT,NINPUT) DIMENSION RINPUT(*) C------------------------------------------------------ C Changes buffer airfoil C thickness and/or camber highpoint C------------------------------------------------------ INCLUDE 'XFOIL.INC' REAL XFN(5), YFN(5), YFNP(5), SFN(5) C C C--- Check chordline direction (should be unrotated for camber routines) C to function correctly 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)) AROT = ATAN2(YLE-YTE,XTE-XLE) / DTOR IF(ABS(AROT).GT.1.0) THEN WRITE(*,*) ' ' WRITE(*,*) 'Warning: HIGH does not work well on rotated foils' WRITE(*,*) 'Current chordline angle: ',AROT WRITE(*,*) 'Proceeding anyway...' ENDIF C C C---- find leftmost point location CALL XLFIND(SBL,XB,XBP,YB,YBP,SB,NB) XBL = SEVAL(SBL,XB,XBP,SB,NB) YBL = SEVAL(SBL,YB,YBP,SB,NB) C 10 CONTINUE C C---- find the current buffer airfoil camber and thickness CALL GETCAM(XCM,YCM,NCM,XTK,YTK,NTK, & XB,XBP,YB,YBP,SB,NB ) C C---- find the max thickness and camber CALL GETMAX(XCM,YCM,YCMP,NCM,CXMAX,CYMAX) CALL GETMAX(XTK,YTK,YTKP,NTK,TXMAX,TYMAX) C C C---- make a picture and get some input specs for mods cc IPLT = 0 cc CALL PLOTC WRITE(*,1010) 2.0*TYMAX,TXMAX, CYMAX,CXMAX 1010 FORMAT(/' Max thickness = ',F8.4,' at x = ',F7.3, & /' Max camber = ',F8.4,' at x = ',F7.3/) C IF (NINPUT .GE. 2) THEN THPNT = RINPUT(1) CHPNT = RINPUT(2) ELSEIF(NINPUT .GE. 1) THEN THPNT = RINPUT(1) IF(LGSYM) THEN WRITE(*,*) 'Symmetry enforced: Maintaining zero camber.' ELSE CHPNT = 0.0 CALL ASKR('Enter new camber highpoint x: ^',CHPNT) ENDIF ELSE THPNT = 0.0 CALL ASKR('Enter new thickness highpoint x: ^',THPNT) IF(LGSYM) THEN WRITE(*,*) 'Symmetry enforced: Maintaining zero camber.' ELSE CHPNT = 0.0 CALL ASKR('Enter new camber highpoint x: ^',CHPNT) ENDIF ENDIF C IF (THPNT.LE.0.0) THPNT = TXMAX IF (CHPNT.LE.0.0) CHPNT = CXMAX C C--- a simple cubic mapping function is used to map x/c to move highpoints C C the assumption is that a smooth function (cubic, given by the old and C new highpoint locations) maps the range 0-1 for x/c C into the range 0-1 for altered x/c distribution for the same y/c C thickness or camber (ie. slide the points smoothly along the x axis) C C--- shift thickness highpoint IF (THPNT .GT. 0.0) THEN XFN(1) = XTK(1) XFN(2) = TXMAX XFN(3) = XTK(NTK) YFN(1) = XTK(1) YFN(2) = THPNT YFN(3) = XTK(NTK) CALL SPLINA(YFN,YFNP,XFN,3) DO I = 1, NTK XTK(I) = SEVAL(XTK(I),YFN,YFNP,XFN,3) ENDDO ENDIF C C--- shift camber highpoint IF (CHPNT .GT. 0.0) THEN XFN(1) = XCM(1) XFN(2) = CXMAX XFN(3) = XCM(NCM) YFN(1) = XCM(1) YFN(2) = CHPNT YFN(3) = XCM(NCM) CALL SPLINA(YFN,YFNP,XFN,3) DO I = 1, NCM XCM(I) = SEVAL(XCM(I),YFN,YFNP,XFN,3) ENDDO ENDIF C cc IPLT = 1 cc CALL PLOTC C C CALL ASKL('Is this acceptable? ^',OK) C IF(.NOT.OK) GO TO 10 C C---- Make new airfoil from thickness and camber C new airfoil points are spaced to match the original C--- HHY 4/24/01 got rid of splining vs X,Y vs S (buggy), now spline Y(X) CALL SEGSPL(YTK,YTKP,XTK,NTK) CALL SEGSPL(YCM,YCMP,XCM,NCM) C C C---- for each orig. airfoil point setup new YB from camber and thickness DO 40 I=1, NB C C------ spline camber and thickness at original xb points YCC = SEVAL(XB(I),YCM,YCMP,XCM,NCM) YTT = SEVAL(XB(I),YTK,YTKP,XTK,NTK) C C------ set new y coordinate from new camber & thickness IF (SB(I) .LE. SBL) THEN YB(I) = YCC + YTT ELSE YB(I) = YCC - YTT ENDIF C---- Add Y-offset for original leftmost (LE) point to camber YB(I) = YB(I) + YBL 40 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 ! HIPNT SUBROUTINE GETCAM (XCM,YCM,NCM,XTK,YTK,NTK, & X,XP,Y,YP,S,N ) C------------------------------------------------------ C Finds camber and thickness C distribution for input airfoil C------------------------------------------------------ REAL XCM(*), YCM(*) REAL XTK(*), YTK(*) REAL X(*),XP(*),Y(*),YP(*),S(*) C CALL XLFIND(SL,X,XP,Y,YP,S,N) XL = SEVAL(SL,X,XP,S,N) YL = SEVAL(SL,Y,YP,S,N) C C---- go over each point, finding opposite points, getting camber and thickness DO 10 I=1, N C------ coordinates of point on the opposite side with the same x value CALL SOPPS(SOPP, S(I), X,XP,Y,YP,S,N,SL) XOPP = SEVAL(SOPP,X,XP,S,N) YOPP = SEVAL(SOPP,Y,YP,S,N) C C------ get camber and thickness XCM(I) = 0.5*(X(I)+XOPP) YCM(I) = 0.5*(Y(I)+YOPP) XTK(I) = 0.5*(X(I)+XOPP) YTK(I) = 0.5*(Y(I)-YOPP) YTK(I) = ABS(YTK(I)) c if (XOPP.gt.0.9) then c write(*,*) 'cm i,x,y ',i,xcm(i),ycm(i) c write(*,*) 'tk i,x,y ',i,xtk(i),ytk(i) c endif 10 CONTINUE C C---- Tolerance for nominally identical points TOL = 1.0E-5 * (S(N)-S(1)) ccc TOL = 1.0E-3 * (S(N)-S(1)) ! Bad bug -- was losing x=1.0 point C C---- Sort the camber points NCM = N+1 XCM(N+1) = XL YCM(N+1) = YL CALL SORTOL(TOL,NCM,XCM,YCM) C C--- Reorigin camber from LE so camberlines start at Y=0 4/24/01 HHY C policy now to generate camber independent of Y-offsets YOF = YCM(1) DO I = 1, NCM YCM(I) = YCM(I) - YOF END DO C C---- Sort the thickness points NTK = N+1 XTK(N+1) = XL YTK(N+1) = 0.0 CALL SORTOL(TOL,NTK,XTK,YTK) C RETURN END ! GETCAM SUBROUTINE GETMAX(X,Y,YP,N,XMAX,YMAX) REAL X(*), Y(*), YP(*) C------------------------------------------------ C Calculates camber or thickness highpoint C and x position C------------------------------------------------ C XLEN = X(N) - X(1) XTOL = XLEN * 1.0E-5 C CALL SEGSPL(Y,YP,X,N) C C---- get approx max point and rough interval size YMAX0 = Y(1) XMAX0 = X(1) DO 5 I = 2, N IF (ABS(Y(I)).GT.ABS(YMAX0)) THEN YMAX0 = Y(I) XMAX0 = 0.5*(X(I-1) + X(I)) DDX = 0.5*ABS(X(I+1) - X(I-1)) ENDIF 5 CONTINUE XMAX = XMAX0 C C---- do a Newton loop to refine estimate DO 10 ITER=1, 10 YMAX = SEVAL(XMAX,Y,YP,X,N) RES = DEVAL(XMAX,Y,YP,X,N) RESP = D2VAL(XMAX,Y,YP,X,N) IF (ABS(XLEN*RESP) .LT. 1.0E-6) GO TO 20 DX = -RES/RESP DX = SIGN( MIN(0.5*DDX,ABS(DX)) , DX) XMAX = XMAX + DX IF(ABS(DX) .LT. XTOL) GO TO 20 10 CONTINUE WRITE(*,*) & 'GETMAX: Newton iteration for max camber/thickness failed.' YMAX = YMAX0 XMAX = XMAX0 C 20 RETURN END ! GETMAX SUBROUTINE CPCAM(N,X,Y,DYDX,P,DPDX) REAL X(*), Y(*), DYDX(*), P(*), DPDX(*) C------------------------------------------------------------------ C Generates y(x) camberline from specified DCp(x) distribution. C C Input: N number of points C X(.) x array C P(.) DCp array C DPDX(.) dDCp/dx array C C Output: Y(.) y(x) array C DYDX(.) dy/dx array C------------------------------------------------------------------ C---- 1 / 4 pi DATA QOPI / 7.9577471545948E-02 / C C---- singular part of camber y(x) due to finite loadings P0,P1 at LE and TE C- dYSING/dX has logarithmic singularity at x=X0,X1 YSING(XT) = QOPI*P1*((XT-X1)*LOG(MAX((X1-XT)/(X1-X0),1.E-6)) - XT) & - QOPI*P0*((XT-X0)*LOG(MAX((XT-X0)/(X1-X0),1.E-6)) - XT) C P0 = P(1) P1 = P(N) C X0 = X(1) X1 = X(N) C C---- calculate Cauchy integral for y'(x) with removed singularity DO I=1, N write(*,'(1x,i4,3f10.4)') i, x(i), p(i), dpdx(i) !###@@@ DYDX(I) = 0.0 J = 1 IF(I.EQ.J) THEN YP1 = DPDX(J) ELSE YP1 = (P(J) - P(I)) / (X(J) - X(I)) ENDIF DO J=2, N IF(I.EQ.J) THEN YP2 = DPDX(J) ELSE YP2 = (P(J) - P(I)) / (X(J) - X(I)) ENDIF DYDX(I) = DYDX(I) + 0.5*(YP1+YP2)*(X(J)-X(J-1)) YP1 = YP2 END DO DYDX(I) = QOPI*DYDX(I) C C------ add on removed part of Cauchy integral, further leaving out the C- possible infinities at LE and TE so that y(x) can be safely splined. C- The infinities are analytically integrated, and added on to y(x) C- with the statement function YSING. IF(I.NE.1) THEN DYDX(I) = DYDX(I) & - QOPI*(P(I) - P0)*LOG(X(I) - X0) ENDIF IF(I.NE.N) THEN DYDX(I) = DYDX(I) & + QOPI*(P(I) - P1)*LOG(X1 - X(I)) ENDIF END DO C C---- integrate regular part of y'(x) from LE Y(1) = 0. DO I=2, N Y(I) = Y(I-1) & + 0.5*(DYDX(I) + DYDX(I-1))*(X(I) - X(I-1)) END DO C C---- add on singular part DO I=1, N Y(I) = Y(I) + YSING(X(I)) END DO C C---- add offset and angle of attack to get y(0) = y(1) = 0 Y0 = Y(1) Y1 = Y(N) DO I=1, N Y(I) = Y(I) & - Y0*(X1 -X(I))/(X1-X0) & - Y1*(X(I)-X0 )/(X1-X0) END DO C RETURN END ! CPCAM