PROGRAM RGX DEMO C This program demonstrates the use of the coefficients calculated C by the programs axitrx or axistx. The program will produce a table of C values for different values of lambda as set in the DATA statements. C DOUBLE PRECISION DLAMDA(9),RGX,RGX11,RGX12,RGX21,RGX22 INTEGER J DATA DLAMDA/.5D0,1D0,2D0,3D0,4D0,5D0,10D0,20D0,100D0/ WRITE(6,5) 5 FORMAT(1X, &'RESISTANCE FUNCTIONS RGX from Jeffrey, 1992, Phys. Fluids, &A4, 16.') WRITE(6,7) 7 FORMAT(1X, &'Tabulation of the functions RGX appearing in eqs. (19),(21).') WRITE(6,10) 10 FORMAT(1X,' LAMBDA RGX11 RGX12 RGX21 RGX22'/) DO 50 J=1,9 RGX11=RGX(1,1,DLAMDA(J)) RGX12=RGX(1,2,DLAMDA(J)) RGX21=RGX(2,1,DLAMDA(J)) RGX22=RGX(2,2,DLAMDA(J)) 50 WRITE(6,101)DLAMDA(J),RGX11,RGX12,RGX21,RGX22 101 FORMAT(1X,F8.3, 4(F10.5) ) STOP END DOUBLE PRECISION FUNCTION RGX(IALPHA,IBETA,DLAMDA) C ------------------ Arguments of function --------------------------- DOUBLE PRECISION DLAMDA INTEGER IALPHA,IBETA C -------------------------------------------------------------------- C X C This subprogram calculates the functions G defined in C ALPHA BETA C Jeffrey, 'The calculation of the low Reynolds number resistance C functions for two unequal spheres', Phys. Fluids (1992). C Coefficients calculated by one of the programs axitrx or axistx are C summed using equations (21). C C ----------------------- Local variables ---------------------------- DOUBLE PRECISION RXGCOF(0:45450), G1,G2,G3,G1SUM,G2SUM,G3SUM DOUBLE PRECISION ZL,COEF, XL,XLP1,XL13 INTEGER I,J,ISTEP,IQ,MMIN,MMAX,M,NM,ISIGN, MAXS LOGICAL LOADED CHARACTER*4 ID DATA LOADED/.FALSE./ C ---------------------- Start of executable code -------------------- IF(.NOT.LOADED) THEN OPEN(1,FILE='rxg300.dat',STATUS='OLD') READ(1,10) ID,MAXS 10 FORMAT(A4,I5) IF(ID.NE.' RXG') THEN WRITE(*,20) 20 FORMAT(1X,'WRONG COEFFICIENTS FOR RXG') STOP ENDIF READ(1,30) RXGCOF 30 FORMAT(D22.16) CLOSE(1) LOADED=.TRUE. ENDIF IF(DLAMDA.GT.1D0) THEN I=3-IALPHA J=3-IBETA XL=1D0/DLAMDA C RXG and RGX change sign when we invert lambda. ISIGN = -1 ELSE I=IALPHA J=IBETA XL=DLAMDA ISIGN = 1 ENDIF XLP1=XL+1D0 XL13=XLP1**3 C These are the g_i for X^G_{11,12} (i=1..3 in paper) IF(I.EQ.1) THEN G1=3D0*XL**2/XL13 G2=3D0*(XL+12D0*XL**2-4D0*XL**3)/(10D0*XL13) G3=(5D0+181D0*XL-453D0*XL**2+566D0*XL**3-65D0*XL**4)/ & (140D0*XL13) ISTEP=1 ELSE C These are the g_i for X^G_{21,22} G1=3D0*XL/XL13 G2=3D0*(-4D0+12D0*XL+XL*XL)/(10D0*XL13) G3=(-65D0+566D0*XL-453D0*XL**2+181D0*XL**3+5D0*XL**4)/ & (140D0*XL*XL13) C These expressions essentially invert lambda. ISIGN=-ISIGN ISTEP=-1 ENDIF IF(I.NE.J) THEN RGX=0D0-G1/4D0 -G3 MMIN=2 MMAX=MAXS ELSE RGX=G1/4D0+G2*DLOG(4D0)-2D0*G3 MMIN=1 MMAX=MAXS-1 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=XLP1**(-M) COEF=0D0 C Some computers will underflow in the next loop. DO 50 IQ=0,M COEF=COEF+RXGCOF(NM)*ZL NM=NM+ISTEP ZL=ZL*XL 50 CONTINUE RGX=RGX+COEF-G1SUM-(G2SUM-G3SUM/DBLE(M+2))/DBLE(M) 100 CONTINUE RGX = DBLE(ISIGN)*RGX IF(I.EQ.J)RETURN RGX=-RGX*4D0/XLP1**2 IF(I.EQ.2)RGX=RGX*XL**2 RETURN END