PROGRAM RXM DEM C This program demonstrates the use of the coefficients calculated C by the programs axistx. The program will produce a table of values C for different values of LAMBDA and S as set in the DATA statements. C INTEGER MAXS,MAXL PARAMETER(MAXS=5,MAXL=6) DOUBLE PRECISION DLAMDA(MAXL),S(MAXS),RXM,VAL(MAXL) INTEGER I,J DATA DLAMDA/1D0,2D0,3D0,4D0,5D0,0.5D0/ DATA S/2.01D0,2.02D0,2.05D0,2.1D0,2.5D0/ WRITE(6,10) 10 FORMAT(1X,'TABLE OF VALUES FOR RESISTANCE FUNCTION XM11'/) WRITE(6,20) 20 FORMAT(' s',6X,25('-'),' lambda ',24('-'),/, 14X, & '1 2 3 4 5 1/2') DO 50 I=1,MAXS DO 40 J=1,MAXL 40 VAL(J)=RXM(1,1,S(I),DLAMDA(J)) 50 WRITE(6,1100)S(I),(VAL(J),J=1,MAXL) WRITE(6,110) 110 FORMAT(/1X,'TABLE OF VALUES FOR RESISTANCE FUNCTION XM12'/) WRITE(6,20) DO 150 I=1,MAXS DO 140 J=1,MAXL 140 VAL(J)=RXM(1,2,S(I),DLAMDA(J)) 150 WRITE(6,1100)S(I),(VAL(J),J=1,MAXL) WRITE(6,210) 210 FORMAT(/1X,'TABLE OF VALUES FOR RESISTANCE FUNCTION XM21'/) WRITE(6,20) DO 250 I=1,MAXS DO 240 J=1,MAXL 240 VAL(J)=RXM(2,1,S(I),DLAMDA(J)) 250 WRITE(6,1100)S(I),(VAL(J),J=1,MAXL) WRITE(6,310) 310 FORMAT(/1X,'TABLE OF VALUES FOR RESISTANCE FUNCTION XM22'/) WRITE(6,20) DO 350 I=1,MAXS DO 340 J=1,MAXL 340 VAL(J)=RXM(2,2,S(I),DLAMDA(J)) 350 WRITE(6,1100)S(I),(VAL(J),J=1,MAXL) 1100 FORMAT(F8.3,8F10.4) STOP END C DOUBLE PRECISION FUNCTION RXM(IALPHA,IBETA,S,DLAMDA) C ---------------------- Arguments of function ------------------------ INTEGER IALPHA, IBETA DOUBLE PRECISION S, DLAMDA C --------------------------------------------------------------------- C M C This routine calculates the functions X (s,lambda) C ALPHA BETA C defined in Jeffrey, 'The calculation of the low-Reynolds-number C resistance functions for two unequal spheres',Phys. Fluids (1991). C Coefficients calculated by one of the programs AXISTx are summed. C It goes forwards or backwards through the coefficients ( ISTEP) C depending upon whether I=1 or 2. C This technique is not explained in the paper. C Starting at equn (5.d) we see that X can be obtained from X C 22 11 C by inverting lambda. Now looking at equation (47) we substitute C 1/lambda and redefine q as k-q. Factoring the lambda**k, we see that C lambda**q now multiplies P(1,q,k-q), which means running 'backwards' C through the coefficients relative to the published case. C C The program reads the coefficients from RXM200.DAT only once and checks C to make certain the signature is correct (cf AXIST1). C C ------------------------ Local variables ---------------------------- DOUBLE PRECISION RXMCOF(0:45450) DOUBLE PRECISION G1, G2, G3, G1SUM, G2SUM, G3SUM DOUBLE PRECISION T, ZT, ZL, COEF, EPS, XL, XLP1, XL13 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='rxm300.dat',STATUS='OLD') READ(1,10) ID,MAXS 10 FORMAT(A4,I5) IF (ID.NE.' RXM') THEN WRITE(*,20) 20 FORMAT(1X,'WRONG COEFFICIENTS FOR RXM') STOP ENDIF READ(1,30) RXMCOF 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 XL13=XLP1**3 IF(I.EQ.J) THEN IF(I.EQ.1) THEN C These are the g_i for X^M_{11} (i=1..3 in paper) G1=6D0*XL**2/(5D0*XL13) G2=3D0*(XL+17D0*XL**2-9D0*XL**3)/(25D0*XL13) G3=(5D0+272D0*XL-831D0*XL**2+1322D0*XL**3-415D0*XL**4) & /(350D0*XL13) ISTEP=1 ELSE C These are the g_i for X^M_{22} G1=6D0*XL/(5D0*XL13) G2=3D0*(XL**2+17D0*XL-9D0)/(25D0*XL13) G3=(5D0*XL**4+272D0*XL**3-831D0*XL**2+1322D0*XL-415D0) & /(350D0*XL*XL13) ISTEP=-1 ENDIF ELSE IF(I.EQ.1) THEN C These are the g_i for X^M_{12} (i=4..6 in paper) G1=6D0*XL**3/(5D0*XL13) G2=3D0*(-4D0*XL**2+17D0*XL**3-4D0*XL**4)/(25D0*XL13) G3=(-65D0*XL+832D0*XL**2-1041D0*XL**3+832D0*XL**4-65D0*XL**5) & /(350D0*XL13) ISTEP=1 ELSE C These are the g_i for X^M_21 G1=6D0/(5D0*XL13) G2=3D0*(-4D0*XL**2+17D0*XL-4D0)/(25D0*XL*XL13) G3=(-65D0+832D0*XL-1041D0*XL**2+832D0*XL**3-65D0*XL**4) & /(350D0*XL*XL*XL13) ISTEP=-1 ENDIF ENDIF T=2D0/S EPS=0.25D0*S*S-1D0 IF(I.EQ.J) THEN ZT=T*T C Notice that because the first coef is a signature we supply the 1. RXM=1D0+G1/EPS+(G2+G3*EPS)*DLOG( (S*S)/(S*S-4D0))-G3 MMIN=2 MMAX=MAXS ELSE ZT=T RXM=G1*S/(2D0*EPS)+(G2+G3*EPS)*DLOG((S+2D0)/(S-2D0))-G3*S MMIN=1 MMAX=MAXS-1 ENDIF T=T*T G1SUM=G1 G2SUM=2D0*G2 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 C This loop adds up the powers of lambda to obtain the term f_m in the C paper. Notice that the 2**(-m) is not needed, because the coefficients C in the file already contain that factor. ZL=1D0 COEF=0D0 C Some computers will underflow in the next loop. DO 50 IQ=0,M COEF=COEF+RXMCOF(NM)*ZL NM=NM+ISTEP ZL=ZL*XL 50 CONTINUE COEF=COEF*XLP1**(-M) RXM=RXM + (COEF-G1SUM-(G2SUM-G3SUM/DBLE(M+2))/DBLE(M) )*ZT ZT=ZT*T 100 CONTINUE IF(I.EQ.J)RETURN RXM=RXM*8D0/XLP1**3 IF(I.EQ.2)RXM=RXM*XL**3 RETURN END C BLOCK DATA RXMINI DOUBLE PRECISION RXMCOF(0:45450) COMMON /XMD200/RXMCOF DATA RXMCOF(0)/0D0/ END