      REAL FUNCTION RF(X,Y,Z,ERRTOL,IERR)
C
C          THIS FUNCTION SUBROUTINE COMPUTES THE INCOMPLETE ELLIPTIC
C          INTEGRAL OF THE FIRST KIND,
C          RF(X,Y,Z) = INTEGRAL FROM ZERO TO INFINITY OF
C
C                                -1/2     -1/2     -1/2
C                      (1/2)(T+X)    (T+Y)    (T+Z)    DT,
C
C          WHERE X, Y, AND Z ARE NONNEGATIVE AND AT MOST ONE OF THEM
C          IS ZERO.  IF ONE OF THEM IS ZERO, THE INTEGRAL IS COMPLETE.
C          THE DUPLICATION THEOREM IS ITERATED UNTIL THE VARIABLES ARE
C          NEARLY EQUAL, AND THE FUNCTION IS THEN EXPANDED IN TAYLOR
C          SERIES TO FIFTH ORDER.  REFERENCE: B. C. CARLSON, COMPUTING
C          ELLIPTIC INTEGRALS BY DUPLICATION, NUMER. MATH. 33 (1979),
C          1-16.  CODED BY B. C. CARLSON AND ELAINE M. NOTIS, AMES
C          LABORATORY-DOE, IOWA STATE UNIVERSITY, AMES, IOWA 50011.
C          MARCH 1, 1980.
C
      INTEGER IERR
      REAL C1,C2,C3,E2,E3,EPSLON,ERRTOL,LAMDA
      REAL LOLIM,MU,S,UPLIM,X,XN,XNDEV,XNROOT
      REAL Y,YN,YNDEV,YNROOT,Z,ZN,ZNDEV,ZNROOT
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, Y, AND Z ARE THE VARIABLES IN THE INTEGRAL RF(X,Y,Z).
C
C          ERRTOL IS SET TO THE DESIRED ERROR TOLERANCE.
C          RELATIVE ERROR DUE TO TRUNCATION IS LESS THAN
C          ERRTOL ** 6 / (4 * (1 - ERRTOL)).
C
C          SAMPLE CHOICES   ERRTOL   RELATIVE TRUNCATION
C                                    ERROR LESS THAN
C                           1.E-3    3.E-19
C                           3.E-3    2.E-16
C                           1.E-2    3.E-13
C                           3.E-2    2.E-10
C                           1.E-1    3.E-7
C
C          OUTPUT ...
C
C          X, Y, Z, 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 (AMIN1(X,Y,Z) .LT. 0.0) GO TO 100
      IF (AMIN1(X+Y,X+Z,Y+Z) .LT. LOLIM) GO TO 100
      IF (AMAX1(X,Y,Z) .LE. UPLIM) GO TO 110
  100 RF = 0.0
      IERR = 1
      RETURN
C
  110 IERR = 0
      XN = X
      YN = Y
      ZN = Z
C
  115 MU = (XN + YN + ZN) / 3.0
      XNDEV = 2.0 - (MU + XN) / MU
      YNDEV = 2.0 - (MU + YN) / MU
      ZNDEV = 2.0 - (MU + ZN) / MU
      EPSLON = AMAX1(ABS(XNDEV),ABS(YNDEV),ABS(ZNDEV))
      IF (EPSLON .LT. ERRTOL) GO TO 120
      XNROOT = SQRT(XN)
      YNROOT = SQRT(YN)
      ZNROOT = SQRT(ZN)
      LAMDA = XNROOT * (YNROOT + ZNROOT) + YNROOT * ZNROOT
      XN = (XN + LAMDA) * 0.25
      YN = (YN + LAMDA) * 0.25
      ZN = (ZN + LAMDA) * 0.25
      GO TO 115
C
  120 C1 = 1.0 / 24.0
      C2 = 3.0 / 44.0
      C3 = 1.0 / 14.0
      E2 = XNDEV * YNDEV - ZNDEV * ZNDEV
      E3 = XNDEV * YNDEV * ZNDEV
      S = 1.0 + (C1 * E2 - 0.1 - C2 * E3) * E2 + C3 * E3
      RF = S / SQRT(MU)
      RETURN
      END
