/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:48 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "srcval.h" #include #include /* PARAMETER translations */ #define C1 (3.0e0/10.0e0) #define C2 (1.0e0/7.0e0) #define C3 (3.0e0/8.0e0) #define C4 (9.0e0/22.0e0) /* end of PARAMETER translations */ void /*FUNCTION*/ srcval( float x, float y, float *rc, long *ierr) { long int _l0; float lamda, mu, sn, w, xn, yn; static float errtol, uplim, xlu225, y2236l; static float 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 SRCVAL Krogh Change -1.0 to -1.e0. *>> 1996-03-30 SRCVAL Krogh Added external statement. *>> 1995-11-17 SRCVAL Krogh Converted SFTRAN to Fortran 77. *>> 1994-10-19 SRCVAL Krogh Changes to use M77CON *>> 1991-10-31 SRCVAL WV Snyder JPL Incorporate changes from Carlson *>> 1990-12-20 SRCVAL WV Snyder JPL Convert from NSWC for Math 77. * * THIS SUBROUTINE COMPUTES THE ELEMENTARY INTEGRAL * * RC(X,Y) = INTEGRAL FROM ZERO TO INFINITY OF * * -1/2 -1 * (1/2)(T+X) (T+Y) DT, * * WHERE X IS NONNEGATIVE AND Y IS NONZERO. IF Y IS NEGATIVE, * THE CAUCHY PRINCIPAL VALUE IS COMPUTED BY USING A PRELIMI- * NARY TRANSFORMATION TO MAKE Y POSITIVE; SEE EQUATION (2.12) * OF THE SECOND REFERENCE BELOW. WHEN Y IS POSITIVE, THE * DUPULICATION THEOREM IS ITERATED UNTIL THE VARIABLES ARE * NEARLY EQUAL, AND THE FUNCTION IS THEN EXPANDED IN TAYLOR * SERIES TO FIFTH ORDER. LOGARITHMIC, INVERSE CIRCULAR, AND * INVERSE HYPERBOLIC FUNCTIONS ARE EXPRESSED IN TERMS OF RC * BY EQUATIONS (4.9)-(4.13) OF THE SECOND REFERENCE BELOW. * * 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 SEPT. 1, 1987. * * CHECK VALUES: RC(0,1/4) = RC(1/16,1/8) = PI, * RC(9/4,2) = LN(2), * RC(1/4,-2) = LN(2)/3. * CHECK BY ADDITION THEOREM: RC(X,X+Z) + RC(Y,Y+Z) = RC(0,Z), * WHERE X, Y, AND Z ARE POSITIVE AND X * Y = Z * Z. * * ***** Formal Arguments *********************************** * * 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. * * INPUT ... * * X AND Y ARE THE VARIABLES IN THE INTEGRAL RC(X,Y). * * OUTPUT ... * * RC IS THE VALUE OF THE INTEGRAL. * * IERR IS THE RETURN ERROR CODE. * 0 FOR NORMAL COMPLETION OF THE SUBROUTINE. * 1 X IS NEGATIVE, OR Y = 0. * 2 X+ABS(Y) IS TOO SMALL. * 3 X OR ABS(Y) IS TOO LARGE, OR X + ABS(Y) IS TOO LARGE. * 4 Y < -2.236/SQRT(LOLIM) AND 0 < X < (LOLIM*UPLIM)**2/25 * */ /*--S replaces "?": ?ERM1, ?ERV1, ?RCVAL * * ***** External References ******************************** * */ /* ***** Local Variables ************************************ * */ /* ERRTOL IS SET TO THE DESIRED ERROR TOLERANCE. * RELATIVE ERROR DUE TO TRUNCATION IS LESS THAN * 16 * ERRTOL ** 6 / (1 - 2 * ERRTOL). * * SAMPLE CHOICES ERRTOL RELATIVE TRUNCATION * ERROR LESS THAN * 1.E-3 2.E-17 * 3.E-3 2.E-14 * 1.E-2 2.E-11 * 3.E-2 2.E-8 * 1.E-1 2.E-5 * * We could put in a Newton iteration to solve for ERRTOL in * terms of R1MACH(4), but it seems good enough to put * ERRTOL = (R1MACH(4)/16.0)**(1.0/6.0) * * ---------------------------------------------------------------------- * WARNING. CHANGES IN THE PROGRAM MAY IMPROVE SPEED AT THE * EXPENSE OF ROBUSTNESS. * ---------------------------------------------------------------------- * */ if (x < 0.0e0 || y == 0.0e0) { *rc = 0.0e0; serm1( "SRCVAL", 1, 0, "X < 0 or Y = 0", "X", x, ',' ); serv1( "Y", y, '.' ); *ierr = 1; return; } if (lolim < 0.0e0) { lolim = 5.0e0*FLT_MIN; uplim = FLT_MAX/5.0e0; xlu225 = powif(lolim*uplim,2)/25.0e0; y2236l = -2.236e0/sqrtf( lolim ); errtol = powf(FLT_EPSILON/16.0e0,1.0e0/6.0e0); } yn = fabsf( y ); if (x > uplim || yn > uplim) { *rc = 0.0e0; serm1( "SRCVAL", 3, 0, "X > UPLIM or ABS(Y) > UPLIM", "X", x, ',' ); serv1( "Y", y, ',' ); serv1( "UPLIM", uplim, '.' ); *ierr = 3; return; } if (x + yn < lolim) { *rc = 0.0e0; serm1( "SRCVAL", 2, 0, "X + ABS(Y) < LOLIM", "X", x, ',' ); serv1( "Y", y, ',' ); serv1( "LOLIM", lolim, '.' ); *ierr = 2; return; } if (x + yn > uplim) { *rc = 0.0e0; serm1( "SRCVAL", 3, 0, "X + ABS(Y) > UPLIM", "X", x, ',' ); serv1( "Y", y, ',' ); serv1( "UPLIM", uplim, '.' ); *ierr = 3; return; } if ((y < y2236l && x > 0.0e0) && x < xlu225) { *rc = 0.0e0; serm1( "SRCVAL", 4, 0, "Y < -2.236/SQRT(LOLIM) AND 0 < X < (LOLIM*UPLIM)**2/25" , "X", x, ',' ); serv1( "Y", y, ',' ); serv1( "LOLIM", lolim, ',' ); serv1( "UPLIM", uplim, '.' ); *ierr = 4; return; } *ierr = 0; if (y > 0.0e0) { xn = x; w = 1.0e0; } else { /* TRANSFORM TO POSITIVE Y */ xn = x - y; w = sqrtf( x )/sqrtf( xn ); } L_20: ; mu = (xn + yn + yn)/3.0e0; sn = (yn + mu)/mu - 2.0e0; if (fabsf( sn ) < errtol) { *rc = w*(1.0e0 + sn*sn*(C1 + sn*(C2 + sn*(C3 + sn*C4))))/sqrtf( mu ); return; } lamda = 2.0*sqrtf( xn )*sqrtf( yn ) + yn; xn = (xn + lamda)*0.25e0; yn = (yn + lamda)*0.25e0; goto L_20; } /* end of function */