/*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 #include "fcrt.h" #include "drfval.h" #include #include void /*FUNCTION*/ drfval( double x, double y, double z, double *rf, long *ierr) { long int _l0; static double 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. *>> 1998-10-29 DRFVAL Krogh Moved external statement up for mangle. *>> 1995-11-17 DRFVAL Krogh Converted SFTRAN to Fortran 77. *>> 1994-11-02 DRFVAL Krogh Changes to use M77CON *>> 1994-08-15 DRFVAL WV Snyder JPL use 2-arg min and max for C conv. *>> 1993-02-23 DRFVAL WV Snyder JPL Split DRFVLX into a separate subr. *>> 1991-10-31 DRFVAL WV Snyder JPL Incorporate changes from Carlson *>> 1990-11-27 DRFVAL WV Snyder JPL Convert Carlson's code for MATH77 * *--D replaces "?": ?ERM1, ?ERV1, ?RFVAL, ?RFVLX * * CHECK ARGUMENT RANGES AND CALL DRFVLX TO DO THE COMPUTATION * OF CARLSON'S ELLIPTIC INTEGRAL RF(X,Y,Z). * * INPUT ... * * X, Y, AND Z ARE THE VARIABLES IN THE INTEGRAL RF(X,Y,Z). * * OUTPUT ... * * RF 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, X+Z, OR Y+Z IS TOO SMALL. * IERR = 3 X, Y, OR Z IS TOO LARGE. * */ /* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- * * MACHINE DEPENDENT PARAMETERS ... * * LOLIM AND UPLIM DETERMINE THE RANGE OF VALID ARGUMENTS. * LOLIM IS NOT LESS THAN THE MACHINE MINIMUM MULTIPLIED BY 5. * UPLIM IS NOT GREATER THAN THE MACHINE MAXIMUM DIVIDED BY 5. * */ /* ---------------------------------------------------------------------- * WARNING. CHANGES IN THE PROGRAM MAY IMPROVE SPEED AT THE * EXPENSE OF ROBUSTNESS. * ---------------------------------------------------------------------- * */ *rf = 0.0e0; if (fmin( x, fmin( y, z ) ) < 0.0e0) { *ierr = 1; derm1( "DRFVAL", 1, 0, "One of X, Y or Z is negative", "X" , x, ',' ); derv1( "Y", y, ',' ); derv1( "Z", z, '.' ); return; } if (lolim < 0.0e0) { lolim = 5.0e0*DBL_MIN; uplim = DBL_MAX/5.0e0; } if (fmax( x, fmax( y, z ) ) > uplim) { *ierr = 3; derm1( "DRFVAL", 3, 0, "One of X, Y or Z > UPLIM", "X", x, ',' ); derv1( "Y", y, ',' ); derv1( "Z", z, ',' ); derv1( "UPLIM", uplim, '.' ); return; } if (fmin( x + y, fmin( x + z, y + z ) ) < lolim) { *ierr = 2; derm1( "DRFVAL", 2, 0, "One of X+Y, X+Z or Y+Z < LOLIM", "X" , x, ',' ); derv1( "Y", y, ',' ); derv1( "Z", z, ',' ); derv1( "LOLIM", lolim, '.' ); return; } *ierr = 0; drfvlx( x, y, z, rf ); return; } /* end of function */ /* PARAMETER translations */ #define C1 (10.0e0/24.0e0) #define C2 (30.0e0/44.0e0) #define C3 (10.0e0/14.0e0) /* end of PARAMETER translations */ void /*FUNCTION*/ drfvlx( double x, double y, double z, double *rf) { long int _l0; double e2, e3, epslon, lamda, mu, s, xn, xndev, xnroot, yn, yndev, ynroot, zn, zndev, znroot; static double errtol = -1.0e0; /*>> 1993-02-23 DRFVLX WV Snyder JPL Split DRFVLX into a separate subr. * ------------------------------------------------------------------ * DRFVLX SHOULD ONLY BE CALLED FROM OTHER LIBRARY ROUTINES, WHERE WE * KNOW X, Y, Z ARE IN RANGE. USERS SHOULD NOT CALL DRFVLX DIRECTLY. * ------------------------------------------------------------------ * * COMPUTE THE INCOMPLETE ELLIPTIC INTEGRAL OF THE FIRST KIND * * RF(X,Y,Z) = INTEGRAL FROM ZERO TO INFINITY OF * * -1/2 -1/2 -1/2 * (1/2)(T+X) (T+Y) (T+Z) DT, * * WHERE X, Y, AND Z ARE NONNEGATIVE AND AT MOST ONE OF THEM * IS ZERO. IF ONE OF THEM 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. * REFERENCES: 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: RF(0,1,2) = 1.31102 87771 46059 90523 24198 * CHECK BY ADDITION THEOREM: RF(X,X+Z,X+W) + RF(Y,Y+Z,Y+W) * = RF(0,Z,W), WHERE X,Y,Z,W ARE POSITIVE AND X * Y = Z * W. * */ /* INPUT ... * * X, Y, AND Z ARE THE VARIABLES IN THE INTEGRAL RF(X,Y,Z). * * OUTPUT ... * * RF IS THE VALUE OF THE INCOMPLETE ELLIPTIC INTEGRAL. * * ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- * MACHINE DEPENDENT PARAMETER ... * ERRTOL IS SET TO THE DESIRED ERROR TOLERANCE. * RELATIVE ERROR DUE TO TRUNCATION IS LESS THAN * ERRTOL ** 6 / (4 * (1 - ERRTOL)). * */ /* ---------------------------------------------------------------------- */ if (errtol < 0.0e0) errtol = pow(3.6e0*DBL_EPSILON,1.0e0/6.0e0); xn = x; yn = y; zn = z; L_20: ; mu = (xn + yn + zn)/3.0e0; xndev = 2.0e0 - (mu + xn)/mu; yndev = 2.0e0 - (mu + yn)/mu; zndev = 2.0e0 - (mu + zn)/mu; epslon = fmax( fabs( xndev ), fmax( fabs( yndev ), fabs( zndev ) ) ); if (epslon < errtol) { e3 = xndev*yndev; e2 = e3 - zndev*zndev; e3 *= zndev; s = 10.0e0 + (C1*e2 - 1.0e0 - C2*e3)*e2 + C3*e3; *rf = s/(10.0e0*sqrt( mu )); return; } xnroot = sqrt( xn ); ynroot = sqrt( yn ); znroot = sqrt( zn ); lamda = xnroot*(ynroot + znroot) + ynroot*znroot; xn = (xn + lamda)*0.25e0; yn = (yn + lamda)*0.25e0; zn = (zn + lamda)*0.25e0; goto L_20; } /* end of function */