aboutsummaryrefslogtreecommitdiff
path: root/orrs/src/mapgen.f
diff options
context:
space:
mode:
authorDimitri Sokolyuk <demon@dim13.org>2009-05-11 00:27:49 +0000
committerDimitri Sokolyuk <demon@dim13.org>2009-05-11 00:27:49 +0000
commit0d4f43d355de79178b1142e9735902cf641670b6 (patch)
tree2ced2323f6351db2a51090b3fd13eb11f69ff53f /orrs/src/mapgen.f
Xfoil 6.97
Diffstat (limited to 'orrs/src/mapgen.f')
-rwxr-xr-xorrs/src/mapgen.f225
1 files changed, 225 insertions, 0 deletions
diff --git a/orrs/src/mapgen.f b/orrs/src/mapgen.f
new file mode 100755
index 0000000..2af02e0
--- /dev/null
+++ b/orrs/src/mapgen.f
@@ -0,0 +1,225 @@
+ PROGRAM MAPGEN
+ PARAMETER (NMAX=257,NRX=101,NWX=101)
+ REAL ETA(NMAX), F(NMAX), U(NMAX), S(NMAX)
+ REAL UTR(NMAX), UTI(NMAX), VTR(NMAX), VTI(NMAX)
+ CHARACTER*48 FNAME
+ REAL AR(NRX,NWX), AI(NRX,NWX)
+ REAL RT(NRX),RTL(NRX), WS(NWX),WSL(NWX)
+ LOGICAL CONV(NRX,NWX)
+C
+ LST = 1
+ LRE = 1
+C
+ WRMAX = 0.25
+C
+ RESMAX = 0.1
+C
+C---- default profile parameters
+ N = 256
+ GEO = 1.02
+ ETAE = 14.0
+C
+ DO 5 IR=1, NRX
+ DO 4 IW=1, NWX
+ CONV(IR,IW) = .FALSE.
+ 4 CONTINUE
+ 5 CONTINUE
+C
+C---- generate or read in profile
+ CALL PFLGET(N,GEO,ETAE,ETA,F,U,S,H)
+C
+C
+ CALL ASKR('Enter lower log10(Rtheta)^',RT1L)
+ CALL ASKR('Enter upper log10(Rtheta)^',RT2L)
+ CALL ASKI('Enter number of log10(Rtheta )intervals^',NR)
+C
+ CALL ASKR('Enter lower log10(Wr*sqrt(Rtheta))^',WS1L)
+ CALL ASKR('Enter upper log10(Wr*sqrt(Rtheta))^',WS2L)
+ CALL ASKI('Enter number of log10(Wr) intervals^',NW)
+C
+ NRP = NR + 1
+ NWP = NW + 1
+C
+ IF(NRP.GT.NRX) STOP 'Array overflow'
+ IF(NWP.GT.NWX) STOP 'Array overflow'
+C
+ RT1 = 10.0 ** RT1L
+ RT2 = 10.0 ** RT2L
+ DO 10 IR=1, NRP
+ RTL(IR) = RT1L + (RT2L-RT1L)*FLOAT(IR-1)/FLOAT(NR)
+ RT(IR) = 10.0 ** RTL(IR)
+ 10 CONTINUE
+C
+ WS1 = 10.0 ** WS1L
+ WS2 = 10.0 ** WS2L
+ DO 15 IW=1, NWP
+ WSL(IW) = WS1L + (WS2L-WS1L)*FLOAT(IW-1)/FLOAT(NW)
+ WS(IW) = 10.0 ** WSL(IW)
+ 15 CONTINUE
+C
+C
+ CALL ASKR('Enter initial ar for lower Rtheta, upper Wr^',AR0)
+ CALL ASKR('Enter initial ai for lower Rtheta, upper Wr^',AI0)
+C
+C
+ CALL ASKS('Enter map output filename^',FNAME)
+ OPEN(19,FILE=FNAME,STATUS='NEW',FORM='UNFORMATTED')
+ WRITE(19) N, H
+ WRITE(19) (ETA(I),I=1, N)
+ WRITE(19) (U(I) ,I=1, N)
+ WRITE(19) (S(I) ,I=1, N)
+ WRITE(19) NRP, NWP
+ WRITE(19) (RTL(IR),IR=1,NRP)
+ WRITE(19) (WSL(IW),IW=1,NWP)
+C
+ IR1 = NRP
+ IR2 = 1
+ IRD = -1
+C
+ DO 100 IW=1, NWP
+ WRITE(6,2010)
+ 2010 FORMAT(/1X,'--------------------')
+ DO 90 IR=IR1, IR2, IRD
+C
+ WR = WS(IW)/SQRT(RT(IR))
+C
+ WRITE(6,2020) IR,IW, RT(IR), WR
+ 2020 FORMAT(/1X,2I4,' Rth =', E12.4, ' Wr =', E12.4)
+C
+ WR0 = WR
+ WI0 = 0.0
+C
+C-------- set initial wavenumber guess
+ IRM1 = IR - IRINCR
+ IRM2 = IR - 2*IRINCR
+ IRM3 = IR - 3*IRINCR
+C
+ IWM1 = IW - IWINCR
+ IWM2 = IW - 2*IWINCR
+ IWM3 = IW - 3*IWINCR
+C
+ IF(IRM2.GE.1 .AND. IRM2.LE.NRP .AND.
+ & IWM1.GE.1 .AND. IWM1.LE.NWP ) THEN
+ AR0 = 2.0*AR(IRM1,IW ) - AR(IRM2,IW )
+ & + AR(IR ,IWM1) - 2.0*AR(IRM1,IWM1) + AR(IRM2,IWM1)
+ AI0 = 2.0*AI(IRM1,IW ) - AI(IRM2,IW )
+ & + AI(IR ,IWM1) - 2.0*AI(IRM1,IWM1) + AI(IRM2,IWM1)
+ ELSE IF(IRM1.GE.1 .AND. IRM1.LE.NRP .AND.
+ & IWM2.GE.1 .AND. IWM2.LE.NWP ) THEN
+ AR0 = AR(IRM1,IW )
+ & + 2.0*AR(IR ,IWM1) - 2.0*AR(IRM1,IWM1)
+ & - AR(IR ,IWM2) + AR(IRM1,IWM2)
+ AI0 = AI(IRM1,IW )
+ & + 2.0*AI(IR ,IWM1) - 2.0*AI(IRM1,IWM1)
+ & - AI(IR ,IWM2) + AI(IRM1,IWM2)
+ ELSE IF(IRM1.GE.1 .AND. IRM1.LE.NRP .AND.
+ & IWM1.GE.1 .AND. IWM1.LE.NWP ) THEN
+ AR0 = AR(IRM1,IW )
+ & + AR(IR ,IWM1) - AR(IRM1,IWM1)
+ AI0 = AI(IRM1,IW )
+ & + AI(IR ,IWM1) - AI(IRM1,IWM1)
+ ELSE IF(IRM2.GE.1 .AND. IRM2.LE.NRP) THEN
+ AR0 = 2.0*AR(IRM1,IW) - AR(IRM2,IW)
+ AI0 = 2.0*AI(IRM1,IW) - AI(IRM2,IW)
+ ELSE IF(IWM2.GE.1 .AND. IWM2.LE.NWP) THEN
+ AR0 = 2.0*AR(IR,IWM1) - AR(IR,IWM2)
+ AI0 = 2.0*AI(IR,IWM1) - AI(IR,IWM2)
+ ELSE IF(IRM1.GE.1 .AND. IRM1.LE.NRP) THEN
+ AR0 = AR(IRM1,IW)
+ AI0 = AI(IRM1,IW)
+ ELSE IF(IWM1.GE.1 .AND. IWM1.LE.NWP) THEN
+ AR0 = AR(IR,IWM1)
+ AI0 = AI(IR,IWM1)
+CCC ELSE
+CCC STOP 'Cannot start in corner and go in'
+ ENDIF
+c
+ AR(IR,IW) = AR0
+ AI(IR,IW) = AI0
+C
+C-------- don't bother with absurdly high frequency
+ IF(WR .GE. WRMAX) THEN
+ DELMAX = 0.0
+ GO TO 89
+ ENDIF
+C
+ ITMAX = 10
+ CALL ORRS(LST,LRE,N,ETA,U,S, RT(IR), ITMAX,
+ & AR0,AI0, WR0,WI0, UTR,UTI,VTR,VTI,DELMAX)
+C
+ 89 IF(DELMAX.LT.RESMAX) CONV(IR,IW) = .TRUE.
+C
+ AR(IR,IW) = AR0
+ AI(IR,IW) = AI0
+C
+ 90 CONTINUE
+C
+ WRITE(19) (AR(IR,IW),IR=1,NRP)
+ WRITE(19) (AI(IR,IW),IR=1,NRP)
+C
+ 100 CONTINUE
+C
+ CLOSE(19)
+C
+ STOP
+ END
+
+
+ SUBROUTINE PFLGET(N,GEO,ETAE,ETA,F,U,S,H)
+ DIMENSION ETA(N),F(N),U(N),S(N)
+ CHARACTER*48 FNAME
+C
+C---- eta coordinate normalized with momentum thickness
+ INORM = 3
+C
+ WRITE(6,*) ' '
+ WRITE(6,*) ' 1 Falkner-Skan parameter m = x/U dU/dx'
+ WRITE(6,*) ' 2 Falkner-Skan parameter beta = 2m/(m+1)'
+ WRITE(6,*) ' 3 Falkner-Skan shape parameter H'
+ WRITE(6,*) ' 4 General profile input file'
+ WRITE(6,*) ' '
+ CALL ASKI('Select profile option^',IOPT)
+C
+ IF(IOPT.NE.4) THEN
+ CALL ASKI('Enter number of BL points^',N)
+ CALL ASKR('Enter geometric stretching factor^',GEO)
+ CALL ASKR('Enter edge eta value^',ETAE)
+ ENDIF
+C
+C
+ IF(IOPT.EQ.1) THEN
+C
+ CALL ASKR('Enter m^',BU)
+ CALL FS(INORM,1,BU,H,N,ETAE,GEO,ETA,F,U,S)
+C
+ ELSE IF(IOPT.EQ.2) THEN
+C
+ CALL ASKR('Enter beta^',BETA)
+ BU = BETA/(2.0-BETA)
+ CALL FS(INORM,1,BU,H,N,ETAE,GEO,ETA,F,U,S)
+C
+ ELSE IF(IOPT.EQ.3) THEN
+C
+ CALL ASKR('Enter H^',H)
+ CALL FS(INORM,2,BU,H,N,ETAE,GEO,ETA,F,U,S)
+C
+ ELSE
+C
+ CALL ASKS('Enter profile filename^',FNAME)
+ OPEN(1,FILE=FNAME,STATUS='OLD')
+ READ(1,*) N, H
+ DO 5 I=1, N
+ READ(1,*) ETA(I), U(I), S(I)
+ 5 CONTINUE
+ CLOSE(1)
+C
+ GEO = (ETA(3)-ETA(2)) / (ETA(2)-ETA(1))
+ ENDIF
+C
+ WRITE(6,1050) N, H, ETA(N), GEO
+ 1050 FORMAT(/' n =', I4,' H =', F7.3,
+ & ' Ye =', F7.3,
+ & ' dYi+1/dYi =',F6.3 /)
+C
+ RETURN
+ END