PROGRAM RZM 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 and s as set in the DATA statements.
C
C The program is set up to write the output to FORTRAN-77 unit *, which
C is standard output, typically a console. For output to a file or printer
C the unit can be changed, or the output re-directed by the operating system.
C
      DOUBLE PRECISION DLAMDA(5),S(5),RZM,RZM11,RZM12,RZM21,RZM22,SUM
      INTEGER I,J
      DATA DLAMDA/1D0,.5D0,.25D0,.125D0,2D0/
      DATA S/2D0,2.01D0,2.1D0,2.5D0,3D0/
C
      WRITE(*,20)
20    FORMAT(1X,'TABLE OF VALUES FOR RESISTANCE FUNCTIONS ZM')
      DO 100 J=1,5
        WRITE(*,30)DLAMDA(J)
30      FORMAT(/1X,'LAMBDA =',F8.3)
        WRITE(*,40)
40      FORMAT(/5X,'S',8X,
     &        'RZM(1,1)',6X,'RZM(1,2)',6X,'RZM(2,1)',6X,'RZM(2,2)')
        DO 100 I=1,5
          RZM11=RZM(1,1,S(I),DLAMDA(J))
          RZM12=RZM(1,2,S(I),DLAMDA(J))
          RZM21=RZM(2,1,S(I),DLAMDA(J))
          RZM22=RZM(2,2,S(I),DLAMDA(J))
          SUM = RZM11 + 0.125D0*(DLAMDA(J)+1D0 )**3*RZM12
          WRITE(*,50)S(I),RZM11,RZM12,RZM21,RZM22,SUM
50        FORMAT(F8.3,5F14.5)
100   CONTINUE
      STOP
      END
C
      DOUBLE PRECISION FUNCTION RZM(IALPHA,IBETA,S,DLAMDA)
C ----------------------------- Arguments of function -----------------
      INTEGER IALPHA, IBETA
      DOUBLE PRECISION S, DLAMDA
C ---------------------------------------------------------------------
C                                        M
C This routine calculates the functions Z          (s,lambda)
C                                        ALPHA BETA
C  defined in Jeffrey, 'The calculation of the low-Reynolds-number
C  resistance functions for two unequal spheres', Phys. Fluids (1991).
C  Coefficients calculated by one of the programs PRPSTx are summed using
C  equations (80a,80b). The coefficients are expected in RZM200.DAT in
C  the format on the  distribution disk.
C ------------------------ Local variables ----------------------------
      DOUBLE PRECISION RZMCOF(0:45450)
      DOUBLE PRECISION T, ZT, ZL, COEF, XL, XLP1,EPS, EPSLOG
      INTEGER I, J, ISTEP, IQ, MMIN, MMAX, M, NM,G3,G3SUM
      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
      T=2D0/S
      EPS=1D0-4D0/(S*S)
      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
      G3SUM=4D0*G3
      IF(I.EQ.J) THEN
        ZT=T*T
C  Notice that because the first coef is a signature we supply the 1.
        IF(EPS.GT.1E-10) THEN
          EPSLOG=EPS*DLOG(EPS)
        ELSE
          EPSLOG=0D0
        ENDIF
        RZM=1D0 - G3*EPSLOG*S*S/4D0 - G3
        MMIN=2
        MMAX=MAXS
      ELSE
        ZT=T
        IF(EPS.GT.1E-10) THEN
          RZM=(G3*EPS*S*S/4D0)*DLOG((S+2D0)/(S-2D0)) - G3*S
        ELSE
          RZM= - G3*S
        ENDIF
        MMIN=1
        MMAX=MAXS-1
      ENDIF
      T=T*T
      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
        DO 50 IQ=0,M
          COEF=COEF+RZMCOF(NM)*ZL
          NM=NM+ISTEP
          ZL=ZL*XL
  50    CONTINUE
        COEF=COEF*XLP1**(-M)
        RZM=RZM + ZT*(COEF+G3SUM/(DBLE(M)*DBLE(M+2)))
        ZT=ZT*T
  100 CONTINUE
      IF(I.EQ.J)RETURN
      RZM=-RZM*8D0/XLP1**3
      IF(I.EQ.2)RZM=RZM*XL**3
      RETURN
      END
C
      BLOCK DATA RZMINI
      DOUBLE PRECISION RZMCOF
      COMMON /ZMD200/RZMCOF(0:45450)
      DATA RZMCOF(1)/0D0/
      END