      REAL FUNCTION RJ(X,Y,Z,P,ERRTOL,IERR)
C
C          THIS FUNCTION SUBROUTINE COMPUTES AN INCOMPLETE ELLIPTIC
C          INTEGRAL OF THE THIRD KIND,
C          RJ(X,Y,Z,P) = INTEGRAL FROM ZERO TO INFINITY OF
C
C                                  -1/2     -1/2     -1/2     -1
C                        (3/2)(T+X)    (T+Y)    (T+Z)    (T+P)  DT,
C
C          WHERE X, Y, AND Z ARE NONNEGATIVE, AT MOST ONE OF THEM IS
C          ZERO, AND P IS POSITIVE.  IF X OR Y OR Z IS ZERO, THE
C          INTEGRAL IS COMPLETE.  THE DUPLICATION THEOREM IS ITERATED
C          UNTIL THE VARIABLES ARE NEARLY EQUAL, AND THE FUNCTION IS
C          THEN EXPANDED IN TAYLOR SERIES TO FIFTH ORDER.  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 ALFA,BETA,C1,C2,C3,C4,EA,EB,EC,E2,E3
      REAL EPSLON,ERRTOL,ETOLRC,LAMDA,LOLIM,MU,P,PN,PNDEV
      REAL POWER4,RC,SIGMA,S1,S2,S3,UPLIM,X,XN,XNDEV
      REAL XNROOT,Y,YN,YNDEV,YNROOT,Z,ZN,ZNDEV,ZNROOT
C
C          RC IS A FUNCTION COMPUTED BY AN EXTERNAL SUBROUTINE.
C
C          LOLIM AND UPLIM DETERMINE THE RANGE OF VALID ARGUMENTS.
C          LOLIM IS NOT LESS THAN THE CUBE ROOT OF THE VALUE
C          OF LOLIM USED IN THE SUBROUTINE FOR RC, AND
C          UPLIM IS NOT GREATER THAN 0.3 TIMES THE CUBE ROOT OF
C          THE VALUE OF UPLIM USED IN THE SUBROUTINE FOR RC.
C          LOLIM AND UPLIM MAY BE OBTAINED BY COMPUTING ...
C
C              LOLIM = CBRT(5.0 * SPMPAR(2))
C              UPLIM = 0.3 * CBRT(0.2 * SPMPAR(3))
C
C          IN PARTICULAR, THE FOLLOWING VALUES ARE APPROPRIATE
C          FOR THE CDC 6000-7000 SERIES COMPUTERS ...
C
      DATA LOLIM/.25E-97/, UPLIM/.40E+107/
C
C          INPUT ...
C
C          X, Y, Z, AND P ARE THE VARIABLES IN THE INTEGRAL RJ(X,Y,Z,P).
C
C          ERRTOL IS SET TO THE DESIRED ERROR TOLERANCE. THE
C          RELATIVE ERROR DUE TO TRUNCATION OF THE SERIES FOR RJ
C          IS LESS THAN 3 * ERRTOL ** 6 / (1 - ERRTOL) ** 3/2.
C          AN ERROR TOLERANCE (ETOLRC) WILL BE PASSED TO THE SUBROUTINE
C          FOR RC TO MAKE THE TRUNCATION ERROR FOR RC LESS THAN FOR RJ.
C
C          SAMPLE CHOICES   ERRTOL   RELATIVE TRUNCATION
C                                    ERROR LESS THAN
C                           1.E-3    4.E-18
C                           3.E-3    3.E-15
C                           1.E-2    4.E-12
C                           3.E-2    3.E-9
C                           1.E-1    4.E-6
C
C          OUTPUT ...
C
C          X, Y, Z, P, 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,P) .LT. LOLIM) GO TO 100
      IF (AMAX1(X,Y,Z,P) .LE. UPLIM) GO TO 110
  100 RJ = 0.0
      IERR = 1
      RETURN
C
  110 IERR = 0
      XN = X
      YN = Y
      ZN = Z
      PN = P
      SIGMA = 0.0
      POWER4 = 1.0
      ETOLRC = 0.5 * ERRTOL
C
  115 MU = (XN + YN + ZN + PN + PN) * 0.2
      XNDEV = (MU - XN) / MU
      YNDEV = (MU - YN) / MU
      ZNDEV = (MU - ZN) / MU
      PNDEV = (MU - PN) / MU
      EPSLON = AMAX1(ABS(XNDEV),ABS(YNDEV),ABS(ZNDEV),ABS(PNDEV))
      IF (EPSLON .LT. ERRTOL) GO TO 120
      XNROOT = SQRT(XN)
      YNROOT = SQRT(YN)
      ZNROOT = SQRT(ZN)
      LAMDA = XNROOT * (YNROOT + ZNROOT) + YNROOT * ZNROOT
      ALFA = PN * (XNROOT + YNROOT + ZNROOT) + XNROOT * YNROOT * ZNROOT
      ALFA = ALFA * ALFA
      BETA = PN * (PN + LAMDA) * (PN + LAMDA)
      SIGMA = SIGMA + POWER4 * RC(ALFA,BETA,ETOLRC,IERR)
      IF (IERR .NE. 0) GO TO 124
      POWER4 = POWER4 * 0.25
      XN = (XN + LAMDA) * 0.25
      YN = (YN + LAMDA) * 0.25
      ZN = (ZN + LAMDA) * 0.25
      PN = (PN + LAMDA) * 0.25
      GO TO 115
C
  120 C1 = 3.0 / 14.0
      C2 = 1.0 / 3.0
      C3 = 3.0 / 22.0
      C4 = 3.0 / 26.0
      EA = XNDEV * (YNDEV + ZNDEV) + YNDEV * ZNDEV
      EB = XNDEV * YNDEV * ZNDEV
      EC = PNDEV * PNDEV
      E2 = EA - 3.0 * EC
      E3 = EB + 2.0 * PNDEV * (EA - EC)
      S1 = 1.0 + E2 * (-C1 + 0.75 * C3 * E2 - 1.5 * C4 * E3)
      S2 = EB * (0.5 * C2 + PNDEV * (- C3 - C3 + PNDEV * C4))
      S3 = PNDEV * EA * (C2 - PNDEV * C3) - C2 * PNDEV * EC
      RJ = 3.0 * SIGMA + POWER4 * (S1 + S2 + S3) / (MU * SQRT(MU))
C
  124 RETURN
      END
