PROGRAM RQX 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),RQX,RQX11,RQX12,RQX21,RQX22 INTEGER J DATA DLAMDA/.5D0,1D0,2D0,3D0,4D0,5D0,10D0,20D0,100D0/ WRITE(6,10) 10 FORMAT(1X,' LAMBDA RQX11 RQX12 RQX21 RQX22'/) DO 50 J=1,9 RQX11=RQX(1,1,DLAMDA(J)) RQX12=RQX(1,2,DLAMDA(J)) RQX21=RQX(2,1,DLAMDA(J)) RQX22=RQX(2,2,DLAMDA(J)) 50 WRITE(6,101)DLAMDA(J),RQX11,RQX12,RQX21,RQX22 101 FORMAT(1X,F8.3, 4(F10.5) ) STOP END DOUBLE PRECISION FUNCTION RQX(IALPHA,IBETA,DLAMDA) C ------------------------- Arguments of function ----------------------- INTEGER IALPHA,IBETA DOUBLE PRECISION DLAMDA C ----------------------------------------------------------------------- C X C This subprogram calculates the functions Q (lambda) C ALPHA BETA C defined in Jeffrey,D.J.'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 The subprogram works generally the same as rxq.f where further C comments can be found. C ------------------------ Local variables ------------------------------ DOUBLE PRECISION RXQCOF(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,MAXS LOGICAL LOADED CHARACTER*4 ID DATA LOADED/.FALSE./ C -------------------- Start of executable code ------------------------- IF(.NOT.LOADED) THEN OPEN(1,FILE='rxq300.dat',STATUS='OLD') READ(1,10) ID,MAXS 10 FORMAT(A4,I5) IF(ID.NE.' RXQ') THEN WRITE(*,20) 20 FORMAT(1X,'WRONG COEFFICIENTS FOR RXQ') STOP ENDIF READ(1,30)RXQCOF 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 RQX=1D0 - G1/4D0 - G3 MMIN=2 MMAX=200 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 RQX=G1/4D0+G2*DLOG(4D0)-2D0*G3 MMIN=1 MMAX=199 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+RXQCOF(NM)*ZL NM=NM+ISTEP ZL=ZL*XL 50 CONTINUE COEF=COEF*XLP1**(-M) RQX=RQX+COEF-G1SUM - (G2SUM-G3SUM/DBLE(M+2))/DBLE(M) 100 CONTINUE IF(I.EQ.J)RETURN RQX=RQX*8D0/XLP1**3 IF(I.EQ.2)RQX=RQX*XL**3 RETURN END C BLOCK DATA RQXINI DOUBLE PRECISION RXQCOF COMMON /XQD200/RXQCOF(0:45450) DATA RXQCOF(1)/0D0/ END