PROGRAM RPX DEMO
C This program demonstrates the use of the coefficients calculated
C  by the programs axitrx or axistx. 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),RPX,RPX11,RPX12,RPX21,RPX22
      INTEGER J
      DATA DLAMDA/.5D0,1D0,2D0,3D0,4D0,5D0,10D0,20D0,100D0/
      WRITE(6,10)
10    FORMAT(1X,'  LAMBDA    RPX11     RPX12     RPX21     RPX22'/)
      DO 50 J=1,9
        RPX11=RPX(1,1,DLAMDA(J))
        RPX12=RPX(1,2,DLAMDA(J))
        RPX21=RPX(2,1,DLAMDA(J))
        RPX22=RPX(2,2,DLAMDA(J))
   50 WRITE(6,101)DLAMDA(J),RPX11,RPX12,RPX21,RPX22
  101 FORMAT(1X,F8.3, 4(F10.5) )
      STOP
      END
      DOUBLE PRECISION FUNCTION RPX(IALPHA,IBETA,DLAMDA)
C ------------------ Arguments of function ---------------------------
      DOUBLE PRECISION DLAMDA
      INTEGER IALPHA,IBETA
C --------------------------------------------------------------------
C                                           X
C This subprogram calculates the functions P            defined in
C                                           ALPHA BETA
C  Jeffrey,   Phys. Fluids (1992).
C ---------------------- Local variables --------------------------------
      DOUBLE PRECISION RXPCOF(0:45450)
      DOUBLE PRECISION G1,G2,G3,G1SUM,G2SUM,G3SUM
      DOUBLE PRECISION ZL,COEF, XL,XLP1,XL13, XL12
      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(1,FILE='rxp300.dat',STATUS='OLD')
        READ(1,10) ID,MAXS
 10     FORMAT(A4,I5)
        IF(ID.NE.' RXP') THEN
           WRITE(*,20)
 20        FORMAT(1X,'WRONG COEFFICIENTS FOR RPX')
           STOP
        ENDIF
        READ(1,30) RYpCOF
 30     FORMAT(D22.16)
        CLOSE(1)
        LOADED=.TRUE.
      ENDIF
      IF(DLAMDA.GT.1D0) THEN
        I=3-IALPHA
        J=3-IBETA
        XL=1D0/DLAMDA
C RXP and RPX change sign when we invert lambda.
        ISIGN = -1
      ELSE
        I=IALPHA
        J=IBETA
        XL=DLAMDA
        ISIGN = 1
      ENDIF
      XLP1=XL+1D0
      XL12=XLP1**2
      XL13=XLP1**3
C These are the g_i for X^G_{11,12} (i=1..3 in paper)
      IF(I.EQ.1) THEN
          G1=3D0*XL**2/XL13
          G2=3d0*(XL-4D0*XL**2)/(10D0*XL12)
          G3=(5D0-97D0*XL+64D0*XL**2-44D0*XL**3+XL**4)/(140D0*XL12)
         ISTEP=1
      ELSE
C These are the g_i for X^G_{21,22}
          G1=3D0*XL/XL13
          G2=3d0*(XL-4D0)/(10D0*XL12)
          G3=(5D0*XL**4-97D0*XL**3+64D0*XL**2-44D0*XL+1D0)
     &     /(140D0*XL12*XL**2)
C These expressions essentially invert lambda.
         ISIGN=-ISIGN
         ISTEP=-1
      ENDIF
      IF(I.NE.J) THEN
        RPX=0D0-G1/4D0 -G3
        MMIN=2
        MMAX=MAXS
      ELSE
        RPX=G1/4D0+G2*DLOG(4D0)-2D0*G3
        MMIN=1
        MMAX=MAXS-1
      ENDIF
      G1SUM=G1
      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+RXPCOF(NM)*ZL
          NM=NM+ISTEP
          ZL=ZL*XL
   50   CONTINUE
        RPX=RPX+COEF-G1SUM-(G2SUM-G3SUM/DBLE(M+2))/DBLE(M)
  100 CONTINUE
      RPX = DBLE(ISIGN)*RPX
      IF(I.EQ.J)RETURN
      RPX=-RPX*4D0/XLP1**2
      IF(I.EQ.2)RPX=RPX*XL**2
      RETURN
      END
C
      BLOCK DATA RPXINI
      DOUBLE PRECISION RXPCOF(45450)
      COMMON /XPD250/RXPCOF
      DATA RXPCOF(1)/0D0/
      END