PROGRAM RMX 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 as set in the DATA statements. C DOUBLE PRECISION DLAMDA(9),RMX,RMX11,RMX12,RMX21,RMX22 INTEGER J DATA DLAMDA/.5D0,1D0,2D0,3D0,4D0,5D0,10D0,20D0,100D0/ WRITE(6,5) 5 FORMAT(1X,'RESISTANCE FUNCTIONS RMX from Jeffrey, 1992,') WRITE(6,6) 6 FORMAT(1X,'Phys. Fluids, A4, 16.') WRITE(6,7) 7 FORMAT(1X, &'Tabulation of the RMX appearing in eqs. (48),(50).') WRITE(6,10) 10 FORMAT(1X,' LAMBDA RMX11 RMX12 RMX21 RMX22'/) DO 50 J=1,9 RMX11=RMX(1,1,DLAMDA(J)) RMX12=RMX(1,2,DLAMDA(J)) RMX21=RMX(2,1,DLAMDA(J)) RMX22=RMX(2,2,DLAMDA(J)) 50 WRITE(6,101)DLAMDA(J),RMX11,RMX12,RMX21,RMX22 101 FORMAT(1X,F8.3, 4(F10.5) ) STOP END DOUBLE PRECISION FUNCTION RMX(IALPHA,IBETA,DLAMDA) C ------------------------- Arguments of function ----------------------- INTEGER IALPHA,IBETA DOUBLE PRECISION DLAMDA C ----------------------------------------------------------------------- C X C This subprogram calculates the functions M (lambda) C ALPHA BETA C defined in Jeffrey, 'The calculation of the low Reynolds number C resistance functions for two unequal spheres', Phys. Fluids (1992). C Coefficients are summed using equations (50). C C -------------------------- Local variables ---------------------------- DOUBLE PRECISION RXMCOF(0:45450) DOUBLE PRECISION XL,XLP1,G1,G2,G3,G1SUM,G2SUM,G3SUM DOUBLE PRECISION ZL, COEF, 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 RMX=1D0 - G1/4D0 - G3 MMIN=2 MMAX=300 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 RMX=G1/4D0+G2*DLOG(4D0)-2D0*G3 MMIN=1 MMAX=299 ENDIF 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 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) RMX=RMX+COEF-G1SUM - (G2SUM-G3SUM/DBLE(M+2))/DBLE(M) 100 CONTINUE IF(I.EQ.J)RETURN RMX=RMX*8D0/XLP1**3 IF(I.EQ.2)RMX=RMX*XL**3 RETURN END C BLOCK DATA RMXINI DOUBLE PRECISION RXMCOF COMMON /XMD300/RXMCOF(0:45450) DATA RXMCOF(1)/0D0/ END