PROGRAM RMY 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 as set in the DATA statements. C DOUBLE PRECISION DLAMDA(9),RMY,RMY11,RMY12,RMY21,RMY22 INTEGER J DATA DLAMDA/.5D0,1D0,2D0,3D0,4D0,5D0,10D0,20D0,100D0/ WRITE(6,5) 5 FORMAT(1X, &'RESISTANCE FUNCTIONS from Jeffrey, Phys. Fluids, 1992' &1X,'vol A4,p16' &/1X,'Tabulation of the RMY appearing in eqs. (64), (66).') WRITE(6,10) 10 FORMAT(1X,' LAMBDA RMY11 RMY12 RMY21 RMY22'/) DO 50 J=1,9 RMY11=RMY(1,1,DLAMDA(J)) RMY12=RMY(1,2,DLAMDA(J)) RMY21=RMY(2,1,DLAMDA(J)) RMY22=RMY(2,2,DLAMDA(J)) 50 WRITE(6,101)DLAMDA(J),RMY11,RMY12,RMY21,RMY22 101 FORMAT(1X,F8.3, 4(F10.5) ) STOP END C DOUBLE PRECISION FUNCTION RMY(IALPHA,IBETA,DLAMDA) C ---------------------- Arguments of function ------------------------- DOUBLE PRECISION DLAMDA INTEGER IALPHA, IBETA C ---------------------------------------------------------------------- C This calculates the function RMY defined in Jeffrey, 'The calculation of C the low-Reynolds-number resistance functions for two unequal spheres', C Phys. Fluids (1991). C The coefficients are summed using equations (66a,66b). C The coefficients must conform to those generated by one of the programs C TRNSTx. The subprogram works by converting any call to one in which C lambda.LT.1. It goes forwards or backwards through the coefficients C (the variable ISTEP) depending upon whether I=1 or 2. C -------------------------- Local Variables --------------------------- DOUBLE PRECISION RYMCOF(0:45450) DOUBLE PRECISION XL,XLP1,G2,G3,G2SUM,G3SUM DOUBLE PRECISION ZL,COEF 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(4,FILE='rym300.dat',STATUS='OLD') READ(4,10) ID,MAXS 10 FORMAT(A4,I5) IF (ID.NE.' RYM') THEN WRITE(*,20) 20 FORMAT(1X,'WRONG COEFFICIENTS FOR RYM') STOP ENDIF READ(4,30) RYMCOF 30 FORMAT(D22.16) CLOSE(4) 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. RMY=1D0-G3 MMIN=2 MMAX=MAXS 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 RMY=2D0*(G2*DLOG(2D0)-G3) MMIN=1 MMAX=MAXS-1 ENDIF 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 DO 50 IQ=0,M COEF=COEF+RYMCOF(NM)*ZL NM=NM+ISTEP ZL=ZL*XL 50 CONTINUE RMY=RMY+COEF-(G2SUM-G3SUM/DBLE(M+2))/DBLE(M) 100 CONTINUE IF(I.EQ.J)RETURN RMY=RMY*8D0/XLP1**3 IF(I.EQ.2)RMY=RMY*XL**3 RETURN END C BLOCK DATA RMYINI DOUBLE PRECISION RYMCOF(0:45450) COMMON /YMD300/RYMCOF DATA RYMCOF(1)/0D0/ END