PROGRAM RMZ DEM
C This program demonstrates the use of the coefficients calculated
C  by the programs prpstx. 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),RMZ,RMZ11,RMZ12,RMZ21,RMZ22
      INTEGER J
      DATA DLAMDA/.5D0,1D0,2D0,3D0,4D0,5D0,10D0,20D0,100D0/
      WRITE(6,5)
 5    FORMAT(1X,
     &'RESISTANCE FUNCTIONS RMZ from Jeffrey,'
     &1X,'Phys. Fluids, (1992), vol A4, p16.'
     &/1X,'Tabulation of RMZ appearing in eqs. (79) & (81).') 
      WRITE(6,10)
10    FORMAT(/1X,'  LAMBDA    RMZ11     RMZ12     RMZ21     RMZ22'/)
      DO 50 J=1,9
        RMZ11=RMZ(1,1,DLAMDA(J))
        RMZ12=RMZ(1,2,DLAMDA(J))
        RMZ21=RMZ(2,1,DLAMDA(J))
        RMZ22=RMZ(2,2,DLAMDA(J))
   50 WRITE(6,101)DLAMDA(J),RMZ11,RMZ12,RMZ21,RMZ22
  101 FORMAT(1X,F8.3, 4(F10.5) )
      STOP
      END
C
      DOUBLE PRECISION FUNCTION RMZ(IALPHA,IBETA,DLAMDA)
C ---------------------- Arguments of function -------------------------
      DOUBLE PRECISION DLAMDA
      INTEGER IALPHA,IBETA
C ----------------------------------------------------------------------
C This program calculates the function RMZ defined in Jeffrey, 'The
C  calculation of the low-Reynolds-number resistance functions for two
C  unequal spheres',Phys. Fluids (1991). The coefficients are summed using
C  equations (81a,81b). The coefficients must conform to those generated by
C  one of the programs PRPSTx. The subprogram works by converting any call
C  to one in which lambda.LT.1. It goes forwards or backwards through the
C  coefficients (the variable ISTEP) depending upon whether I=1 or 2.
C ------------------------- Local variables ----------------------------
      DOUBLE PRECISION RZMCOF(0:45450)
      DOUBLE PRECISION XL,XLP1,G3,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(1,FILE='rzm300.dat',STATUS='OLD')
        READ(1,10) ID,MAXS
 10     FORMAT(A4,I5)
        IF(ID.NE.' RZM') THEN
           WRITE(*,20)
 20        FORMAT(1X,'WRONG COEFFICIENTS FOR RZM')
           STOP
        ENDIF
        READ(1,30) RZMCOF
 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.1) THEN
        G3=-3D0*(XL**2+XL**4)/(10D0*XLP1**3)
        ISTEP=1
      ELSE
        G3=-3D0*(XL**2+1D0)/(10D0*XL*XLP1**3)
        ISTEP=-1
      ENDIF
      IF(I.EQ.J) THEN
C       Notice that because the first coef is a signature we supply the 1.
        RMZ=1D0-G3
        MMIN=2
        MMAX=MAXS
      ELSE
        RMZ=2D0*(-G3)
        MMIN=1
        MMAX=MAXS-1
      ENDIF
      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+RZMCOF(NM)*ZL
          NM=NM+ISTEP
          ZL=ZL*XL
   50   CONTINUE
        RMZ=RMZ+COEF+(G3SUM/DBLE(M+2))/DBLE(M)
  100 CONTINUE
      IF(I.EQ.J)RETURN
      RMZ=-RMZ*8D0/XLP1**3
      IF(I.EQ.2)RMZ=RMZ*XL**3
      RETURN
      END
C
      BLOCK DATA RMZINI
      DOUBLE PRECISION RZMCOF(0:45450)
      COMMON /ZMD300/RZMCOF
      DATA RZMCOF(1)/0D0/
      END