PROGRAM RYG DEM
      INTEGER MAXS,MAXL
      PARAMETER(MAXL=6,MAXS=5)
C This program demonstrates the use of the coefficients calculated
C  by the programs trntrx or trnstx. The program will produce a table of
C  values for different values of lambda and s as set in the DATA statements.
C
      DOUBLE PRECISION DLAMDA(MAXL),S(MAXS),RYG,VAL(MAXL)
      INTEGER I,J
      DATA DLAMDA/1D0,2D0,3D0,4D0,5D0,0.5D0/
      DATA S/2.01D0,2.02D0,2.05D0,2.1D0,2.5D0/
      WRITE(6,10)
10    FORMAT(1X,'TABLE OF VALUES FOR RESISTANCE FUNCTION YG11'/
     & '   s        l=1       l=2        l=3        l=4       l=5')
      DO 50 I=1,MAXS
        DO 40 J=1,MAXL
40      VAL(J)=RYG(1,1,S(I),DLAMDA(J))
50    WRITE(6,1100)S(I),(VAL(J),J=1,MAXL)
      WRITE(6,110)
110    FORMAT(/1X,'TABLE OF VALUES FOR RESISTANCE FUNCTION YG12'/
     & '   s        l=1       l=2        l=3        l=4       l=5')
      DO 150 I=1,MAXS
        DO 140 J=1,MAXL
140     VAL(J)=RYG(1,2,S(I),DLAMDA(J))
150   WRITE(6,1100)S(I),(VAL(J),J=1,MAXL)
      WRITE(6,210)
210    FORMAT(/1X,'TABLE OF VALUES FOR RESISTANCE FUNCTION YG21'/
     & '   s        l=1       l=2        l=3        l=4       l=5')
      DO 250 I=1,MAXS
        DO 240 J=1,MAXL
240     VAL(J)=RYG(2,1,S(I),DLAMDA(J))
250   WRITE(6,1100)S(I),(VAL(J),J=1,MAXL)
      WRITE(6,310)
310    FORMAT(/1X,'TABLE OF VALUES FOR RESISTANCE FUNCTION YG22'/
     & '   s        l=1       l=2        l=3        l=4       l=5')
      DO 350 I=1,MAXS
        DO 340 J=1,MAXL
340     VAL(J)=RYG(2,2,S(I),DLAMDA(J))
350   WRITE(6,1100)S(I),(VAL(J),J=1,MAXL)
 1100 FORMAT(F8.3,8F10.4)
      STOP
      END
      DOUBLE PRECISION FUNCTION RYG(IALPHA,IBETA,S,DLAMDA)
C ------------------------- Arguments of function ---------------------
      DOUBLE PRECISION DLAMDA,S
      INTEGER IALPHA,IBETA
C ---------------------------------------------------------------------
C This calculates the function RYG defined in Jeffrey, 'The calculation
C  of the low Reynolds number resistance functions for two unequal spheres',
C  Phys. Fluids (1992) .
C  Coefficients calculated by one of the programs TRNTRx or TRNSTx are summed
C  using equations (28a,28b). The subprogram works by converting any call
C  to one in which lambda.LT.1. It goes forwards or backwards through
C  the coefficients (the variable ISTEP) depending upon whether I=1 or 2.
C  This version was written during August 1988. Changes 1997.
C ------------------------  Local variables -----------------------------
      DOUBLE PRECISION RYGCOF(0:45450)
      DOUBLE PRECISION XL,XLP1,G2,G3,G2SUM,G3SUM
      DOUBLE PRECISION T,ZT,ZL,COEF
      INTEGER I,J,ISTEP,IQ,MMIN,MMAX,M,NM,ISIGN
      LOGICAL LOADED
      CHARACTER*4 ID
      DATA LOADED/.FALSE./
C -------------------- Start of executable code -------------------------
      IF( .NOT.LOADED) THEN
        OPEN(3,FILE='ryg300.dat',STATUS='OLD')
        READ(3,10) ID,MAXS
10     FORMAT(A4,I5)
        IF( ID.NE.' RYG') THEN
           WRITE(*,20)
20         FORMAT( 1X,'WRONG COEFFICIENTS FOR RYG')
            STOP
        ENDIF
        READ(3,30)RYGCOF
 30     FORMAT(D22.16)
        CLOSE(3)
        LOADED=.TRUE.
      ENDIF
      IF(DLAMDA.GT.1D0) THEN
        I=3-IALPHA
        J=3-IBETA
        XL=1D0/DLAMDA
C RYG and RGY change sign when we invert lambda
        ISIGN = -1
      ELSE
        I=IALPHA
        J=IBETA
        XL=DLAMDA
        ISIGN = 1
      ENDIF
      XLP1=XL+1D0
C
C  These Gx are the functions g<x> from Jeffrey (1991) (i=1,2 in paper)
      IF(I.EQ.1) THEN
        G2=(4*XL-XL**2+7D0*XL**3)/(10D0*XLP1**3)
        G3=(32D0-179*XL+532D0*XL**2-356*XL**3+221D0*XL**4)
     &       /(500D0*XLP1**3)
        ISTEP=1
      ELSE
        G2=(7D0-XL+4D0*XL**2)/(10D0*XLP1**3)
        G3=(221D0-356D0*XL+532D0*XL**2-179D0*XL**3+32D0*XL**4)
     &       /(500D0*XL*XLP1**3)
C These expressions essentially invert lambda
        ISIGN =-ISIGN
        ISTEP=-1
      ENDIF
      T=2D0/S
      IF(I.NE.J) THEN
        ZT=T*T
        RYG=-(G2+G3*(S*S/4D0-1D0))*DLOG(1D0-4D0/(S*S))-G3
        MMIN=2
        MMAX=MAXS
      ELSE
        ZT=T
        RYG=(G2+G3*(S*S/4D0-1D0))*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
        ZL=XLP1**(-M)
        COEF=0D0
C Some computers will underflow in the next loop.
        DO 50 IQ=0,M
          COEF=COEF+RYGCOF(NM)*ZL
          NM=NM+ISTEP
          ZL=ZL*XL
   50   CONTINUE
        RYG=RYG+ZT*(COEF-((G2SUM-G3SUM/DBLE(M+2))/DBLE(M)))
        ZT=ZT*T
  100 CONTINUE
      RYG=DBLE(ISIGN)*RYG
      IF(I.EQ.J) RETURN
      RYG=-RYG*4D0/(XLP1**2)
      IF(I.EQ.2)RYG=RYG*XL**2
      RETURN
      END