/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:59 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "srjval.h" #include #include /* PARAMETER translations */ #define C1 (3.0e0/14.0e0) #define C2 (1.0e0/3.0e0) #define C3 (3.0e0/22.0e0) #define C4 (3.0e0/26.0e0) #define FIFTH (1.0e0/5.0e0) #define THIRD (1.0e0/3.0e0) #define Z1 (3.0e0/10.0e0) #define Z2 (1.0e0/7.0e0) #define Z3 (3.0e0/8.0e0) #define Z4 (9.0e0/22.0e0) /* end of PARAMETER translations */ void /*FUNCTION*/ srjval( float x, float y, float z, float r, float *rj, long *ierr) { long int _l0; float a, alfa, b, beta, delta, e1, e2, e3, ea, eb, ec, epslon, gamma, lamda, mu, power4, rc, rcx, rf, rn, rndev, s1, s2, s3, sigma, sn, xn, xndev, xnroot, yn, yndev, ynroot, yy, zn, zndev, znroot; static float errtol, etolrc, uplim; 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 SRJVAL Krogh Change -1.0 to -1.e0. *>> 1996-03-30 SRJVAL Krogh Added external statement. *>> 1995-11-17 SRJVAL Krogh Converted SFTRAN to Fortran 77. *>> 1994-10-19 SRJVAL Krogh Changes to use M77CON *>> 1994-08-15 SRJVAL WV Snyder JPL use 2-arg min and max for C conv. *>> 1992-04-17 SRJVAL WV Snyder JPL *>> 1992-04-13 SRJVAL WV Snyder JPL Give RJ a value even if ierr.ne.0 *>> 1992-04-08 SRJVAL WV Snyder JPL Declare RCX, spell YN correctly. *>> 1990-12-20 SRJVAL WV Snyder JPL Convert from NSWC for Math 77. * * THIS SUBROUTINE COMPUTES THE INCOMPLETE ELLIPTIC INTEGRAL * OF THE THIRD KIND * * RJ(X,Y,Z,R) = INTEGRAL FROM ZERO TO INFINITY OF * * -1/2 -1/2 -1/2 -1 * (3/2)(T+X) (T+Y) (T+Z) (T+R) DT, * * WHERE X, Y, AND Z ARE NONNEGATIVE, AT MOST ONE OF THEM IS * ZERO, AND R IS NONZERO. IF X OR Y OR Z IS ZERO, THE INTE- * GRAL IS COMPLETE. IF R IS NEGATIVE, THE CAUCHY PRINCIPAL * VALUE IS COMPUTED BY USING A PRELIMINARY TRANSFORMATION * TO MAKE R POSITIVE; SEE EQUATION (2.22) OF THE SECOND REF- * ERENCE BELOW. WHEN R IS POSITIVE, 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 SEPT. 1, 1987. * * CHECK VALUES: RJ(2,3,4,5) = 0.14297 57966 71567 53833 23308 * RJ(2,3,4,-5) = - 0.12711 23004 29638 83590 80083 * CHECK BY ADDITION THEOREM: RJ(X,X+Z,X+W,X+P) * + RJ(Y,Y+Z,Y+W,Y+P) + (A-B) * RJ(A,B,B,A) + 3 / DSQRT(A) * = RJ(0,Z,W,P), WHERE X,Y,Z,W,P ARE POSITIVE AND X * Y * = Z * W, A = P * P * (X+Y+Z+W), B = P * (P+X) * (P+Y), * AND B - A = P * (P-Z) * (P-W). THE SUM OF THE THIRD AND * FOURTH TERMS ON THE LEFT SIDE IS 3 * RC(A,B). * * ***** Formal Arguments *********************************** * * INPUT ... * * X, Y, Z, AND R ARE THE VARIABLES IN THE INTEGRAL RJ(X,Y,Z,R). * * OUTPUT ... * * RJ 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, OR R = 0 * IERR = 2 X+Y, X+Z, Y+Z, OR ABS(R) IS TOO SMALL. * IERR = 3 X, Y, Z, OR ABS(R) IS TOO LARGE. * */ /*--S replaces "?": ?ERM1, ?ERV1, ?RCVAL, ?RFVLX, ?RJVAL * * ***** External References ******************************** * * SRCVAL COMPUTES RC. IT CHECKS THE ARGUMENTS, BUT IT'S NOT * CALLED FROM WITHIN THE INNER LOOP. THERE, WE INLINE SRCVAL. * SRFVLX COMPUTES RF, BUT ASSUMES ARGUMENTS ARE IN RANGE. * SRFVLX IS INSIDE SRFVAL. * */ /* ***** Local Variables ************************************ * */ /* MACHINE DEPENDENT PARAMETERS ... * * LOLIM AND UPLIM DETERMINE THE RANGE OF VALID ARGUMENTS. * LOLIM IS NOT LESS THAN THE CUBE ROOT OF THE VALUE * OF LOLIM USED IN THE CODE TO COMPUTE RC, AND * UPLIM IS NOT GREATER THAN 0.3 TIMES THE CUBE ROOT OF * THE VALUE OF UPLIM USED IN THE CODE TO COMPUTE RC. * * ERRTOL IS SET TO THE DESIRED ERROR TOLERANCE. THE * RELATIVE ERROR DUE TO TRUNCATION OF THE SERIES FOR RJ * IS LESS THAN 3 * ERRTOL ** 6 / (1 - ERRTOL) ** 3/2. * ERRTOL / 2 IS USED IN THE CODE TO COMPUTE RC TO MAKE * THE TRUNCATION ERROR FOR RC LESS THAN FOR RJ. * */ /* ---------------------------------------------------------------------- * WARNING. CHANGES IN THE PROGRAM MAY IMPROVE SPEED AT THE * EXPENSE OF ROBUSTNESS. * ---------------------------------------------------------------------- * */ if (fminf( x, fminf( y, z ) ) < 0.0e0 || r == 0.0e0) { *rj = 0.0e0; *ierr = 1; serm1( "SRJVAL", 1, 0, "ONE OF X, Y, or Z is negative, or R is zero" , "X", x, ',' ); serv1( "Y", y, ',' ); serv1( "Z", z, ',' ); serv1( "R", r, '.' ); return; } if (lolim < 0.0) { lolim = 1.0001e0*powf(5.0e0*FLT_MIN,THIRD); uplim = .29999e0*powf(FLT_MAX/5.0e0,THIRD); errtol = powf(0.28e0*FLT_EPSILON,1.0e0/6.0e0); etolrc = 0.5e0*errtol; } if (fmaxf( x, fmaxf( y, fmaxf( z, fabsf( r ) ) ) ) > uplim) { *rj = 0.0e0; *ierr = 3; serm1( "SRJVAL", 3, 0, "ONE OF X, Y, Z or ABS(R) > UPLIM", "X", x, ',' ); serv1( "Y", y, ',' ); serv1( "Z", z, ',' ); serv1( "R", r, ',' ); serv1( "UPLIM", uplim, '.' ); return; } if (fminf( x + y, fminf( x + z, fminf( y + z, fabsf( r ) ) ) ) < lolim) { *rj = 0.0e0; *ierr = 2; serm1( "SRJVAL", 2, 0, "ONE OF X+Y, X+Z, Y+Z or ABS(R) < LOLIM" , "X", x, ',' ); serv1( "Y", y, ',' ); serv1( "Z", z, ',' ); serv1( "R", r, ',' ); serv1( "LOLIM", lolim, '.' ); return; } *ierr = 0; if (r > 0.0e0) { xn = x; yn = y; zn = z; rn = r; } else { /* ORDER X,Y,Z AND TRANSFORM TO POSITIVE R */ xn = fminf( x, y ); yy = fmaxf( x, y ); zn = fmaxf( yy, z ); yy = fminf( yy, z ); yn = fmaxf( xn, yy ); xn = fminf( xn, yy ); a = 1.0e0/(yn - r); b = (zn - yn)*a*(yn - xn); rn = yn + b; alfa = xn*zn/yn; beta = r*rn/yn; srcval( alfa, beta, &rcx, ierr ); if (*ierr != 0) { *rj = 0.0e0; return; } } sigma = 0.0e0; power4 = 1.0e0; L_20: ; mu = (xn + yn + zn + rn + rn)*FIFTH; xndev = (mu - xn)/mu; yndev = (mu - yn)/mu; zndev = (mu - zn)/mu; rndev = (mu - rn)/mu; epslon = fmaxf( fabsf( xndev ), fmaxf( fabsf( yndev ), fmaxf( fabsf( zndev ), fabsf( rndev ) ) ) ); if (epslon < errtol) goto L_80; xnroot = sqrtf( xn ); ynroot = sqrtf( yn ); znroot = sqrtf( zn ); lamda = xnroot*(ynroot + znroot) + ynroot*znroot; alfa = rn*(xnroot + ynroot + znroot) + xnroot*ynroot*znroot; alfa *= alfa; beta = rn*(rn + lamda)*(rn + lamda); /* CALL SRCVAL (ALFA, BETA, RC, IERR) * IF (IERR .NE. 0) RETURN * We use the following instead of calling SRCVAL to avoid the * tests for argument range, which we know will succeed. */ L_40: ; delta = (alfa + beta + beta)*THIRD; sn = (beta + delta)/delta - 2.0e0; if (fabsf( sn ) < etolrc) goto L_60; gamma = 2.0*sqrtf( alfa )*sqrtf( beta ) + beta; alfa = (alfa + gamma)*0.25e0; beta = (beta + gamma)*0.25e0; goto L_40; L_60: ; rc = (1.0e0 + sn*sn*(Z1 + sn*(Z2 + sn*(Z3 + sn*Z4))))/sqrtf( delta ); sigma += power4*rc; power4 *= 0.25e0; xn = (xn + lamda)*0.25e0; yn = (yn + lamda)*0.25e0; zn = (zn + lamda)*0.25e0; rn = (rn + lamda)*0.25e0; goto L_20; L_80: ; e1 = yndev*zndev; ea = xndev*(yndev + zndev) + e1; eb = xndev*e1; ec = rndev*rndev; e2 = ea - 3.0e0*ec; e3 = eb + 2.0e0*rndev*(ea - ec); s1 = 1.0e0 + e2*(-C1 + 0.75e0*C3*e2 - 1.5e0*C4*e3); s2 = eb*(0.5e0*C2 + rndev*(-C3 - C3 + rndev*C4)); s3 = rndev*ea*(C2 - rndev*C3) - C2*rndev*ec; *rj = 3.0e0*sigma + power4*(s1 + s2 + s3)/(mu*sqrtf( mu )); if (r > 0.e0) return; srfvlx( xn, yn, zn, &rf ); *rj = a*(b**rj + 3.0e0*(rcx - rf)); return; } /* end of function */