PROGRAM RYM DEM C This program demonstrates the use of the coefficients calculated C by the programs trnstx. The program will produce a table of values C for different values of LAMBDA and S as set in the DATA statements. C DOUBLE PRECISION DLAMDA(5),S(4),RYM,RYM11,RYM12,RYM21,RYM22 INTEGER I,J DATA DLAMDA/1D0,.5D0,.25D0,.125D0,2D0/ DATA S/2.01D0,2.02D0,2.1D0,2.5D0/ C WRITE(*,20) 20 FORMAT(1X,'TABLE OF VALUES FOR RESISTANCE FUNCTIONS YM') DO 100 J=1,5 WRITE(*,30)DLAMDA(J) 30 FORMAT(/1X,'LAMBDA =',F8.3) WRITE(*,40) 40 FORMAT(/7X,'S',10X, & 'RYM(1,1)',8X,'RYM(1,2)',8X,'RYM(2,1)',8X,'RYM(2,2)') DO 100 I=1,4 RYM11=RYM(1,1,S(I),DLAMDA(J)) RYM12=RYM(1,2,S(I),DLAMDA(J)) RYM21=RYM(2,1,S(I),DLAMDA(J)) RYM22=RYM(2,2,S(I),DLAMDA(J)) WRITE(*,50)S(I),RYM11,RYM12,RYM21,RYM22 50 FORMAT(F10.3,4F16.5) 100 CONTINUE STOP END C DOUBLE PRECISION FUNCTION RYM(IALPHA,IBETA,S,DLAMDA) C ----------------------------- Arguments of function ----------------- INTEGER IALPHA,IBETA DOUBLE PRECISION S,DLAMDA C --------------------------------------------------------------------- C M C This routine calculates the functions Y (s,lambda) C ALPHA BETA C defined in 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 TRNSTx are summed using C equations (65). C The technique used is the same as for the XA functions, defined C in Jeffrey and Onishi (1984) J. Fluid Mech. Vol 139, 261-290. C It goes forwards or backwards through the coefficients ( ISTEP) C depending upon whether I=1 or 2. This technique is not explained in the C paper. C Starting at equn (1.9) (J&O,1984) we see that Y can be obtained from Y C 22 11 C by inverting lambda. Now looking at equn (3.15) (J&O,1984) we substitute C 1/lambda and redefine q as k-q. Factoring the lambda**k, we see C that 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 rym300.dat only once and checks C to make certain the signature is correct (cf axist1). C C ------------------------ Local variables ---------------------------- DOUBLE PRECISION RYMCOF(0:45450) DOUBLE PRECISION G2, G3, G2SUM, G3SUM DOUBLE PRECISION T, ZT, ZL, COEF, EPS, XL, XLP1 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='rym300.dat',STATUS='OLD') READ(1,10) ID,MAXS 10 FORMAT(A4,I5) IF(ID.NE.' RYM') THEN WRITE(*,20) 20 FORMAT(1X,'WRONG COEFFICIENTS FOR RYM') STOP ENDIF READ(1,30) RYMCOF 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 IF(I.EQ.J) THEN IF(I.EQ.1) THEN G2=6D0*(XL-XL**2+4D0*XL**3)/(25D0*XLP1**3) G3=(24D0-201D0*XL+882D0*XL**2-1182D0*XL**3+591D0*XL**4) & /(625D0*XLP1**3) ISTEP=1 ELSE G2=6D0*(4D0-XL+XL**2)/(25D0*XLP1**3) G3=(24D0*XL**4-201D0*XL**3+882D0*XL**2-1182D0*XL+591D0) & /(625D0*XL*XLP1**3) ISTEP=-1 ENDIF C Notice that because the first coef is a signature we supply the 1. ELSE IF(I.EQ.1) THEN G2=3D0*(7D0*XL*XL-10D0*XL**3+7D0*XL**4)/(50D0*XLP1**3) G3=(663D0*XL-2244D0*XL**2+5706D0*XL**3-2244D0*XL**4 & +663D0*XL**5)/(2500D0*XLP1**3) ISTEP=1 ELSE G2=3D0*(7D0*XL**4-10D0*XL**3+7D0*XL*XL)/(50D0*XL**3*XLP1**3) G3=(663D0*XL**5-2244D0*XL**4+5706D0*XL**3-2244D0*XL**2 & +663D0*XL)/(2500D0*XLP1**3*XL**3) 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. RYM=1D0 - (G2+G3*EPS)*DLOG( 1D0-4D0/(S*S) ) - G3 MMIN=2 MMAX=MAXS ELSE ZT=T RYM=(G2+G3*EPS)*DLOG((S+2D0)/(S-2D0))-G3*S MMIN=1 MMAX=MAXS-1 ENDIF T=T*T 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 on the next statements. DO 50 IQ=0,M COEF=COEF+RYMCOF(NM)*ZL NM=NM+ISTEP ZL=ZL*XL 50 CONTINUE COEF=COEF*XLP1**(-M) RYM=RYM + (COEF-(G2SUM-G3SUM/DBLE(M+2))/DBLE(M) )*ZT ZT=ZT*T 100 CONTINUE IF(I.EQ.J)RETURN RYM=RYM*8D0/XLP1**3 IF(I.EQ.2)RYM=RYM*XL**3 RETURN END C BLOCK DATA RYMINI DOUBLE PRECISION RYMCOF(0:45450) COMMON /YMD200/RYMCOF DATA RYMCOF(1)/0D0/ END