SLATEC Routines --- DRJ ---
*DECK DRJ
DOUBLE PRECISION FUNCTION DRJ (X, Y, Z, P, IER)
C***BEGIN PROLOGUE DRJ
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 DOUBLE 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. DRJ
C Standard FORTRAN function routine
C Double precision version
C The routine calculates an approximation result to
C DRJ(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 DRJ( X, Y, Z, P, IER )
C
C Parameters on Entry
C Values assigned by the calling routine
C
C X - Double precision, nonnegative variable
C
C Y - Double precision, nonnegative variable
C
C Z - Double precision, nonnegative variable
C
C P - Double precision, positive variable
C
C
C On Return (values assigned by the DRJ routine)
C
C DRJ - Double 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 DRJ routine
C
C Value assigned Error Message printed
C IER = 1 MIN(X,Y,Z) .LT. 0.0D0
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 DRC.
C
C UPLIM is not greater than 0.3 times the cube root of
C the value of UPLIM used in the routine for DRC.
C
C
C Acceptable values for: LOLIM UPLIM
C IBM 360/370 SERIES : 2.0D-26 3.0D+24
C CDC 6000/7000 SERIES : 5.0D-98 3.0D+106
C UNIVAC 1100 SERIES : 5.0D-103 6.0D+101
C CRAY : 1.32D-822 1.4D+821
C VAX 11 SERIES : 2.5D-13 9.0D+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 DRJ
C is less than 3 * ERRTOL ** 6 / (1 - ERRTOL) ** 3/2.
C
C
C
C The accuracy of the computed approximation to the integral
C 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 truncation
C error there will be round-off error, but in practice the
C total error from both sources is usually less than the
C amount given in the table.
C
C
C
C Sample choices: ERRTOL Relative truncation
C error less than
C 1.0D-3 4.0D-18
C 3.0D-3 3.0D-15
C 1.0D-2 4.0D-12
C 3.0D-2 3.0D-9
C 1.0D-1 4.0D-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 DRJ Special Comments
C
C
C Check by addition theorem: DRJ(X,X+Z,X+W,X+P)
C + DRJ(Y,Y+Z,Y+W,Y+P) + (A-B) * DRJ(A,B,B,A) + 3.0D0 / SQRT(A)
C = DRJ(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.0D0 * DRC(A,B).
C
C
C On Input:
C
C X, Y, Z, and P are the variables in the integral DRJ(X,Y,Z,P).
C
C
C On Output:
C
C
C X, Y, Z, 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 double precision functions via DRJ and DRF
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
C 2 2 -1/2
C *(1-K SIN (THETA) ) D THETA
C
C
C 2 2 2
C = SIN (PHI) DRF(COS (PHI), 1-K SIN (PHI),1)
C
C 3 2 2 2
C -(N/3) SIN (PHI) DRJ(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 DRF(1,1+KC X ,1+X ) +
C
C 3 2 2 2 2
C +(1/3)(1-P) X DRJ(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
C 2
C +(1/3)(B-PA) DRJ(0,KC ,1,P)
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) DRF(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 *DRJ(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 (PI/2) LAMBDA0(A,B) =L(A,B,PI/2) =
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 *DRF(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 *DRJ(0,COS (A),1,COS (A) COS (B)/(1-COS (A) SIN (B)))
C
C
C Jacobi ZETA function
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 *DRJ(0,1-K ,1,1-K SIN (B)) / DRF (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 D1MACH, DRC, 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 DRJ