/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:43 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include "drdval.h"
#include <float.h>
#include <stdlib.h>
		/* PARAMETER translations */
#define	C1	(3.0e0/14.0e0)
#define	C2	(1.0e0/6.0e0)
#define	C3	(9.0e0/22.0e0)
#define	C4	(3.0e0/26.0e0)
#define	EPWR	(-2.0e0/3.0e0)
		/* end of PARAMETER translations */
 
void /*FUNCTION*/ drdval(
double x,
double y,
double z,
double *rd,
long *ierr)
{
	long int _l0;
	double ea, eb, ec, ed, ef, epslon, lamda, mu, power4,
	 s1, s2, sigma, xn, xndev, xnroot, yn, yndev, ynroot, zn, zndev,
	 znroot;
	static double errtol, uplim;
	static double lolim = -1.0e0;
 
	/* Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
	 * ALL RIGHTS RESERVED.
	 * Based on Government Sponsored Research NAS7-03001.
	 *>>   2001-07-16 DRDVAL Krogh  Change -1.0 to -1.d0.
	 *>>   1998-10-29 DRDVAL Krogh  Moved external statement up for mangle.
	 *>>   1995-11-17 DRDVAL Krogh  Converted SFTRAN to Fortran 77.
	 *>>   1994-10-19 DRDVAL Krogh  Changes to use M77CON
	 *>>   1994-08-31 DRDVAL WV Snyder JPL use 2-arg min and max for C conv.
	 *>>   1991-10-31 DRDVAL WV Snyder JPL Incorporate changes from Carlson
	 *>>   1990-11-27 DRDVAL WV Snyder JPL Convert Carlson's code for MATH77
	 *
	 *     COMPUTE THE INCOMPLETE ELLIPTIC INTEGRAL OF THE SECOND KIND
	 *
	 *     RD(X,Y,Z) = INTEGRAL FROM ZERO TO INFINITY OF
	 *
	 *                           -1/2     -1/2     -3/2
	 *                 (3/2)(T+X)    (T+Y)    (T+Z)    DT,
	 *
	 *     WHERE X AND Y ARE NONNEGATIVE, X + Y IS POSITIVE, AND Z IS
	 *     POSITIVE.  IF X OR Y IS ZERO, THE INTEGRAL IS COMPLETE.
	 *     THE DUPLICATION THEOREM IS ITERATED UNTIL THE VARIABLES ARE
	 *     NEARLY EQUAL, AND THE FUNCTION IS THEN EXPANDED IN TAYLOR
	 *     SERIES TO FIFTH ORDER.
	 *     REFERENCE: B. C. CARLSON AND E. M. NOTIS, ALGORITHMS FOR
	 *     INCOMPLETE ELLIPTIC INTEGRALS, ACM TRANSACTIONS ON MATHEMA-
	 *     TICAL SOFTWARE, 7 (1981), 398-403; B. C. CARLSON, COMPUTING
	 *     ELLIPTIC INTEGRALS BY DUPLICATION, NUMER. MATH. 33 (1979),
	 *     1-16.
	 *     AUTHORS: B. C. CARLSON AND ELAINE M. NOTIS, AMES LABORATORY-
	 *     DOE, IOWA STATE UNIVERSITY, AMES, IA 50011, AND R. L. PEXTON,
	 *     LAWRENCE LIVERMORE NATIONAL LABORATORY, LIVERMORE, CA 94550.
	 *     AUG. 1, 1979, REVISED JAN. 15, 1987.
	 *
	 *     CHECK VALUE: RD(0,2,1) = 1.79721 03521 03388 31115 98837
	 *     CHECK: RD(X,Y,Z) + RD(Y,Z,X) + RD(Z,X,Y)
	 *     = 3 / DSQRT(X * Y * Z), WHERE X, Y, AND Z ARE POSITIVE.
	 * */
 
	/*     INPUT ...
	 *
	 *     X, Y, AND Z ARE THE VARIABLES IN THE INTEGRAL RD(X,Y,Z).
	 *
	 *     OUTPUT ...
	 *
	 *     RD IS THE VALUE OF THE INCOMPLETE ELLIPTIC INTEGRAL.
	 *
	 *     IERR IS THE RETURN ERROR CODE.
	 *          IERR = 0  FOR NORMAL COMPLETION OF THE SUBROUTINE.
	 *          IERR = 1  X, Y, OR Z IS NEGATIVE.
	 *          IERR = 2  X+Y OR Z IS TOO SMALL.
	 *          IERR = 3  X, Y, OR Z IS TOO LARGE.
	 *
	 * ----------------------------------------------------------------------
	 *--D replaces "?": ?ERM1, ?ERV1, ?RDVAL
	 * ---------------------------------------------------------------------- */
	/* ----------------------------------------------------------------------
	 *
	 *     MACHINE DEPENDENT PARAMETERS ...
	 *
	 *     ERRTOL IS SET TO THE DESIRED ERROR TOLERANCE.
	 *     RELATIVE ERROR DUE TO TRUNCATION IS LESS THAN
	 *     3 * ERRTOL ** 6 / (1 - ERRTOL) ** 3/2.
	 *
	 *     LOLIM AND UPLIM DETERMINE THE RANGE OF VALID ARGUMENTS.
	 *     LOLIM IS NOT LESS THAN 2 / (MACHINE MAXIMUM) ** (2/3).
	 *     UPLIM IS NOT GREATER THAN (0.1 * ERRTOL / MACHINE
	 *     MINIMUM) ** (2/3).
	 * */
	/* ----------------------------------------------------------------------
	 *          WARNING. CHANGES IN THE PROGRAM MAY IMPROVE SPEED AT THE
	 *          EXPENSE OF ROBUSTNESS.
	 * ----------------------------------------------------------------------
	 * */
	*rd = 0.0;
	if (fmin( x, fmin( y, z ) ) < 0.0e0)
	{
		*ierr = 1;
		derm1( "DRDVAL", 1, 0, "One of X, Y or Z is negative", "X"
		 , x, ',' );
		derv1( "Y", y, ',' );
		derv1( "Z", z, '.' );
		return;
	}
	if (lolim < 0.0)
	{
		errtol = pow(.28e0*DBL_EPSILON,1.0e0/6.0e0);
		lolim = 2.0001e0*pow(DBL_MAX,EPWR);
		uplim = pow(10.0e0*DBL_MIN/errtol,EPWR);
	}
	if (fmin( x + y, z ) < lolim)
	{
		*ierr = 2;
		derm1( "DRDVAL", 2, 0, "MIN(X+Y,Z) < LOLIM", "X", x, ',' );
		derv1( "Y", y, ',' );
		derv1( "Z", z, ',' );
		derv1( "LOLIM", lolim, '.' );
		return;
	}
	if (fmax( x, fmax( y, z ) ) > uplim)
	{
		*ierr = 3;
		derm1( "DRDVAL", 3, 0, "One of X, Y or Z > UPLIM", "X", x,
		 ',' );
		derv1( "Y", y, ',' );
		derv1( "Z", z, ',' );
		derv1( "UPLIM", uplim, '.' );
		return;
	}
 
	*ierr = 0;
	xn = x;
	yn = y;
	zn = z;
	sigma = 0.0e0;
	power4 = 1.0e0;
 
L_20:
	;
	mu = (xn + yn + 3.0e0*zn)/5.0e0;
	xndev = (mu - xn)/mu;
	yndev = (mu - yn)/mu;
	zndev = (mu - zn)/mu;
	epslon = fmax( fabs( xndev ), fmax( fabs( yndev ), fabs( zndev ) ) );
	if (epslon < errtol)
	{
		ea = xndev*yndev;
		eb = zndev*zndev;
		ec = ea - eb;
		ed = ea - 6.0e0*eb;
		ef = ed + ec + ec;
		eb = C4*zndev;
		s1 = ed*(0.25e0*C3*ed - 1.5e0*eb*ef - C1);
		s2 = zndev*(C2*ef + zndev*(eb*ea - C3*ec));
		*rd = 3.0e0*sigma + power4*(1.0e0 + s1 + s2)/(mu*sqrt( mu ));
		return;
	}
	xnroot = sqrt( xn );
	ynroot = sqrt( yn );
	znroot = sqrt( zn );
	lamda = xnroot*(ynroot + znroot) + ynroot*znroot;
	sigma += power4/(znroot*(zn + lamda));
	power4 *= 0.25e0;
	xn = (xn + lamda)*0.25e0;
	yn = (yn + lamda)*0.25e0;
	zn = (zn + lamda)*0.25e0;
	goto L_20;
} /* end of function */