aboutsummaryrefslogtreecommitdiff
path: root/src/dplot1.f
diff options
context:
space:
mode:
Diffstat (limited to 'src/dplot1.f')
-rw-r--r--src/dplot1.f288
1 files changed, 288 insertions, 0 deletions
diff --git a/src/dplot1.f b/src/dplot1.f
new file mode 100644
index 0000000..1a15038
--- /dev/null
+++ b/src/dplot1.f
@@ -0,0 +1,288 @@
+C***********************************************************************
+C Module: dplot.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***********************************************************************
+
+ SUBROUTINE DPLOT(NPR1,XPR,YPR)
+ INCLUDE 'XFOIL.INC'
+C-----------------------------------------------------------
+C Plots analytical profiles at specified points.
+C If NPR=0, then cursor-selected points are requested.
+C-----------------------------------------------------------
+ DIMENSION XPR(*), YPR(*)
+C
+ CHARACTER*1 KCHAR
+ LOGICAL LCRS, TURB
+ LOGICAL LGUI
+C
+ CALL GETCOLOR(ICOL0)
+C
+ LCRS = NPR1 .LE. 0
+C
+ IF(LCRS) THEN
+ KDONE = 1
+ XDWIN = XPAGE - 2.0*XMARG
+ YDWIN = YPAGE - 2.0*YMARG
+ X1 = XMARG + 0.91*XDWIN
+ X2 = XMARG + 0.99*XDWIN
+ Y1 = YMARG + 0.01*YDWIN
+ Y2 = YMARG + 0.05*YDWIN
+ CALL NEWPEN(5)
+ CALL GUIBOX(KDONE, X1,X2,Y1,Y2, 'GREEN' , ' Done ')
+C
+ WRITE(*,*) ' '
+ WRITE(*,*) 'Locate profiles with cursor, type "D" when done...'
+ NPR = 12345
+C
+ ELSE
+ NPR = NPR1
+C
+ ENDIF
+C
+C---- go over profiles ...
+ DO 50 IPR=1, NPR
+C
+ IF(LCRS) THEN
+C------- get cursor plot coordinates
+ CALL GETCURSORXY(XC,YC,KCHAR)
+ IF(INDEX('Dd',KCHAR).NE.0 .OR. LGUI(KDONE,XC,YC)) THEN
+ RETURN
+ ENDIF
+C
+C------- transform to airfoil coordinates
+ XC = XC/FACA - XOFA
+ YC = YC/FACA - YOFA
+C
+ ELSE
+ XC = XPR(IPR)
+ YC = YPR(IPR)
+C
+ ENDIF
+C
+C------ find nearest airfoil surface point
+ RSQMIN = 1.0E23
+ ISMIN = 0
+ IBLMIN = 0
+ DOFF = 0.00001*(S(N)-S(1))
+ DO IS = 1, 2
+ DO IBL = 2, IBLTE(IS)
+ I = IPAN(IBL,IS)
+ XSURF = X(I) + DOFF*YP(I)
+ YSURF = Y(I) - DOFF*XP(I)
+ RSQ = (XC-XSURF)**2 + (YC-YSURF)**2
+ IF(RSQ .LE. RSQMIN) THEN
+ RSQMIN = RSQ
+ ISMIN = IS
+ IBLMIN = IBL
+ ENDIF
+ ENDDO
+ ENDDO
+C
+ IS = ISMIN
+ IBL = IBLMIN
+C
+ I = IPAN(IBL,IS)
+ CRSP = (XC-X(I))*NY(I) - (YC-Y(I))*NX(I)
+ IF(IS.EQ.2) CRSP = -CRSP
+C
+ IF(CRSP.GT.0.0) THEN
+ IBLP = IBL+1
+ IBLO = IBL
+ ELSE
+ IBLP = IBL
+ IBLO = IBL-1
+ ENDIF
+ ISP = IS
+ ISO = IS
+C
+ IF(IBLP.GT.IBLTE(IS)) THEN
+ IBLP = IBLTE(IS)
+ IBLO = IBLP-1
+ IBL = IBLTE(IS)
+ ELSEIF(IBLO.LT.2) THEN
+ IBLO = 2
+ IF(ISO.EQ.1) THEN
+ ISO = 2
+ ELSE
+ ISO = 1
+ ENDIF
+ ENDIF
+C
+ IP = IPAN(IBLP,ISP)
+ IO = IPAN(IBLO,ISO)
+C
+C------ set interpolation fraction at profile location
+ DX = X(IP) - X(IO)
+ DY = Y(IP) - Y(IO)
+ VX = XC - X(IO)
+ VY = YC - Y(IO)
+ FRAC = (DX*VX + DY*VY)/(DX*DX+DY*DY)
+ FRAC = MIN( MAX( FRAC , 0.0 ) , 1.0 )
+C
+C------ set averaged displacement vector at profile location
+ CA = FRAC*NY(IP) + (1.0-FRAC)*NY(IO)
+ SA = FRAC*NX(IP) + (1.0-FRAC)*NX(IO)
+ CSMOD = SQRT(CA**2 + SA**2)
+ CA = CA/CSMOD
+ SA = SA/CSMOD
+C
+ X0 = FRAC*X(IP) + (1.0-FRAC)*X(IO)
+ Y0 = FRAC*Y(IP) + (1.0-FRAC)*Y(IO)
+C
+ DS = FRAC*DSTR(IBLP,ISP) + (1.0-FRAC)*DSTR(IBLO,ISO)
+ TH = FRAC*THET(IBLP,ISP) + (1.0-FRAC)*THET(IBLO,ISO)
+ UE = FRAC*UEDG(IBLP,ISP) + (1.0-FRAC)*UEDG(IBLO,ISO)
+C
+ XI = FRAC*XSSI(IBLP,ISP) + (1.0-FRAC)*XSSI(IBLO,ISO)
+ TURB = XI .GT. XSSITR(IS)
+C
+C------ 1 / (total enthalpy)
+ HSTINV = GAMM1*(MINF/QINF)**2 / (1.0 + 0.5*GAMM1*MINF**2)
+C
+C------ Sutherland's const./To (assumes stagnation conditions are at STP)
+ HVRAT = 0.35
+C
+C------ fill Rtheta arrays
+ UEC = UE * (1.0-TKLAM) / (1.0 - TKLAM*(UE/QINF)**2)
+ HERAT = (1.0 - 0.5*HSTINV*UEC **2)
+ & / (1.0 - 0.5*HSTINV*QINF**2)
+ RHOE = HERAT ** (1.0/GAMM1)
+ AMUE = SQRT(HERAT**3) * (1.0+HVRAT)/(HERAT+HVRAT)
+ RTHETA = REINF * RHOE*UE*TH/AMUE
+C
+ AMSQ = UEC*UEC*HSTINV / (GAMM1*(1.0 - 0.5*UEC*UEC*HSTINV))
+ CALL HKIN( DS/TH, AMSQ, HK, DUMMY, DUMMY)
+C
+ WRITE(*,9100) X0,Y0, DS, RTHETA, HK
+ 9100 FORMAT(1X,'x y =', 2F8.4,' Delta* =', G12.4,
+ & ' Rtheta =', F10.2,' Hk =', F9.4)
+C
+ IF(IS.EQ.1) THEN
+ UDIR = 1.0
+ ELSE
+ UDIR = -1.0
+ ENDIF
+C
+ UEI = UE/QINF
+ UN = 0.0
+ CALL NEWCOLORNAME('green')
+ UPRWTS = UPRWT*0.5*(S(N)-S(1))
+ CALL PRPLOT(X0,Y0,TH,UEI,UN,HK,RTHETA,AMSQ,TURB,
+ & -XOFA,-YOFA,FACA,UPRWTS,SA,CA,UDIR)
+ 50 CONTINUE
+C
+ CALL NEWCOLOR(ICOL0)
+ CALL PLFLUSH
+C
+ RETURN
+ END ! DPLOT
+
+
+
+ SUBROUTINE PRPLOT(X0,Y0,TH,UE,UN,HK,RET,MSQ,TURB,
+ & XOFA,YOFA,FACA,UWT,SINA,COSA,UDIR)
+C-----------------------------------------------------------------
+C Plots velocity profile taken from flow solution.
+C
+C X0,Y0 coordinates of point through which profile axis passes
+C SA,CA sin,cos of profile axis angle (cw from vertical)
+C-----------------------------------------------------------------
+ REAL MSQ
+ LOGICAL TURB
+C
+ PARAMETER (KPRX=129)
+ DIMENSION XX(KPRX), YY(KPRX), FFS(KPRX), SFS(KPRX)
+c
+ XMOD(XTMP) = FACA * (XTMP - XOFA)
+ YMOD(YTMP) = FACA * (YTMP - YOFA)
+C
+ NN = KPRX
+ UO = 1.0
+ DK = HK*TH
+ CT = 0.
+C
+ IF(TURB) THEN
+C------ set Spalding + power-law turbulent profile
+ CALL PRWALL(DK,TH,UO,RET,MSQ,CT, BB,
+ & DE, DE_DS, DE_TH, DE_UO, DE_RT, DE_MS,
+ & US, US_DS, US_TH, US_UO, US_RT, US_MS,
+ & HS, HS_DS, HS_TH, HS_UO, HS_RT, HS_MS,
+ & CF, CF_DS, CF_TH, CF_UO, CF_RT, CF_MS,
+ & CD, CD_DS, CD_TH, CD_UO, CD_RT, CD_MS,
+ & CD_CT )
+c
+ CALL UWALL(TH,UO,DE,US,RET,CF,BB, YY,XX,NN)
+C
+C------ limit profile height
+ DECORR = 1.5 * (3.15 + 1.72/(HK-1.0) + HK) * TH
+ DO 422 K=NN, 1, -1
+ IF(YY(K) .LE. DECORR) GO TO 423
+ 422 CONTINUE
+ 423 NN = K
+ DE = YY(K)
+C
+ ELSE
+C------ set Falkner-Skan profile
+ INORM = 3
+ ISPEC = 2
+ HSPEC = HK
+ ETAE = 1.5*(3.15 + 1.72/(HK-1.0) + HK)
+ GEO = 1.0
+ CALL FS(INORM,ISPEC,BU,HSPEC,NN,ETAE,GEO,YY,FFS,XX,SFS,DEFS)
+ DE = ETAE*TH
+C
+ DO 425 K=1, NN
+ YY(K) = YY(K)*TH
+ 425 CONTINUE
+C
+ ENDIF
+C
+ YAX = 1.1*DE
+C
+ X1 = X0
+ Y1 = Y0
+ X2 = X0 + YAX*SINA
+ Y2 = Y0 + YAX*COSA
+C
+C---- plot axis
+ CALL NEWPEN(1)
+ CALL PLOT(XMOD(X1),YMOD(Y1),3)
+ CALL PLOT(XMOD(X2),YMOD(Y2),2)
+C
+ DO K=1, NN
+ ULOC = UE + UN*(YY(K)-DK)
+ XX(K) = XX(K)*UE * UWT * UDIR
+CCC YY(K) = YY(K)
+ ENDDO
+C
+C---- rotate and position profile
+ DO K=1, NN
+ XBAR = XX(K)
+ YBAR = YY(K)
+ XROT = XBAR*COSA + YBAR*SINA + X0
+ YROT = YBAR*COSA - XBAR*SINA + Y0
+ XX(K) = XMOD(XROT)
+ YY(K) = YMOD(YROT)
+ ENDDO
+C
+ CALL NEWPEN(2)
+ CALL XYLINE(NN,XX,YY,0.0,1.0,0.0,1.0,1)
+C
+ RETURN
+ END ! PRPLOT
+