PROGRAM RCX DEM
C This program demonstrates the use of the coefficients calculated
C  by the programs axirox. The program will produce a table of values
C  for different values of lambda as set in the DATA statements.
C
      DOUBLE PRECISION DLAMDA(7),RCX,RCX11,RCX12,RCX22,ZETA3
      INTEGER J
      DATA DLAMDA/1D0,.5D0,.25D0,.2D0,.125D0,.1D0,2D0/
      DO 50 J=1,7
        WRITE(6,1000)DLAMDA(J)
 1000   FORMAT(/1X,'LAMBDA =',F8.3)
        RCX11=RCX(1,1,DLAMDA(J))
        RCX12=RCX(1,2,DLAMDA(J))
        RCX22=RCX(2,2,DLAMDA(J))
        ZETA3=RCX12*(1D0+DLAMDA(J))**6/(8D0*DLAMDA(J)**3)
   50 WRITE(6,101)RCX11,RCX12,RCX22,ZETA3
  101 FORMAT(1X,'  RCX(1,1) =',F10.5,'  RCX(1,2) =',F10.5,
     &    '  RCX(2,2) =',F10.5,' ZETA3 =',F10.5)
      STOP
      END
C
      DOUBLE PRECISION FUNCTION RCX(IALPHA,IBETA,DLAMDA)
C --------------------- Arguments of function ----------------------------
      DOUBLE PRECISION DLAMDA
      INTEGER IALPHA,IBETA
C ------------------------------------------------------------------------
C                                           X
C This subprogram calculates the functions C  defined in Jeffrey and Onishi
C (1984) J. Fluid Mech vol 139, 261-290. The coefficients are generated
C by one of the programs axirox. The subprogram works by converting any
C call 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---------------------------- Local Variables -----------------------------
      DOUBLE PRECISION RXCCOF(0:45450),G3,G3SUM
      DOUBLE PRECISION XL,XLP1,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='rxc300.dat',STATUS='OLD')
        READ(1,10)  ID,MAXS
 10     FORMAT(A4,I5)
        IF (ID.NE.' RXC') THEN
           WRITE(*,20)
 20        FORMAT(1X,'WRONG COEFFICIENTS FOR RXC')
           STOP
        ENDIF
        READ(1,30) RXCCOF
 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=-XL*XL/(4D0*XLP1)
        ISTEP=1
      ELSE
        G3=-1D0/(4D0*XL*XLP1)
        ISTEP=-1
      ENDIF
      IF(I.EQ.J) THEN
C    Notice that because the first coef is a signature we supply the 1.
        RCX=1D0-G3
        MMIN=2
        MMAX=MAXS
      ELSE
        RCX=-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+RXCCOF(NM)*ZL
          NM=NM+ISTEP
          ZL=ZL*XL
   50   CONTINUE
        RCX=RCX+COEF+(G3SUM/DBLE(M+2))/DBLE(M)
  100 CONTINUE
      IF(I.EQ.J)RETURN
      RCX=-RCX*8D0/XLP1**3
      IF(I.EQ.2)RCX=RCX*XL**3
      RETURN
      END
C
      BLOCK DATA RCXINI
      DOUBLE PRECISION RXCCOF(0:45450)
      COMMON /XCD200/RXCCOF
      DATA RXCCOF(1)/0D0/
      END