PROGRAM RYC 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 and s as set in the DATA statements.
C
      DOUBLE PRECISION DLAMDA(5),S(4),RYC,RYC11,RYC12,RYC21,RYC22
      INTEGER I,J
      DATA DLAMDA/1D0,.5D0,.25D0,.125D0,2D0/
      DATA S/2.01D0,2.1D0,2.5D0,3D0/
      WRITE(6,200)
  200 FORMAT(1X,'Table of values for Resistance Functions YC')
      DO 50 J=1,5
        WRITE(6,1000)DLAMDA(J)
 1000   FORMAT(/1X,'LAMBDA =',F8.3)
        WRITE(6,100)
  100   FORMAT(/10X,'S',10X,'RYC(1,1)',8X,'RYC(1,2)',8X,'RYC(2,1)',
     &   8X,'RYC(2,2)')
        DO 50 I=1,4
          RYC11=RYC(1,1,S(I),DLAMDA(J))
          RYC12=RYC(1,2,S(I),DLAMDA(J))
          RYC21=RYC(2,1,S(I),DLAMDA(J))
          RYC22=RYC(2,2,S(I),DLAMDA(J))
   50 WRITE(6,1100)S(I),RYC11,RYC12,RYC21,RYC22
 1100 FORMAT(F13.3,4F16.5)
      STOP
      END
      DOUBLE PRECISION FUNCTION RYC(IALPHA,IBETA,S,DLAMDA)
C ----------------------- Arguments of function -------------------------
      DOUBLE PRECISION DLAMDA,S
      INTEGER IALPHA,IBETA
C -----------------------------------------------------------------------
C                                        C
C This routine calculates the functions Y          (s,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 trnrox are summed using
C  equations (7.14)-(7.15). The subprogram uses equation (1.9a) to convert
C  any call to one in which lambda.LT.1. It goes forwards or backwards
C  through the coefficients ( ISTEP) depending upon whether I=1 or 2.
C  This technique is not explained in the paper; I have forgotten why.
C  Starting at equn 1.9 we see that Y   can be obtained from Y
C                                    22                       11
C  by inverting lambda. Now looking at equn 3.15 we substitute 1/lambda and
C  redefine q as k-q. Factoring the lambda**k, we see that lambda**q now
C  multiplies P(1,q,k-q), which means running 'backwards' through the
C  coefficients relative to the published case.
C
C The coefficients are read in from rycxxx.dat where the xxx refers 
C  to the maximum summation index. This file has varied in format, so
C  the program checks first.
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
      T=2D0/S
      EPS=1D0-4D0/(S*S)
      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
        ZT=T*T
C    Notice that because the first coef is a signature we supply the 1.
        RYC=1D0-(G2+G3*EPS*S*S/4D0)*DLOG(EPS)-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
        ZT=T
        RYC=(G2+G3*EPS*S*S/4D0)*DLOG((S+2D0)/(S-2D0))
        RYC=RYC-G3*S
        MMIN=1
        MMAX=MAXS-1
      ENDIF
      T=T*T
      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
C Some computers may underflow in the next loop.
          COEF=COEF+RYCCOF(NM)*ZL
          NM=NM+ISTEP
          ZL=ZL*XL
   50   CONTINUE
        RYC=RYC+ZT*(COEF-((G2SUM-G3SUM/DBLE(M+2))/DBLE(M)))
        ZT=ZT*T
  100 CONTINUE
      IF(I.EQ.J)RETURN
      RYC=RYC*8D0/XLP1**3
      IF(I.EQ.2)RYC=RYC*XL**3
      RETURN
      END