PROGRAM RXP DEMO
C This program demonstrates the use of the coefficients calculated
C  by the programs axitrx. The program will produce a table of values
C  for different values of lambda and s as set in the DATA statements.
      INTEGER MAXS,MAXL
      PARAMETER(MAXS=5,MAXL=6)
      DOUBLE PRECISION DLAMDA(MAXL),S(MAXS),RXP,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 XP11'/)
      WRITE(6,20)
20    FORMAT('    s',6X,25('-'),' lambda ',24('-'),/, 14X,
     & '1         2         3         4         5        1/2')
      DO 50 I=1,MAXS
        DO 40 J=1,MAXL
40      VAL(J)=RXP(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 XP12'/)
      WRITE(6,20)
      DO 150 I=1,MAXS
        DO 140 J=1,MAXL
140     VAL(J)=RXP(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 XP21'/)
      WRITE(6,20)
      DO 250 I=1,MAXS
         DO 240 J=1,MAXL
240      VAL(J)=RXP(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 XP22'/)
      WRITE(6,20)
      DO 350 I=1,MAXS
        DO 340 J=1,MAXL
340     VAL(J)=RXP(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 RXP(IALPHA,IBETA,S,DLAMDA)
C -------------------- Arguments of function ---------------------------
      DOUBLE PRECISION DLAMDA
      INTEGER IALPHA,IBETA
C ----------------------------------------------------------------------
C                                        P
C This routine calculates the functions X          (s,lambda)
C                                        ALPHA BETA
C  defined in Jeffrey, Phys. Fluids(1992).
C  Coefficients calculated by one of the programs AXITRx or AXISTx
C  are summed using equations (20a,20b).
C  The coefficients are expected in RXP200.DAT in the format on the
C  distribution disk.
C -------------------------- Local variables  -----------------------------
      DOUBLE PRECISION RXPCOF(0:45450)
      DOUBLE PRECISION G1,G2,G3,G1SUM,G2SUM,G3SUM
      DOUBLE PRECISION T,ZT,ZL,COEF,EPS,S
      DOUBLE PRECISION  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 RYC')
           STOP
        ENDIF
        READ(1,30) RXPCOF
 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 RGX 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
      T=2D0/S
      EPS=0.25D0*S*S-1D0
      IF(I.NE.J) THEN
        ZT=T*T
C Notice that because the first coef is a signature we supply the 0.
        RXP=0D0+G1/EPS-(G2+G3*EPS)*DLOG( 1D0-4D0/(S*S) )-G3
        MMIN=2
        MMAX=MAXS
      ELSE
        ZT=T
        RXP=G1*S/(2D0*EPS)+(G2+G3*EPS)*DLOG((S+2D0)/(S-2D0))-G3*S
        MMIN=1
        MMAX=MAXS-1
      ENDIF
      T=T*T
      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+2
         ENDIF
         ZL=XLP1**(-M)
         COEF=0D0
C Some computers will underflow in the next loop.
         DO 50 IQ=0,M
               NM=NM+ISTEP
               COEF=COEF+RXPCOF(NM)*ZL
               ZL=ZL*XL
   50    CONTINUE
         RXP=RXP+ZT*(COEF-G1SUM-(G2SUM-G3SUM/DBLE(M+2))/DBLE(M) )
         ZT=ZT*T
  100 CONTINUE
      RXP = DBLE(ISIGN)*RXP
      IF(I.EQ.J)RETURN
      RXP=-RXP*4D0/XLP1**2
      IF(I.EQ.2)RXP=RXP*XL**2
      RETURN
      END
C
      BLOCK DATA RXPINI
      DOUBLE PRECISION RXPCOF(45450)
      COMMON /XPD250/RXPCOF
      DATA RXPCOF(1)/0D0/
      END