PROGRAM RBY 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(7),RBY,RBY11,RBY12,RBY21,RBY22
      INTEGER J
      DATA DLAMDA/1D0,.5D0,.25D0,.2D0,.125D0,.1D0,2D0/
      WRITE(6,5)
 5    FORMAT(1X,
     &'RESISTANCE FUNCTIONS RBY from Jeffrey & Onishi, 1982')
      WRITE(6,7)
 7    FORMAT(1X,'Table of functions RBY for 0.1<lambda<1')
      WRITE(6,10)
10    FORMAT(/1X,'LAMBDA     RBY11       RBY12       RBY21      RBY22'/)
      DO 50 J=1,7
 1000   FORMAT(/1X,'LAMBDA =',F8.3)
        RBY11=RBY(1,1,DLAMDA(J))
        RBY12=RBY(1,2,DLAMDA(J))
        RBY21=RBY(2,1,DLAMDA(J))
        RBY22=RBY(2,2,DLAMDA(J))
   50 WRITE(6,101)DLAMDA(J),RBY11,RBY12,RBY21,RBY22
  101 FORMAT(1X,F6.4,4F12.5)
      STOP
      END
C
      DOUBLE PRECISION FUNCTION RBY(IALPHA,IBETA,DLAMDA)
C ------------------------- Arguments of function --------------------------
      DOUBLE PRECISION DLAMDA
      INTEGER IALPHA,IBETA
C --------------------------------------------------------------------------
C                                           Y
C This subprogram calculates the functions B  defined in
C  Jeffrey & Onishi (1984) J. Fluid Mech Vol 139, 261-290.
C  Equations 5.7 and 5.8 are evaluated using coefficients read from unit 3,
C  and conforming to those generated by one of the programs TRNTRx or 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
C  The equations (5.7)-(5.8) in the paper contain errors. The g3 coefficients
C  have been corrected as with the RYB subroutine.
C ----------------------------- Local variables ----------------------------
      DOUBLE PRECISION RYBCOF(0:45450),XL,XLP1,G2,G3,G2SUM,G3SUM
      DOUBLE PRECISION ZL,COEF
      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 equn 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
      ISIGN=ISIGN*ISTEP
      IF(I.NE.J) THEN
        RBY=-G3
        MMIN=2
        MMAX=MAXS
      ELSE
        RBY=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+RYBCOF(NM)*ZL
          NM=NM+ISTEP
          ZL=ZL*XL
   50   CONTINUE
        RBY=RBY+COEF-(G2SUM-G3SUM/DBLE(M+2))/DBLE(M)
  100 CONTINUE
      RBY=RBY*DBLE(ISIGN)
      IF(I.EQ.J)RETURN
      RBY=-RBY*4D0/XLP1**2
      IF(I.EQ.2)RBY=RBY*XL**2
      RETURN
      END