From 0d4f43d355de79178b1142e9735902cf641670b6 Mon Sep 17 00:00:00 2001 From: Dimitri Sokolyuk Date: Mon, 11 May 2009 00:27:49 +0000 Subject: Xfoil 6.97 --- osrc/getosfile.c | 42 +++++ osrc/osmap.f | 492 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 534 insertions(+) create mode 100644 osrc/getosfile.c create mode 100755 osrc/osmap.f (limited to 'osrc') diff --git a/osrc/getosfile.c b/osrc/getosfile.c new file mode 100644 index 0000000..a1d0358 --- /dev/null +++ b/osrc/getosfile.c @@ -0,0 +1,42 @@ +#include +#include +#include + +/* Handle various system requirements for trailing underscores, or other + fortran-to-C interface shenanigans thru defines for routine names + The provided set gives the option of setting a compile flag -DUNDERSCORE + to include underscores on C routine name symbols */ + +#ifdef UNDERSCORE + +#define GETOSFILE getosfile_ + +#endif + +void +GETOSFILE(osfile,len) + char *osfile; + int *len; +{ + char *bufp; + int l; + +/* get environment variable OSMAP for location of OS map data file */ + bufp = getenv("OSMAP"); + + /* printf("bufp: %s\n",bufp); + printf("osfile: %s\n",osfile); */ + + if(bufp){ + strcpy(osfile,bufp); + l = strlen(bufp); + } + else { + l = 0; + } + + *len = l; + + /* printf("len %d\n",*len); */ +} + diff --git a/osrc/osmap.f b/osrc/osmap.f new file mode 100755 index 0000000..7b9c471 --- /dev/null +++ b/osrc/osmap.f @@ -0,0 +1,492 @@ + SUBROUTINE OSMAP(RSP,WSP,HSP, + & ALFR, + & ALFR_R, ALFR_W, ALFR_H, + & ALFRW_R,ALFRW_W,ALFRW_H , + & ALFI, + & ALFI_R, ALFI_W, ALFI_H, + & ALFIW_R,ALFIW_W,ALFIW_H , OK) +C--------------------------------------------------------------------- +C +C Returns real and imaginary parts of complex wavenumber (Alpha) +C eigenvalue from Orr-Sommerfeld spatial-stability solution +C with mean profiles characterized by shape parameter H. +C Also returns the sensitivities of Alpha with respect to the +C input parameters. +C +C The eigenvalue Alpha(Rtheta,W,H) is stored as a 3-D array at +C discrete points, which is then interpolated to any (Rtheta,W,H) +C via a tricubic spline. The spline coordinates actually used are: +C +C RL = log10(Rtheta) +C WL = log10(W) + 0.5 log10(Rtheta) +C HL = H +C +C +C Input: +C ------ +C RSP momentum thickness Reynolds number Rtheta = Theta Ue / v +C WSP normalized disturbance frequency W = w Theta/Ue +C HSP shape parameter of mean profile H = Dstar/Theta +C +C Output: +C ------- +C ALFR real part of complex wavenumber * Theta +C ALFR_R d(ALFR)/dRtheta +C ALFR_W d(ALFR)/dW +C ALFR_H d(ALFR)/dH +C ALFRW_R d(dALFR/dW)/dRtheta +C ALFRW_W d(dALFR/dW)/dW +C ALFRW_H d(dALFR/dW)/dH +C +C ALFI imag part of complex wavenumber * Theta +C ALFI_R d(ALFI)/dRtheta +C ALFI_W d(ALFI)/dW +C ALFI_H d(ALFI)/dH +C ALFIW_R d(dALFI/dW)/dRtheta +C ALFIW_W d(dALFI/dW)/dW +C ALFIW_H d(dALFI/dW)/dH +C +C OK T if look up was successful; all values returned are valid +C F if point fell outside (RL,WL) spline domain limits; +C all values (ALFR, ALFR_R, etc.) are returned as zero. +C Exception: If points only falls outside HL spline limits, +C then the HL limit is used and an ALFR value is calculated, +C but OK is still returned as F. +C +C--------------------------------------------------------------------- + LOGICAL OK +C +C + REAL B(2,2), BR(2,2), BW(2,2), BH(2,2), + & BRW(2,2),BRH(2,2),BWH(2,2),BRWH(2,2) + REAL C(2) , CR(2) , CW(2) , CH(2) , + & CRW(2) ,CRH(2) ,CWH(2) ,CRWH(2) +C + REAL AINT(2), + & AINT_R(2), AINT_W(2), AINT_H(2), + & AINTW_R(2),AINTW_W(2),AINTW_H(2) +C + PARAMETER (NRX=31, NWX=41, NHX=21) + COMMON /AICOM_I/ NR, NW, NH, + & IC1, IC2, + & IW1(NHX), IW2(NHX), IR1(NHX),IR2(NHX) +C +C--------------------------------------------------------------- +C---- single-precision OS data file + REAL*4 RLSP, WLSP, HLSP, + & RINCR, WINCR, RL, WL, HL, + & A, AR, AW, AH, ARW, ARH, AWH, ARWH +C +C---- native-precision OS data file +c REAL RLSP, WLSP, HLSP, +c & RINCR, WINCR, RL, WL, HL, +c & A, AR, AW, AH, ARW, ARH, AWH, ARWH +C--------------------------------------------------------------- +C + COMMON /AICOM_R/ RINCR, WINCR, RL(NRX), WL(NWX), HL(NHX), + & A(NRX,NWX,NHX,2), + & AR(NRX,NWX,NHX,2), + & AW(NRX,NWX,NHX,2), + & AH(NRX,NWX,NHX,2), + & ARW(NRX,NWX,NHX,2), + & ARH(NRX,NWX,NHX,2), + & AWH(NRX,NWX,NHX,2), + & ARWH(NRX,NWX,NHX,2) +C + LOGICAL LOADED, NOFILE + SAVE LOADED, NOFILE +C +C---- set OSFILE to match the absolute OS database filename + CHARACTER*128 OSFILE + INTEGER LOSF + + DATA OSFILE / '/var/local/codes/orrs/osmap.dat' / +c DATA OSFILE +c &/'/afs/athena.mit.edu/course/16/16_d0006/Codes/orrs/osmap_lx.dat'/ +C + DATA LOADED, NOFILE / .FALSE. , .FALSE. / +C +C---- set ln(10) for derivatives of log10 function + DATA AL10 /2.302585093/ +C +C +C---- set default returned variables in case of error, or OS map not available + ALFR = 0.0 + ALFR_R = 0.0 + ALFR_W = 0.0 + ALFR_H = 0.0 + ALFRW_R = 0.0 + ALFRW_W = 0.0 + ALFRW_H = 0.0 +C + ALFI = 0.0 + ALFI_R = 0.0 + ALFI_W = 0.0 + ALFI_H = 0.0 + ALFIW_R = 0.0 + ALFIW_W = 0.0 + ALFIW_H = 0.0 +C + OK = .FALSE. +C + IF(NOFILE) RETURN +C + IF(LOADED) GO TO 9 +C-------------------------------------------------------------------- +C---- first time OSMAP is called ... load in 3-D spline data +C + CALL GETOSFILE(OSFILE,LOSF) + IF(LOSF.EQ.0) GO TO 800 +C + NR = 0 + NW = 0 + NH = 0 +C + LU = 31 + OPEN(UNIT=LU,FILE=OSFILE(1:LOSF),STATUS='OLD', + & FORM='UNFORMATTED',ERR=900) +C + READ(LU) NR, NW, NH +C + IF(NR.GT.NRX .OR. + & NW.GT.NWX .OR. + & NH.GT.NHX ) THEN + WRITE(*,*) 'OSMAP: Array limit exceeded.' + IF(NR.GT.NRX) WRITE(*,*) ' Increase NRX to', NR + IF(NW.GT.NWX) WRITE(*,*) ' Increase NWX to', NW + IF(NH.GT.NHX) WRITE(*,*) ' Increase NHX to', NH + STOP + ENDIF +C + READ(LU) (RL(IR), IR=1,NR) + READ(LU) (WL(IW), IW=1,NW) + READ(LU) (HL(IH), IH=1,NH) + READ(LU) (IR1(IH),IR2(IH),IW1(IH),IW2(IH), IH=1,NH) +C + DO IC = 2, 1, -1 + DO IH=1, NH + DO IW=IW1(IH), IW2(IH) + READ(LU,END=5) + & ( A(IR,IW,IH,IC), IR=IR1(IH),IR2(IH)) + READ(LU) ( AR(IR,IW,IH,IC), IR=IR1(IH),IR2(IH)) + READ(LU) ( AW(IR,IW,IH,IC), IR=IR1(IH),IR2(IH)) + READ(LU) ( AH(IR,IW,IH,IC), IR=IR1(IH),IR2(IH)) + READ(LU) ( ARW(IR,IW,IH,IC), IR=IR1(IH),IR2(IH)) + READ(LU) ( ARH(IR,IW,IH,IC), IR=IR1(IH),IR2(IH)) + READ(LU) ( AWH(IR,IW,IH,IC), IR=IR1(IH),IR2(IH)) + READ(LU) (ARWH(IR,IW,IH,IC), IR=IR1(IH),IR2(IH)) + ENDDO + ENDDO + ENDDO +C + 5 CONTINUE + IF(IH.LT.NH) THEN +C----- only imaginary part is available + IC1 = 2 + IC2 = 2 + ELSE +C----- both real and imaginary parts available + IC1 = 1 + IC2 = 2 + ENDIF + CLOSE(LU) +C +C + RINCR = (RL(NR) - RL(1))/FLOAT(NR-1) + WINCR = (WL(NW) - WL(1))/FLOAT(NW-1) + LOADED = .TRUE. +C +C-------------------------------------------------------------------- + 9 CONTINUE +C +C + IF(NR.EQ.0 .OR. NW.EQ.0 .OR. NH.EQ.0) THEN +C----- map not available for some reason (OPEN or READ error on osmap.dat?) + OK = .FALSE. + RETURN + ENDIF +C +C---- define specified spline coordinates + RLSP = ALOG10(RSP) + WLSP = ALOG10(WSP) + 0.5*RLSP + HLSP = HSP +C +C---- assume map limits will not be exceeded + OK = .TRUE. +C +C---- find H interval + DO IH = 2, NH + IF(HL(IH) .GE. HLSP) GO TO 11 + ENDDO + IH = NH + 11 CONTINUE +C + IF(HLSP.LT.HL(1) .OR. HLSP.GT.HL(NH)) THEN +CCC OK = .FALSE. +CCC WRITE(*,*) 'Over H limits. R w H:', RSP,WSP,HSP +CCC RETURN + HLSP = MAX( HL(1) , MIN( HL(NH) , HLSP ) ) + ENDIF +C +C---- find R interval + IR = INT((RLSP-RL(1))/RINCR + 2.001) + IR1X = MAX( IR1(IH) , IR1(IH-1) ) + IR2X = MIN( IR2(IH) , IR2(IH-1) ) + IF(IR-1.LT.IR1X .OR. IR.GT.IR2X) THEN + OK = .FALSE. +CCC WRITE(*,*) 'Over R limits. R w H:', RSP,WSP,HSP +CCC RETURN + IR = MAX( IR1X+1 , MIN( IR2X , IR ) ) + RLSP = MAX( RL(1) , MIN( RL(NR) , RLSP ) ) + ENDIF +C +C---- find W interval + IW = INT((WLSP-WL(1))/WINCR + 2.001) + IW1X = MAX( IW1(IH) , IW1(IH-1) ) + IW2X = MIN( IW2(IH) , IW2(IH-1) ) + IF(IW-1.LT.IW1X .OR. IW.GT.IW2X) THEN + OK = .FALSE. +CCC WRITE(*,*) 'Over w limits. R w H:', RSP,WSP,HSP +CCC RETURN + IW = MAX( IW1X+1 , MIN( IW2X , IW ) ) + WLSP = MAX( WL(1) , MIN( WL(NW) , WLSP ) ) + ENDIF +C + DRL = RL(IR) - RL(IR-1) + DWL = WL(IW) - WL(IW-1) + DHL = HL(IH) - HL(IH-1) + TR = (RLSP - RL(IR-1)) / DRL + TW = (WLSP - WL(IW-1)) / DWL + TH = (HLSP - HL(IH-1)) / DHL +C + TR = MAX( 0.0 , MIN( 1.0 , TR ) ) + TW = MAX( 0.0 , MIN( 1.0 , TW ) ) + TH = MAX( 0.0 , MIN( 1.0 , TH ) ) +C +C---- compute real and imaginary parts + DO 1000 IC = IC1, IC2 +C +C---- evaluate spline in Rtheta at the corners of HL,WL cell + DO 20 KH=1, 2 + JH = IH + KH-2 + DO 205 KW=1, 2 + JW = IW + KW-2 + A1 = A (IR-1,JW,JH,IC) + AR1 = AR (IR-1,JW,JH,IC) + AW1 = AW (IR-1,JW,JH,IC) + AH1 = AH (IR-1,JW,JH,IC) + ARW1 = ARW (IR-1,JW,JH,IC) + ARH1 = ARH (IR-1,JW,JH,IC) + AWH1 = AWH (IR-1,JW,JH,IC) + ARWH1 = ARWH(IR-1,JW,JH,IC) +C + A2 = A (IR ,JW,JH,IC) + AR2 = AR (IR ,JW,JH,IC) + AW2 = AW (IR ,JW,JH,IC) + AH2 = AH (IR ,JW,JH,IC) + ARW2 = ARW (IR ,JW,JH,IC) + ARH2 = ARH (IR ,JW,JH,IC) + AWH2 = AWH (IR ,JW,JH,IC) + ARWH2 = ARWH(IR ,JW,JH,IC) +C + DA1 = DRL*AR1 - A2 + A1 + DA2 = DRL*AR2 - A2 + A1 + DAW1 = DRL*ARW1 - AW2 + AW1 + DAW2 = DRL*ARW2 - AW2 + AW1 + DAH1 = DRL*ARH1 - AH2 + AH1 + DAH2 = DRL*ARH2 - AH2 + AH1 + DAWH1 = DRL*ARWH1 - AWH2 + AWH1 + DAWH2 = DRL*ARWH2 - AWH2 + AWH1 +C +C-------- set ALFI, dALFI/dWL, dALFI/dHL, d2ALFI/dHLdWL + B(KW,KH) = (1.0-TR)* A1 + TR* A2 + & + ((1.0-TR)*DA1 - TR*DA2 )*(TR-TR*TR) + BW(KW,KH) = (1.0-TR)* AW1 + TR* AW2 + & + ((1.0-TR)*DAW1 - TR*DAW2 )*(TR-TR*TR) + BH(KW,KH) = (1.0-TR)* AH1 + TR* AH2 + & + ((1.0-TR)*DAH1 - TR*DAH2 )*(TR-TR*TR) + BWH(KW,KH) = (1.0-TR)* AWH1 + TR* AWH2 + & + ((1.0-TR)*DAWH1 - TR*DAWH2)*(TR-TR*TR) +C +C-------- also, the RL derivatives of the quantities above + BR(KW,KH) = (A2 - A1 + & + (1.0-4.0*TR+3.0*TR*TR)*DA1 + (3.0*TR-2.0)*TR*DA2 )/DRL + BRW(KW,KH) = (AW2 - AW1 + & + (1.0-4.0*TR+3.0*TR*TR)*DAW1 + (3.0*TR-2.0)*TR*DAW2 )/DRL + BRH(KW,KH) = (AH2 - AH1 + & + (1.0-4.0*TR+3.0*TR*TR)*DAH1 + (3.0*TR-2.0)*TR*DAH2 )/DRL + BRWH(KW,KH) = (AWH2 - AWH1 + & + (1.0-4.0*TR+3.0*TR*TR)*DAWH1 + (3.0*TR-2.0)*TR*DAWH2)/DRL +C + 205 CONTINUE + 20 CONTINUE +C +C---- evaluate spline in HL at the two WL-interval endpoints + DO 30 KW=1, 2 + B1 = B (KW,1) + BR1 = BR (KW,1) + BW1 = BW (KW,1) + BH1 = BH (KW,1) + BRW1 = BRW (KW,1) + BRH1 = BRH (KW,1) + BWH1 = BWH (KW,1) + BRWH1 = BRWH(KW,1) +C + B2 = B (KW,2) + BR2 = BR (KW,2) + BW2 = BW (KW,2) + BH2 = BH (KW,2) + BRW2 = BRW (KW,2) + BRH2 = BRH (KW,2) + BWH2 = BWH (KW,2) + BRWH2 = BRWH(KW,2) +C + DB1 = DHL*BH1 - B2 + B1 + DB2 = DHL*BH2 - B2 + B1 + DBR1 = DHL*BRH1 - BR2 + BR1 + DBR2 = DHL*BRH2 - BR2 + BR1 + DBW1 = DHL*BWH1 - BW2 + BW1 + DBW2 = DHL*BWH2 - BW2 + BW1 + DBRW1 = DHL*BRWH1 - BRW2 + BRW1 + DBRW2 = DHL*BRWH2 - BRW2 + BRW1 +C +C------ set ALFI, dALFI/dRL, dALFI/dWL + C(KW) = (1.0-TH)* B1 + TH* B2 + & + ((1.0-TH)*DB1 - TH*DB2 )*(TH-TH*TH) + CR(KW) = (1.0-TH)* BR1 + TH* BR2 + & + ((1.0-TH)*DBR1 - TH*DBR2 )*(TH-TH*TH) + CW(KW) = (1.0-TH)* BW1 + TH* BW2 + & + ((1.0-TH)*DBW1 - TH*DBW2 )*(TH-TH*TH) + CRW(KW) = (1.0-TH)* BRW1 + TH* BRW2 + & + ((1.0-TH)*DBRW1 - TH*DBRW2)*(TH-TH*TH) +C +C------ also, the HL derivatives of the quantities above + CH(KW) = (B2 - B1 + & + (1.0-4.0*TH+3.0*TH*TH)*DB1 + (3.0*TH-2.0)*TH*DB2 )/DHL + CRH(KW) = (BR2 - BR1 + & + (1.0-4.0*TH+3.0*TH*TH)*DBR1 + (3.0*TH-2.0)*TH*DBR2 )/DHL + CWH(KW) = (BW2 - BW1 + & + (1.0-4.0*TH+3.0*TH*TH)*DBW1 + (3.0*TH-2.0)*TH*DBW2 )/DHL + CRWH(KW) = (BRW2 - BRW1 + & + (1.0-4.0*TH+3.0*TH*TH)*DBRW1 + (3.0*TH-2.0)*TH*DBRW2)/DHL +C + 30 CONTINUE +C +C---- evaluate cubic in WL + C1 = C (1) + CR1 = CR (1) + CW1 = CW (1) + CH1 = CH (1) + CRW1 = CRW (1) + CRH1 = CRH (1) + CWH1 = CWH (1) + CRWH1 = CRWH(1) +C + C2 = C (2) + CR2 = CR (2) + CW2 = CW (2) + CH2 = CH (2) + CRW2 = CRW (2) + CRH2 = CRH (2) + CWH2 = CWH (2) + CRWH2 = CRWH(2) +C + DC1 = DWL*CW1 - C2 + C1 + DC2 = DWL*CW2 - C2 + C1 + DCH1 = DWL*CWH1 - CH2 + CH1 + DCH2 = DWL*CWH2 - CH2 + CH1 + DCR1 = DWL*CRW1 - CR2 + CR1 + DCR2 = DWL*CRW2 - CR2 + CR1 +CC DCRH1 = DWL*CRWH1 - CRH2 + CRH1 +CC DCRH2 = DWL*CRWH2 - CRH2 + CRH1 +C +C---- set AINT, dAINT/dRL, dAINT/dHL + AINT(IC) = (1.0-TW)* C1 + TW* C2 + & + ((1.0-TW)*DC1 - TW*DC2 )*(TW-TW*TW) + AINT_RL = (1.0-TW)* CR1 + TW* CR2 + & + ((1.0-TW)*DCR1 - TW*DCR2 )*(TW-TW*TW) + AINT_HL = (1.0-TW)* CH1 + TW* CH2 + & + ((1.0-TW)*DCH1 - TW*DCH2 )*(TW-TW*TW) +C +C---- also, the WL derivatives of the quantities above + AINT_WL = (C2 - C1 + & + (1.0-4.0*TW+3.0*TW*TW)*DC1 + (3.0*TW-2.0)*TW*DC2 )/DWL + AINTW_RL = (CR2 - CR1 + & + (1.0-4.0*TW+3.0*TW*TW)*DCR1 + (3.0*TW-2.0)*TW*DCR2 )/DWL + AINTW_HL = (CH2 - CH1 + & + (1.0-4.0*TW+3.0*TW*TW)*DCH1 + (3.0*TW-2.0)*TW*DCH2 )/DWL +C + AINTW_WL = ((6.0*TW-4.0)*DC1 + (6.0*TW-2.0)*DC2 )/DWL**2 +C +C +C---- convert derivatives wrt to spline coordinates (RL,WL,HL) into +C- derivatives wrt input variables (Rtheta,f,H) + AINT_R(IC) = (AINT_RL + 0.5*AINT_WL) / (AL10 * RSP) + AINT_W(IC) = (AINT_WL ) / (AL10 * WSP) + AINT_H(IC) = AINT_HL +C + AINTW_R(IC) = (AINTW_RL + 0.5*AINTW_WL) / (AL10**2 * WSP*RSP) + AINTW_W(IC) = (AINTW_WL - AL10*AINT_WL) / (AL10**2 * WSP*WSP) + AINTW_H(IC) = AINTW_HL / (AL10 * WSP ) +C + 1000 CONTINUE +C + ALFR = AINT(1) + ALFR_R = AINT_R(1) + ALFR_W = AINT_W(1) + ALFR_H = AINT_H(1) + ALFRW_R = AINTW_R(1) + ALFRW_W = AINTW_W(1) + ALFRW_H = AINTW_H(1) +C + ALFI = AINT(2) + ALFI_R = AINT_R(2) + ALFI_W = AINT_W(2) + ALFI_H = AINT_H(2) + ALFIW_R = AINTW_R(2) + ALFIW_W = AINTW_W(2) + ALFIW_H = AINTW_H(2) +C +C---- if we're within the spline data space, the derivatives are valid + IF(OK) RETURN +C +C---- if not, the ai value is clamped, and its derivatives are zero + ALFR_R = 0.0 + ALFR_W = 0.0 + ALFR_H = 0.0 + ALFRW_R = 0.0 + ALFRW_W = 0.0 + ALFRW_H = 0.0 +C + ALFI_R = 0.0 + ALFI_W = 0.0 + ALFI_H = 0.0 + ALFIW_R = 0.0 + ALFIW_W = 0.0 + ALFIW_H = 0.0 +C + RETURN +C +C-------------------------------------------------------- +C---- pick up here if OS file not given + 800 CONTINUE + WRITE(*,*)'OSMAP: Environment variable OSMAP not defined' + WRITE(*,*)' Must be set to Orr-Sommerfeld database filename' +C---- don't try again + NOFILE = .TRUE. + RETURN +C +C-------------------------------------------------------- +C---- pick up here for file open error + 900 CONTINUE + WRITE(*,*) + WRITE(*,*)'OSMAP: Orr-Sommerfeld database file not found: ', + & OSFILE(1:LOSF) + WRITE(*,*)' Will return zero amplification rates' +C---- don't try again + NOFILE = .TRUE. + OK = .FALSE. +C + RETURN + END ! OSMAP + -- cgit v1.2.3