      REAL FUNCTION RC(X,Y,ERRTOL,IERR)
C
C          THIS FUNCTION SUBROUTINE COMPUTES THE ELEMENTARY INTEGRAL
C          RC(X,Y) = INTEGRAL FROM ZERO TO INFINITY OF
C
C                              -1/2     -1
C                    (1/2)(T+X)    (T+Y)  DT,
C
C          WHERE X IS NONNEGATIVE AND Y IS POSITIVE.  THE DUPLICATION
C          THEOREM IS ITERATED UNTIL THE VARIABLES ARE NEARLY EQUAL,
C          AND THE FUNCTION IS THEN EXPANDED IN TAYLOR SERIES TO FIFTH
C          ORDER.  LOGARITHMIC, INVERSE CIRCULAR, AND INVERSE HYPER-
C          BOLIC FUNCTIONS CAN BE EXPRESSED IN TERMS OF RC.  REFERENCE:
C          B. C. CARLSON, COMPUTING ELLIPTIC INTEGRALS BY DUPLICATION,
C          NUMER. MATH. 33 (1979), 1-16.  CODED BY B. C. CARLSON AND
C          ELAINE M. NOTIS, AMES LABORATORY-DOE, IOWA STATE UNIVERSITY,
C          AMES, IOWA 50011.  MARCH 1, 1980.
C
      INTEGER IERR
      REAL C1,C2,ERRTOL,LAMDA,LOLIM
      REAL MU,S,SN,UPLIM,X,XN,Y,YN
C
C          LOLIM AND UPLIM DETERMINE THE RANGE OF VALID ARGUMENTS.
C          LOLIM IS NOT LESS THAN THE MACHINE MINIMUM MULTIPLIED BY 5.
C          UPLIM IS NOT GREATER THAN THE MACHINE MAXIMUM DIVIDED BY 5.
C
      LOLIM = 5.0 * SPMPAR(2)
      UPLIM = 0.2 * SPMPAR(3)
C
C          INPUT ...
C
C          X AND Y ARE THE VARIABLES IN THE INTEGRAL RC(X,Y).
C
C          ERRTOL IS SET TO THE DESIRED ERROR TOLERANCE.
C          RELATIVE ERROR DUE TO TRUNCATION IS LESS THAN
C          16 * ERRTOL ** 6 / (1 - 2 * ERRTOL).
C
C          SAMPLE CHOICES   ERRTOL   RELATIVE TRUNCATION
C                                    ERROR LESS THAN
C                           1.E-3    2.E-17
C                           3.E-3    2.E-14
C                           1.E-2    2.E-11
C                           3.E-2    2.E-8
C                           1.E-1    2.E-5
C
C          OUTPUT ...
C
C          X, Y, AND ERRTOL ARE UNALTERED.
C
C          IERR IS THE RETURN ERROR CODE
C               IERR = 0 FOR NORMAL COMPLETION OF THE SUBROUTINE,
C               IERR = 1 FOR ABNORMAL TERMINATION.
C
C          ********************************************************
C          WARNING: CHANGES IN THE PROGRAM MAY IMPROVE SPEED AT THE
C          EXPENSE OF ROBUSTNESS.
C
      IF (X .LT. 0.0 .OR. Y .LE. 0.0) GO TO 100
      IF ((X + Y) .LT. LOLIM) GO TO 100
      IF (AMAX1(X,Y) .LE. UPLIM) GO TO 110
  100 RC = 0.0
      IERR = 1
      RETURN
C
  110 IERR = 0
      XN = X
      YN = Y
C
  115 MU = (XN + YN + YN) / 3.0
      SN = (YN + MU) / MU - 2.0
      IF (ABS(SN) .LT. ERRTOL) GO TO 120
      LAMDA = 2.0 * SQRT(XN) * SQRT(YN) + YN
      XN = (XN + LAMDA) * 0.25
      YN = (YN + LAMDA) * 0.25
      GO TO 115
C
  120 C1 = 1.0 / 7.0
      C2 = 9.0 / 22.0
      S = SN * SN * (0.3 + SN * (C1 + SN * (0.375 + SN * C2)))
      RC = (1.0 + S) / SQRT(MU)
      RETURN
      END
