PROGRAM rax DEM
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 as set in the DATA statements.
C
C The program writes its output to FORTRAN-77  i/o 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(9),RAX,RAX11,RAX12,RAX21,RAX22,W1
      INTEGER J
      DATA DLAMDA/.5D0,1D0,2D0,3D0,4D0,5D0,10D0,20D0,100D0/
      WRITE(6,10)
10    FORMAT(1X,
     & ' LAMBDA     RAX11      RAX12      RAX21      RAX22       W1'/)
      DO 100 J=1,9
        RAX11=RAX(1,1,DLAMDA(J))
        RAX12=RAX(1,2,DLAMDA(J))
        RAX21=RAX(2,1,DLAMDA(J))
        RAX22=RAX(2,2,DLAMDA(J))
        W1  = RAX11 +(1D0+DLAMDA(J))*RAX12 + DLAMDA(J)*RAX22
        WRITE(*,50)DLAMDA(J),RAX11,RAX12,RAX21,RAX22,W1
50      FORMAT(1X,F7.3, 5(F11.5) )
100   CONTINUE
      STOP
      END
C
      DOUBLE PRECISION FUNCTION RAX(IALPHA,IBETA,DLAMDA)
C ----------------------------- Arguments of function -----------------
      INTEGER IALPHA, IBETA
      DOUBLE PRECISION DLAMDA
C ---------------------------------------------------------------------
C                                           X
C This subprogram calculates the functions A           (lambda)
C                                           ALPHA BETA
C  defined in Jeffrey & Onishi (1984) J. Fluid Mech. Vol 139, 261-290.
C  Coefficients calculated by one of the programs axitrx are summed using
C  essentially equations (3.22)-(3.23), although the form in which the
C  singularities are expanded follows Jeffrey & Corless (1988)
C  PhysicoChemical Hydrodynamics, Vol 10, 461-470.
C  The subprogram works generally the same as rxa.f where further comments
C  can be found.
C
C In the program the coefficients are read in only once and checked
C  to make certain the signature is correct (cf axitr1).
C  The coefficients are expected in rxa300.dat in the 95 format.
C     Thus MAXS = 200, RXACOF(0:20300)
C          MAXS = 250, RXACOF(0:31625)
C          MAXS = 300, RXACOF(0:45450)
C
C ------------------------ Local variables ----------------------------
      DOUBLE PRECISION RXACOF(0:45450), G1, G2, G3, G1SUM, G2SUM, G3SUM
      DOUBLE PRECISION ZL, COEF, XL, XLP1, XL13
      INTEGER I, J, ISTEP, IQ, MMIN, MMAX, M, NM, MAXS
      LOGICAL LOADED
      CHARACTER*4 ID
      DATA LOADED/.FALSE./
C ------------------------ Start of executable code -------------------
      IF(.NOT.LOADED) THEN
        OPEN(1,FILE='rxa300.dat',STATUS='OLD')
        READ(1,10)  ID,MAXS
10      FORMAT(A4,I5)
        IF(ID.NE.' RXA') THEN
          WRITE(*,20)
20        FORMAT(1X,'WRONG COEFFICIENTS FOR RXA')
          STOP
        ENDIF
        READ(1,30) RXACOF
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
      XL13=XLP1**3
      IF(I.EQ.1) THEN
C  These are the g_i for X^A_11 and X^A_12 (i=1..3 in paper)
        G1=2D0*XL**2/XL13
        G2=(XL+7D0*XL**2+XL**3)/(5D0*XL13)
        G3=(1D0+18D0*XL-29D0*XL**2+18D0*XL**3+XL**4)/(42D0*XL13)
        ISTEP=1
      ELSE
        G1=2D0*XL/XL13
        G2=(1D0+7D0*XL+XL**2)/(5D0*XL13)
        G3=(1D0+18D0*XL-29D0*XL**2+18D0*XL**3+XL**4)/(42D0*XL*XL13)
        ISTEP=-1
      ENDIF
      IF(I.EQ.J) THEN
        RAX=1D0 - G1/4D0 - G3
        MMIN=2
        MMAX=MAXS
      ELSE
        RAX=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=1D0
        COEF=0D0
C Some computers will underflow in the next statements.
        DO 50 IQ=0,M
          COEF=COEF+RXACOF(NM)*ZL
          NM=NM+ISTEP
          ZL=ZL*XL
  50    CONTINUE
        COEF=COEF*XLP1**(-M)
        RAX=RAX+COEF-G1SUM-(G2SUM-G3SUM/DBLE(M+2))/DBLE(M)
  100 CONTINUE
      IF(I.EQ.J)RETURN
      RAX=-RAX*2D0/XLP1
      IF(I.EQ.2)RAX=RAX*XL
      RETURN
      END