PROGRAM RXC DEM C This program demonstrates the use of the coefficients calculated C by the programs axirox. The program will produce a table of values C for different values of lambda and s as set in the DATA statements. C DOUBLE PRECISION DLAMDA(5),S(4),RXC,RXC11,RXC12,RXC21,RXC22 INTEGER I,J DATA DLAMDA/1D0,.5D0,.25D0,.125D0,2D0/ DATA S/2.01D0,2.1D0,2.5D0,3D0/ WRITE(*,20) 20 FORMAT(1X,'Table of values for Resistance Functions XC') DO 100 J=1,5 WRITE(*,30)DLAMDA(J) 30 FORMAT(/1X,'LAMBDA =',F8.3) WRITE(*,40) 40 FORMAT(/7X,'S',10X, & 'RXC(1,1)',8X,'RXC(1,2)',8X,'RXC(2,1)',8X,'RXC(2,2)') DO 100 I=1,4 RXC11=RXC(1,1,S(I),DLAMDA(J)) RXC12=RXC(1,2,S(I),DLAMDA(J)) RXC21=RXC(2,1,S(I),DLAMDA(J)) RXC22=RXC(2,2,S(I),DLAMDA(J)) 100 WRITE(*,50)S(I),RXC11,RXC12,RXC21,RXC22 50 FORMAT(F10.3,4F16.5) STOP END C DOUBLE PRECISION FUNCTION RXC(IALPHA,IBETA,S,DLAMDA) C ---------------------- Arguments of function -------------------------- DOUBLE PRECISION DLAMDA,S INTEGER IALPHA,IBETA C ----------------------------------------------------------------------- C C C This routine calculates the functions X (s,lambda) C ALPHA BETA C defined in Jeffrey & Onishi (1984) J. Fluid Mech. Vol 139, 261-290. C Coefficients calculated by one of the programs axirox are summed using C equations (6.12)-(6.13). The subprogram uses equation (1.9a) to convert C any call to one in which lambda.LT.1. It goes forwards or backwards C through the coefficients ( ISTEP) depending upon whether I=1 or 2. C This technique is not explained in the paper; I have forgotten why. C Starting at equn 1.9 we see that X can be obtained from X C 22 11 C by inverting lambda. Now looking at equn 3.15 we substitute 1/lambda and C redefine q as k-q. Factoring the lambda**k, we see that lambda**q now C multiplies P(1,q,k-q), which means running 'backwards' through the C coefficients relative to the published case. C C In the program the coefficients are read in only once and checked C to make certain the signature is correct (cf axiro1). C The coefficients are expected in rxc300.dat in the format on the C distribution disk. C C -----------------------------Local Variables--------------------------- DOUBLE PRECISION RXCCOF(0:45450),G3,G3SUM DOUBLE PRECISION T,ZT,ZL,COEF,EPS,XL,XLP1 INTEGER I,J,ISTEP,IQ,MMIN,MMAX,M,NM LOGICAL LOADED CHARACTER*4 ID DATA LOADED/.FALSE./ C ---------------------------- Start of executable code------------------ IF(.NOT.LOADED) THEN OPEN(1,FILE='rxc300.dat', STATUS='OLD') READ(1,10) ID,MAXS 10 FORMAT(A4,I5) IF (ID.NE.' RXC') THEN WRITE(*,20) 20 FORMAT(1X,'WRONG COEFFICIENTS FOR RXC') STOP ENDIF READ(1,30) RXCCOF 30 FORMAT(D22.16) CLOSE(1) LOADED=.TRUE. ENDIF IF(DLAMDA.GT.1D0) THEN I=3-IALPHA J=3-IBETA XL=1D0/DLAMDA ELSE I=IALPHA J=IBETA XL=DLAMDA ENDIF XLP1=XL+1D0 IF(I.EQ.1) THEN G3=-XL*XL/XLP1 ISTEP=1 ELSE G3=-1D0/(XL*XLP1) ISTEP=-1 ENDIF T=2D0/S EPS=1D0-4D0/(S*S) IF(I.EQ.J) THEN ZT=T*T C Notice that because the first coef is a signature we supply the 1. RXC=1D0-(G3*EPS*S*S/4D0)*DLOG(EPS)-G3 MMIN=2 MMAX=MAXS ELSE ZT=T RXC=(G3*EPS*S*S/4D0)*DLOG((S+2D0)/(S-2D0)) RXC=RXC-G3*S MMIN=1 MMAX=MAXS-1 ENDIF T=T*T G3SUM=4D0*G3 DO 100 M=MMIN,MMAX,2 IF (ISTEP.GT.0) THEN NM = (M*(M+1))/2 ELSE NM = (M*(M+1))/2 +M ENDIF ZL=XLP1**(-M) COEF=0D0 DO 50 IQ=0,M COEF=COEF+RXCCOF(NM)*ZL NM=NM+ISTEP ZL=ZL*XL 50 CONTINUE RXC=RXC+ZT*(COEF+(G3SUM/(DBLE(M+2)*DBLE(M)))) ZT=ZT*T 100 CONTINUE IF(I.EQ.J)RETURN RXC=-RXC*8D0/XLP1**3 IF(I.EQ.2)RXC=RXC*XL**3 RETURN END C BLOCK DATA RXCINI DOUBLE PRECISION RXCCOF(0:20300) COMMON /RXC200/RXCCOF DATA RXCCOF(0)/0D0/ END