PROGRAM RXG DEMO C This program demonstrates the use of the coefficients calculated C by the programs axitrx. The program will produce a table of values C for different values of lambda and s as set in the DATA statements. INTEGER MAXS,MAXL PARAMETER(MAXS=5,MAXL=6) DOUBLE PRECISION DLAMDA(MAXL),S(MAXS),RXG,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 XG11'/) 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)=RXG(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 XG12'/) WRITE(6,20) DO 150 I=1,MAXS DO 140 J=1,MAXL 140 VAL(J)=RXG(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 XG21'/) WRITE(6,20) DO 250 I=1,MAXS DO 240 J=1,MAXL 240 VAL(J)=RXG(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 XG22'/) WRITE(6,20) DO 350 I=1,MAXS DO 340 J=1,MAXL 340 VAL(J)=RXG(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 DOUBLE PRECISION FUNCTION RXG(IALPHA,IBETA,S,DLAMDA) C -------------------- Arguments of function --------------------------- DOUBLE PRECISION DLAMDA,S INTEGER IALPHA,IBETA C ---------------------------------------------------------------------- C G 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 axitrx or axistx C are summed using equations (20a,20b). C ----------------------- Local variables ---------------------------- DOUBLE PRECISION RXGCOF(0:45450), G1,G2,G3,G1SUM,G2SUM,G3SUM DOUBLE PRECISION ZL,COEF, XL,XLP1,XL13,EPS,T,ZT INTEGER I,J,ISTEP,IQ,MMIN,MMAX,M,NM,ISIGN, Nterms 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,Nterms 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 T=2D0/S EPS=0.25D0*S*S-1D0 IF(I.NE.J) THEN ZT=T*T C Notice that because the first coef is a signature we supply the 0. RXG=0D0+G1/EPS-(G2+G3*EPS)*DLOG( 1D0-4D0/(S*S) )-G3 MMIN=2 MMAX=Nterms ELSE ZT=T RXG=G1*S/(2D0*EPS)+(G2+G3*EPS)*DLOG((S+2D0)/(S-2D0))-G3*S MMIN=1 MMAX=Nterms-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 ZL=XLP1**(-M) COEF=0D0 C Some computers will underflow in the next loop. DO 50 IQ=0,M COEF=COEF+RXGCOF(NM)*ZL ZL=ZL*XL NM=NM+ISTEP 50 CONTINUE RXG=RXG+ZT*(COEF-G1SUM-(G2SUM-G3SUM/DBLE(M+2))/DBLE(M) ) ZT=ZT*T 100 CONTINUE RXG = DBLE(ISIGN)*RXG IF(I.EQ.J)RETURN RXG=-RXG*4D0/XLP1**2 IF(I.EQ.2)RXG=RXG*XL**2 RETURN END