/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:56 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "drjval.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*/ drjval( double x, double y, double z, double r, double *rj, long *ierr) { long int _l0; double 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 double errtol, etolrc, 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 DRJVAL Krogh Change -1.0 to -1.d0. *>> 1996-03-30 DRJVAL Krogh Added external statement. *>> 1995-11-17 DRJVAL Krogh Converted SFTRAN to Fortran 77. *>> 1994-10-19 DRJVAL Krogh Changes to use M77CON *>> 1994-08-15 DRJVAL WV Snyder JPL use 2-arg min and max for C conv. *>> 1992-04-17 DRJVAL WV Snyder JPL *>> 1992-04-13 DRJVAL WV Snyder JPL Give RJ a value even if ierr.ne.0 *>> 1992-04-08 DRJVAL WV Snyder JPL Declare RCX, spell YN correctly. *>> 1990-12-20 DRJVAL 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. * */ /*--D replaces "?": ?ERM1, ?ERV1, ?RCVAL, ?RFVLX, ?RJVAL * * ***** External References ******************************** * * DRCVAL COMPUTES RC. IT CHECKS THE ARGUMENTS, BUT IT'S NOT * CALLED FROM WITHIN THE INNER LOOP. THERE, WE INLINE DRCVAL. * DRFVLX COMPUTES RF, BUT ASSUMES ARGUMENTS ARE IN RANGE. * DRFVLX 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 (fmin( x, fmin( y, z ) ) < 0.0e0 || r == 0.0e0) { *rj = 0.0e0; *ierr = 1; derm1( "DRJVAL", 1, 0, "ONE OF X, Y, or Z is negative, or R is zero" , "X", x, ',' ); derv1( "Y", y, ',' ); derv1( "Z", z, ',' ); derv1( "R", r, '.' ); return; } if (lolim < 0.0) { lolim = 1.0001e0*pow(5.0e0*DBL_MIN,THIRD); uplim = .29999e0*pow(DBL_MAX/5.0e0,THIRD); errtol = pow(0.28e0*DBL_EPSILON,1.0e0/6.0e0); etolrc = 0.5e0*errtol; } if (fmax( x, fmax( y, fmax( z, fabs( r ) ) ) ) > uplim) { *rj = 0.0e0; *ierr = 3; derm1( "DRJVAL", 3, 0, "ONE OF X, Y, Z or ABS(R) > UPLIM", "X", x, ',' ); derv1( "Y", y, ',' ); derv1( "Z", z, ',' ); derv1( "R", r, ',' ); derv1( "UPLIM", uplim, '.' ); return; } if (fmin( x + y, fmin( x + z, fmin( y + z, fabs( r ) ) ) ) < lolim) { *rj = 0.0e0; *ierr = 2; derm1( "DRJVAL", 2, 0, "ONE OF X+Y, X+Z, Y+Z or ABS(R) < LOLIM" , "X", x, ',' ); derv1( "Y", y, ',' ); derv1( "Z", z, ',' ); derv1( "R", r, ',' ); derv1( "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 = fmin( x, y ); yy = fmax( x, y ); zn = fmax( yy, z ); yy = fmin( yy, z ); yn = fmax( xn, yy ); xn = fmin( xn, yy ); a = 1.0e0/(yn - r); b = (zn - yn)*a*(yn - xn); rn = yn + b; alfa = xn*zn/yn; beta = r*rn/yn; drcval( 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 = fmax( fabs( xndev ), fmax( fabs( yndev ), fmax( fabs( zndev ), fabs( rndev ) ) ) ); if (epslon < errtol) goto L_80; xnroot = sqrt( xn ); ynroot = sqrt( yn ); znroot = sqrt( zn ); lamda = xnroot*(ynroot + znroot) + ynroot*znroot; alfa = rn*(xnroot + ynroot + znroot) + xnroot*ynroot*znroot; alfa *= alfa; beta = rn*(rn + lamda)*(rn + lamda); /* CALL DRCVAL (ALFA, BETA, RC, IERR) * IF (IERR .NE. 0) RETURN * We use the following instead of calling DRCVAL 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 (fabs( sn ) < etolrc) goto L_60; gamma = 2.0*sqrt( alfa )*sqrt( 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))))/sqrt( 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*sqrt( mu )); if (r > 0.e0) return; drfvlx( xn, yn, zn, &rf ); *rj = a*(b**rj + 3.0e0*(rcx - rf)); return; } /* end of function */