PROGRAM RHY DEM C This program demonstrates the use of the coefficients calculated C by the programs trnrox or trnstx. 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),RHY,RHY11,RHY12,RHY21,RHY22 INTEGER J DATA DLAMDA/.5D0,1D0,2D0,3D0,4D0,5D0,10D0,20D0,100D0/ WRITE(6,5) 5 FORMAT(1X, &'RESISTANCE FUNCTIONS RHY from Jeffrey, Phys. Fluids, 1992') WRITE(6,7) 7 FORMAT(1X,'vol A4, p16.') WRITE(6,8) 8 FORMAT(/1X,'Tabulation of the RHY appearing in eqs. (35),(37)') WRITE(6,10) 10 FORMAT(1X,' LAMBDA RHY11 RHY12 RHY21 RHY22'/) DO 50 J=1,9 RHY11=RHY(1,1,DLAMDA(J)) RHY12=RHY(1,2,DLAMDA(J)) RHY21=RHY(2,1,DLAMDA(J)) RHY22=RHY(2,2,DLAMDA(J)) 50 WRITE(6,101)DLAMDA(J),RHY11,RHY12,RHY21,RHY22 101 FORMAT(1X,F8.3,4F10.5) STOP END DOUBLE PRECISION FUNCTION RHY(IALPHA,IBETA,DLAMDA) C ----------------------- Arguments of function ----------------------- DOUBLE PRECISION DLAMDA INTEGER IALPHA, IBETA C --------------------------------------------------------------------- C This calculates the functions RHY defined in Jeffrey, Phys. Fluids (1992), C 'The calculation of the low Reynolds number resistance C functions for two unequal spheres'. The coefficients are read from the C file ryh300.dat and are summed using equations (37). C 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 RYHCOF(0:45450) DOUBLE PRECISION XL,XLP1,G2,G3,G2SUM,G3SUM DOUBLE PRECISION ZL,COEF 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(4,FILE='ryh300.dat',STATUS='OLD') READ(4,10) ID, MAXS 10 FORMAT( A4,I5) IF( ID.NE.' RYH') THEN WRITE(*,20) 20 FORMAT( 1X,'WRONG COEFFICIENTS FOR RHY') STOP ENDIF READ(4,2000)RYHCOF 2000 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 C These are the g_i for Y^H_{11,22} (i=1,2 in paper) G2=(2D0*XL-XL**2)/(10D0*XLP1**2) G3=(16D0-61D0*XL+180D0*XL**2+2D0*XL**3)/(500D0*XLP1**2) ISTEP=1 ELSE G2=(2D0*XL-1D0)/(10D0*XLP1**2) G3=(16D0*XL**3-61D0*XL**2+180D0*XL+2D0)/(500D0*XL*XLP1**2) ISTEP=-1 ENDIF RHY=-G3 MMIN=2 MMAX=MAXS ELSE C These are the g_i for Y^H_{12,21} (i=5,6 in paper) IF(I.EQ.1) THEN G2=(XL**2+7D0*XL**3)/(20D0*XLP1**2) G3=(43D0*XL+147D0*XL**2-185D0*XL**3+221D0*XL**4)/ & (1000D0*XLP1**2) ISTEP=1 ELSE G2=(7D0*XL+XL**2)/(20D0*XL**2*XLP1**2) G3=(221D0-185D0*XL+147D0*XL**2+43D0*XL**3)/ & (1000D0*XL**2*XLP1**2) ISTEP=-1 ENDIF RHY=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+RYHCOF(NM)*ZL NM=NM+ISTEP ZL=ZL*XL 50 CONTINUE RHY=RHY+COEF-(G2SUM-G3SUM/DBLE(M+2))/DBLE(M) 100 CONTINUE IF(I.EQ.J)RETURN RHY=RHY*8D0/XLP1**3 IF(I.EQ.2)RHY=RHY*XL**3 RETURN END