PROGRAM RCY DEM
C This program demonstrates the use of the coefficients calculated
C  by the programs trnrox. 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),RCY,RCY11,RCY12,RCY21,RCY22
      INTEGER J
      DATA DLAMDA/1D0,.5D0,.25D0,.2D0,.125D0,.1D0,2D0/
      WRITE(6,5)
 5    FORMAT(1X,
     &'RESISTANCE FUNCTIONS RCY from Jeffrey & Onishi, J.F.M., 1982')
      WRITE(6,7)
 7    FORMAT(1X,'vol 139, p261')
      WRITE(6,8)
 8    FORMAT(/1X,'Table of Functions RCY for 0.1<lambda<1')
      WRITE(6,10)
10    FORMAT(3X,'LAMBDA',7X,'RCY(1,1)',7X,'RCY(1,2)',7X,'RCY(2,1)',7X,
     &       'RCY(2,2)' )
      DO 50 J=1,7
        RCY11=RCY(1,1,DLAMDA(J))
        RCY12=RCY(1,2,DLAMDA(J))
        RCY21=RCY(2,1,DLAMDA(J))
        RCY22=RCY(2,2,DLAMDA(J))
   50 WRITE(6,101)DLAMDA(J),RCY11,RCY12,RCY21,RCY22
  101 FORMAT(1X,F7.3,4F15.5)
      STOP
      END
      DOUBLE PRECISION FUNCTION RCY(IALPHA,IBETA,DLAMDA)
C ---------------------- Arguments of function -------------------------
      INTEGER IALPHA,IBETA
      DOUBLE PRECISION DLAMDA
C ----------------------------------------------------------------------
C                                           Y
C This subprogram calculates the functions C  defined in
C  Jeffrey & Onishi (1984) J. Fluid Mech Vol 139, 261-290.
C  Equations (7.11, 7.12) are evaluated using coefficients read from unit 4,
C  and conforming to those generated by one of the programs TRNROx.
C  The subprogram works by converting any call to one in which lambda.LT.1
C  It goes forwards or backwards through the coefficients (the variable ISTEP)
C  depending upon whether I=1 or 2.
C ------------------------  Local variables -----------------------------
      DOUBLE PRECISION RYCCOF(0:45450),XL,XLP1,G2,G3,G2SUM,G3SUM
      DOUBLE PRECISION T,ZT,ZL,COEF,EPS
      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='ryc300.dat',STATUS='OLD')
        READ(1,10)  ID,MAXS
10      FORMAT(A4,I5)
        IF(ID.NE.' RYC') THEN
          WRITE(*,20)
20        FORMAT(1X,'WRONG COEFFICIENTS FOR RYC')
          STOP
        ENDIF
        READ(1,30) RYCCOF
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.J) THEN
        IF(I.EQ.1) THEN
          G2=2D0*XL/(5D0*XLP1)
          G3=(8D0+6D0*XL+33D0*XL**2)/(125D0*XLP1)
          ISTEP=1
        ELSE
          G2=2D0/(5D0*XLP1)
          G3=(33D0+6D0*XL+8D0*XL**2)/(125D0*XL*XLP1)
          ISTEP=-1
        ENDIF
C   Notice that because the first coef is a signature we supply the 1.
        RCY=1D0-G3
        MMIN=2
        MMAX=MAXS
      ELSE
        IF(I.EQ.1) THEN
C    These factors are (XLP1**3)/8 times g<4> and g<5> in the paper.
C    The g<5> seems to have a misprint of 4 instead of 2.
          G2=1D0*XL**2/(10D0*XLP1)
          G3=(43D0*XL-24D0*XL**2+43D0*XL**3)/(500D0*XLP1)
          ISTEP=1
        ELSE
          G2=1D0/(10D0*XL*XLP1)
          G3=(43D0-24D0*XL+43D0*XL**2)/(500D0*XL**2*XLP1)
          ISTEP=-1
        ENDIF
        RCY=2D0*(G2*DLOG(2D0)-G3)
        MMIN=1
        MMAX=MAXS-1
      ENDIF
      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
        DO 50 IQ=0,M        
          COEF=COEF+RYCCOF(NM)*ZL
          NM=NM+ISTEP
          ZL=ZL*XL
   50   CONTINUE
        RCY=RCY+COEF-(G2SUM-G3SUM/DBLE(M+2))/DBLE(M)
  100 CONTINUE
      IF(I.EQ.J)RETURN
      RCY=RCY*8D0/XLP1**3
      IF(I.EQ.2)RCY=RCY*XL**3
      RETURN
      END