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