/*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 <math.h>
#include "fcrt.h"
#include "srcval.h"
#include <float.h>
#include <stdlib.h>
/* 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 */