# rj.f

          REAL FUNCTION RJ (X, Y, Z, P, IER)
    C***BEGIN PROLOGUE  RJ
    C***PURPOSE  Compute the incomplete or complete (X or Y or Z is zero)
    C            elliptic integral of the 3rd kind.  For X, Y, and Z non-
    C            negative, at most one of them zero, and P positive,
    C             RJ(X,Y,Z,P) = Integral from zero to infinity of
    C                                  -1/2     -1/2     -1/2     -1
    C                        (3/2)(t+X)    (t+Y)    (t+Z)    (t+P)  dt.
    C***LIBRARY   SLATEC
    C***CATEGORY  C14
    C***TYPE      SINGLE PRECISION (RJ-S, DRJ-D)
    C***KEYWORDS  COMPLETE ELLIPTIC INTEGRAL, DUPLICATION THEOREM,
    C             INCOMPLETE ELLIPTIC INTEGRAL, INTEGRAL OF THE THIRD KIND,
    C             TAYLOR SERIES
    C***AUTHOR  Carlson, B. C.
    C             Ames Laboratory-DOE
    C             Iowa State University
    C             Ames, IA  50011
    C           Notis, E. M.
    C             Ames Laboratory-DOE
    C             Iowa State University
    C             Ames, IA  50011
    C           Pexton, R. L.
    C             Lawrence Livermore National Laboratory
    C             Livermore, CA  94550
    C***DESCRIPTION
    C
    C   1.     RJ
    C          Standard FORTRAN function routine
    C          Single precision version
    C          The routine calculates an approximation result to
    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.
    C
    C
    C   2.     Calling Sequence
    C          RJ( X, Y, Z, P, IER )
    C
    C          Parameters On Entry
    C          Values assigned by the calling routine
    C
    C          X      - Single precision, nonnegative variable
    C
    C          Y      - Single precision, nonnegative variable
    C
    C          Z      - Single precision, nonnegative variable
    C
    C          P      - Single precision, positive variable
    C
    C
    C          On  Return     (values assigned by the RJ routine)
    C
    C          RJ     - Single precision approximation to the integral
    C
    C          IER    - Integer
    C
    C                   IER = 0 Normal and reliable termination of the
    C                           routine.  It is assumed that the requested
    C                           accuracy has been achieved.
    C
    C                   IER >  0 Abnormal termination of the routine
    C
    C
    C          X, Y, Z, P are unaltered.
    C
    C
    C   3.    Error Messages
    C
    C         Value of IER assigned by the RJ routine
    C
    C                  Value Assigned        Error Message Printed
    C                  IER = 1               MIN(X,Y,Z) .LT. 0.0E0
    C                      = 2               MIN(X+Y,X+Z,Y+Z,P) .LT. LOLIM
    C                      = 3               MAX(X,Y,Z,P) .GT. UPLIM
    C
    C
    C
    C   4.     Control Parameters
    C
    C                  Values of LOLIM, UPLIM, and ERRTOL are set by the
    C                  routine.
    C
    C
    C          LOLIM and UPLIM determine the valid range of X Y, Z, and P
    C
    C          LOLIM is not less than the cube root of the value
    C          of LOLIM used in the routine for RC.
    C
    C          UPLIM is not greater than 0.3 times the cube root of
    C          the value of UPLIM used in the routine for RC.
    C
    C
    C                     Acceptable Values For:   LOLIM      UPLIM
    C                     IBM 360/370 SERIES   :   2.0E-26     3.0E+24
    C                     CDC 6000/7000 SERIES :   5.0E-98     3.0E+106
    C                     UNIVAC 1100 SERIES   :   5.0E-13     6.0E+11
    C                     CRAY                 :   1.32E-822   1.4E+821
    C                     VAX 11 SERIES        :   2.5E-13     9.0E+11
    C
    C
    C
    C          ERRTOL determines the accuracy of the answer
    C
    C                 The value assigned by the routine will result
    C                 in solution precision within 1-2 decimals of
    C                 "machine precision".
    C
    C
    C
    C
    C          Relative error due to truncation of the series for RJ
    C          is less than 3 * ERRTOL ** 6 / (1 - ERRTOL) ** 3/2.
    C
    C
    C
    C              The accuracy of the computed approximation to the inte-
    C              gral can be controlled by choosing the value of ERRTOL.
    C              Truncation of a Taylor series after terms of fifth order
    C              Introduces an error less than the amount shown in the
    C              second column of the following table for each value of
    C              ERRTOL in the first column.  In addition to the trunca-
    C              tion error there will be round-off error, but in prac-
    C              tice the total error from both sources is usually less
    C              than the amount given in the table.
    C
    C
    C
    C          Sample choices:  ERRTOL   Relative Truncation
    C                                    error less than
    C                           1.0E-3    4.0E-18
    C                           3.0E-3    3.0E-15
    C                           1.0E-2    4.0E-12
    C                           3.0E-2    3.0E-9
    C                           1.0E-1    4.0E-6
    C
    C                    Decreasing ERRTOL by a factor of 10 yields six more
    C                    decimal digits of accuracy at the expense of one or
    C                    two more iterations of the duplication theorem.
    C
    C *Long Description:
    C
    C   RJ Special Comments
    C
    C
    C          Check by addition theorem: RJ(X,X+Z,X+W,X+P)
    C          + RJ(Y,Y+Z,Y+W,Y+P) + (A-B) * RJ(A,B,B,A) + 3 / SQRT(A)
    C          = RJ(0,Z,W,P), where X,Y,Z,W,P are positive and X * Y
    C          = Z * W,  A = P * P * (X+Y+Z+W),  B = P * (P+X) * (P+Y),
    C          and B - A = P * (P-Z) * (P-W).  The sum of the third and
    C          fourth terms on the left side is 3 * RC(A,B).
    C
    C
    C          On Input:
    C
    C          X, Y, Z, and P are the variables in the integral RJ(X,Y,Z,P).
    C
    C
    C          On Output:
    C
    C
    C          X, Y, Z, and P are unaltered.
    C
    C          ********************************************************
    C
    C          Warning: Changes in the program may improve speed at the
    C                   expense of robustness.
    C
    C ------------------------------------------------------------
    C
    C
    C   Special Functions via RJ and RF
    C
    C
    C                  Legendre form of ELLIPTIC INTEGRAL of 3rd kind
    C                  ----------------------------------------------
    C
    C
    C                               PHI         2         -1
    C                  P(PHI,K,N) = INT (1+N SIN (THETA) )   *
    C                                0
    C
    C                                      2    2         -1/2
    C                                 *(1-K  SIN (THETA) )     D THETA
    C
    C
    C                                         2          2   2
    C                       = SIN (PHI) RF(COS (PHI), 1-K SIN (PHI),1)
    C
    C                                  3            2         2   2
    C                        -(N/3) SIN (PHI) RJ(COS (PHI),1-K SIN (PHI),
    C
    C                                 2
    C                        1,1+N SIN (PHI))
    C
    C
    C
    C                  Bulirsch form of ELLIPTIC INTEGRAL of 3rd kind
    C                  ----------------------------------------------
    C
    C
    C                                           2 2    2
    C                  EL3(X,KC,P) = X RF(1,1+KC X ,1+X ) +
    C
    C                                            3          2 2    2     2
    C                               +(1/3)(1-P) X  RJ(1,1+KC X ,1+X ,1+PX )
    C
    C
    C                                           2
    C                  CEL(KC,P,A,B) = A RF(0,KC ,1) +
    C
    C                                                     2
    C                                 +(1/3)(B-PA) RJ(0,KC ,1,P)
    C
    C
    C
    C
    C                  Heuman's LAMBDA function
    C                  ------------------------
    C
    C
    C                                 2                     2      2    1/2
    C                  L(A,B,P) = (COS(A)SIN(B)COS(B)/(1-COS (A)SIN (B))   )
    C
    C                                           2         2       2
    C                            *(SIN(P) RF(COS (P),1-SIN (A) SIN (P),1)
    C
    C                                 2       3            2       2
    C                            +(SIN (A) SIN (P)/(3(1-COS (A) SIN (B))))
    C
    C                                   2         2       2
    C                            *RJ(COS (P),1-SIN (A) SIN (P),1,1-
    C
    C                                2       2          2       2
    C                            -SIN (A) SIN (P)/(1-COS (A) SIN (B))))
    C
    C
    C
    C
    C                  (PI/2) LAMBDA0(A,B) =L(A,B,PI/2) =
    C
    C
    C                    2                         2       2    -1/2
    C               = COS (A)  SIN(B) COS(B) (1-COS (A) SIN (B))
    C
    C                           2                  2       2
    C                  *RF(0,COS (A),1) + (1/3) SIN (A) COS (A)
    C
    C                                       2       2    -3/2
    C                  *SIN(B) COS(B) (1-COS (A) SIN (B))
    C
    C                           2         2       2          2       2
    C                  *RJ(0,COS (A),1,COS (A) COS (B)/(1-COS (A) SIN (B)))
    C
    C
    C
    C                  Jacobi ZETA function
    C                  --------------------
    C
    C
    C                             2                     2   2    1/2
    C                  Z(B,K) = (K/3) SIN(B) COS(B) (1-K SIN (B))
    C
    C
    C                                      2      2   2                2
    C                             *RJ(0,1-K ,1,1-K SIN (B)) / RF (0,1-K ,1)
    C
    C
    C    -------------------------------------------------------------------
    C
    C***REFERENCES  B. C. Carlson and E. M. Notis, Algorithms for incomplete
    C                 elliptic integrals, ACM Transactions on Mathematical
    C                 Software 7, 3 (September 1981), pp. 398-403.
    C               B. C. Carlson, Computing elliptic integrals by
    C                 duplication, Numerische Mathematik 33, (1979),
    C                 pp. 1-16.
    C               B. C. Carlson, Elliptic integrals of the first kind,
    C                 SIAM Journal of Mathematical Analysis 8, (1977),
    C                 pp. 231-242.
    C***ROUTINES CALLED  R1MACH, RC, XERMSG
    C***REVISION HISTORY  (YYMMDD)
    C   790801  DATE WRITTEN
    C   890531  Changed all specific intrinsics to generic.  (WRB)
    C   891009  Removed unreferenced statement labels.  (WRB)
    C   891009  REVISION DATE from Version 3.2
    C   891214  Prologue converted to Version 4.0 format.  (BAB)
    C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
    C   900326  Removed duplicate information from DESCRIPTION section.
    C           (WRB)
    C   900510  Changed calls to XERMSG to standard form, and some
    C           editorial changes.  (RWC)).
    C   920501  Reformatted the REFERENCES section.  (WRB)
    C***END PROLOGUE  RJ