PROGRAM RYB DEM
C This program demonstrates the use of the coefficients calculated
C  by either the programs trntrx or trnrox. The program will produce
C  a table of values for different values of lambda and s as set in
C  the DATA statements.
C
      DOUBLE PRECISION DLAMDA(5),S(4),RYB,RYB11,RYB12,RYB21,RYB22
      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 YB')
      DO 50 J=1,5
        WRITE(6,1000)DLAMDA(J)
 1000   FORMAT(/1X,'LAMBDA =',F8.3)
        WRITE(6,100)
  100   FORMAT(/10X,'S',10X,'RYB(1,1)',8X,'RYB(1,2)',
     &                               8X,'RYB(2,1)',8X,'RYB(2,2)')
        DO 50 I=1,4
          RYB11=RYB(1,1,S(I),DLAMDA(J))
          RYB12=RYB(1,2,S(I),DLAMDA(J))
          RYB21=RYB(2,1,S(I),DLAMDA(J))
          RYB22=RYB(2,2,S(I),DLAMDA(J))
   50 WRITE(6,1100)S(I),RYB11,RYB12,RYB21,RYB22
 1100 FORMAT(F13.3,4F16.5)
      STOP
      END
C
      DOUBLE PRECISION FUNCTION RYB(IALPHA,IBETA,S,DLAMDA)
C ---------------------- Arguments of function --------------------------
      DOUBLE PRECISION S,DLAMDA
      INTEGER IALPHA,IBETA
C -----------------------------------------------------------------------
C                                        B
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 TRNTRx are summed using
C  equations (5.9)-(5.10). 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.
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 In the program the coefficients are read in and checked
C  to make certain the signature is correct (cf TRNTR1).
C  The coefficients are expected in RYBxxx.DAT
C
C ------------------------ Local variables -----------------------------
      DOUBLE PRECISION RYBCOF(0:45450),XL,XLP1,G2,G3,G2SUM,G3SUM
      DOUBLE PRECISION T,ZT,ZL,COEF,EPS
      INTEGER I,J,ISIGN,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(3,FILE='ryb300.dat',STATUS='OLD')
        READ(3,10) ID,MAXS
 10     FORMAT(A4,I5)
        IF(ID.NE.' RYB') THEN
           WRITE(*,20)
 20        FORMAT(1X,'WRONG COEFFICIENTS FOR RYB')
           STOP
        ENDIF
        READ(3,30)RYBCOF
 30     FORMAT(D22.16)
        CLOSE(3)
        LOADED=.TRUE.
      ENDIF
      IF(DLAMDA.GT.1D0) THEN
        I=3-IALPHA
        J=3-IBETA
        XL=1D0/DLAMDA
        ISIGN=-1
      ELSE
        I=IALPHA
        J=IBETA
        XL=DLAMDA
        ISIGN=1
      ENDIF
      XLP1=XL+1D0
C
C  These Gx are the functions g<x> from Jeffrey & Onishi eqn 5.6.
      IF(I.EQ.1) THEN
        G2=-XL*(4D0+XL)/(5D0*XLP1**2)
        G3=-(32D0-33D0*XL+83D0*XL**2+43D0*XL**3)/(250D0*XLP1**2)
        ISTEP=1
      ELSE
        G2=-(4D0*XL+1D0)/(5D0*XLP1**2)
        G3=-(43D0+83D0*XL-33D0*XL**2+32D0*XL**3)/(250D0*XL*XLP1**2)
        ISTEP=-1
      ENDIF
      T=2D0/S
      EPS=1D0-4D0/(S*S)
      ISIGN=ISIGN*ISTEP
      IF(I.NE.J) THEN
        ZT=T*T
        RYB=-(G2+G3*EPS*S*S/4D0)*DLOG(EPS)-G3
        MMIN=2
        MMAX=MAXS
      ELSE
        ZT=T
        RYB=(G2+G3*EPS*S*S/4D0)*DLOG((S+2D0)/(S-2D0))
        RYB=RYB-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
          COEF=COEF+RYBCOF(NM)*ZL
          NM = NM + ISTEP
          ZL=ZL*XL
   50   CONTINUE
        RYB=RYB+ZT*(COEF-((G2SUM-G3SUM/DBLE(M+2))/DBLE(M)))
        ZT=ZT*T
 100  CONTINUE
      RYB=RYB*DBLE(ISIGN)
      IF(I.EQ.J)RETURN
      RYB=-RYB*4D0/XLP1**2
      IF(I.EQ.2)RYB=RYB*XL**2
      RETURN
      END