1 SUBROUTINE DRCVAL (X, Y, RC, IERR) 2c Copyright (c) 1996 California Institute of Technology, Pasadena, CA. 3c ALL RIGHTS RESERVED. 4c Based on Government Sponsored Research NAS7-03001. 5c>> 2001-07-16 DRCVAL Krogh Change -1.0 to -1.d0. 6c>> 1996-03-30 DRCVAL Krogh Added external statement. 7C>> 1995-11-17 DRCVAL Krogh Converted SFTRAN to Fortran 77. 8C>> 1994-10-19 DRCVAL Krogh Changes to use M77CON 9C>> 1991-10-31 DRCVAL WV Snyder JPL Incorporate changes from Carlson 10C>> 1990-12-20 DRCVAL WV Snyder JPL Convert from NSWC for Math 77. 11C 12C THIS SUBROUTINE COMPUTES THE ELEMENTARY INTEGRAL 13C 14C RC(X,Y) = INTEGRAL FROM ZERO TO INFINITY OF 15C 16C -1/2 -1 17C (1/2)(T+X) (T+Y) DT, 18C 19C WHERE X IS NONNEGATIVE AND Y IS NONZERO. IF Y IS NEGATIVE, 20C THE CAUCHY PRINCIPAL VALUE IS COMPUTED BY USING A PRELIMI- 21C NARY TRANSFORMATION TO MAKE Y POSITIVE; SEE EQUATION (2.12) 22C OF THE SECOND REFERENCE BELOW. WHEN Y IS POSITIVE, THE 23C DUPULICATION THEOREM IS ITERATED UNTIL THE VARIABLES ARE 24C NEARLY EQUAL, AND THE FUNCTION IS THEN EXPANDED IN TAYLOR 25C SERIES TO FIFTH ORDER. LOGARITHMIC, INVERSE CIRCULAR, AND 26C INVERSE HYPERBOLIC FUNCTIONS ARE EXPRESSED IN TERMS OF RC 27C BY EQUATIONS (4.9)-(4.13) OF THE SECOND REFERENCE BELOW. 28C 29C REFERENCES: B. C. CARLSON AND E. M. NOTIS, ALGORITHMS FOR 30C INCOMPLETE ELLIPTIC INTEGRALS, ACM TRANSACTIONS ON MATHEMA- 31C TICAL SOFTWARE, 7 (1981), 398-403; 32C B. C. CARLSON, COMPUTING ELLIPTIC INTEGRALS BY DUPLICATION, 33C NUMER. MATH. 33 (1979), 1-16. 34C AUTHORS: B. C. CARLSON AND ELAINE M. NOTIS, AMES LABORATORY- 35C DOE, IOWA STATE UNIVERSITY, AMES, IA 50011, AND R. L. PEXTON, 36C LAWRENCE LIVERMORE NATIONAL LABORATORY, LIVERMORE, CA 94550. 37C AUG. 1, 1979, REVISED SEPT. 1, 1987. 38C 39C CHECK VALUES: RC(0,1/4) = RC(1/16,1/8) = PI, 40C RC(9/4,2) = LN(2), 41C RC(1/4,-2) = LN(2)/3. 42C CHECK BY ADDITION THEOREM: RC(X,X+Z) + RC(Y,Y+Z) = RC(0,Z), 43C WHERE X, Y, AND Z ARE POSITIVE AND X * Y = Z * Z. 44C 45C ***** Formal Arguments *********************************** 46C 47C LOLIM AND UPLIM DETERMINE THE RANGE OF VALID ARGUMENTS. 48C LOLIM IS NOT LESS THAN THE MACHINE MINIMUM MULTIPLIED BY 5. 49C UPLIM IS NOT GREATER THAN THE MACHINE MAXIMUM DIVIDED BY 5. 50C 51C INPUT ... 52C 53C X AND Y ARE THE VARIABLES IN THE INTEGRAL RC(X,Y). 54C 55C OUTPUT ... 56C 57C RC IS THE VALUE OF THE INTEGRAL. 58C 59C IERR IS THE RETURN ERROR CODE. 60C 0 FOR NORMAL COMPLETION OF THE SUBROUTINE. 61C 1 X IS NEGATIVE, OR Y = 0. 62C 2 X+ABS(Y) IS TOO SMALL. 63C 3 X OR ABS(Y) IS TOO LARGE, OR X + ABS(Y) IS TOO LARGE. 64C 4 Y < -2.236/SQRT(LOLIM) AND 0 < X < (LOLIM*UPLIM)**2/25 65C 66 DOUBLE PRECISION X, Y, RC 67 INTEGER IERR 68C 69C--D replaces "?": ?ERM1, ?ERV1, ?RCVAL 70C 71C ***** External References ******************************** 72C 73 EXTERNAL D1MACH 74 DOUBLE PRECISION D1MACH 75C 76C ***** Local Variables ************************************ 77C 78 DOUBLE PRECISION C1,C2,C3,C4,ERRTOL,LAMDA,LOLIM 79 DOUBLE PRECISION MU,SN,UPLIM,XLU225,XN,Y2236L,YN,W 80 PARAMETER (C1 = 3.0d0 / 10.0d0) 81 PARAMETER (C2 = 1.0d0 / 7.0d0) 82 PARAMETER (C3 = 3.0d0 / 8.0d0) 83 PARAMETER (C4 = 9.0d0 / 22.0d0) 84 SAVE ERRTOL, LOLIM, UPLIM, XLU225, Y2236L 85 DATA LOLIM /-1.0d0/ 86C 87C ERRTOL IS SET TO THE DESIRED ERROR TOLERANCE. 88C RELATIVE ERROR DUE TO TRUNCATION IS LESS THAN 89C 16 * ERRTOL ** 6 / (1 - 2 * ERRTOL). 90C 91C SAMPLE CHOICES ERRTOL RELATIVE TRUNCATION 92C ERROR LESS THAN 93C 1.D-3 2.D-17 94C 3.D-3 2.D-14 95C 1.D-2 2.D-11 96C 3.D-2 2.D-8 97C 1.D-1 2.D-5 98C 99C We could put in a Newton iteration to solve for ERRTOL in 100C terms of D1MACH(4), but it seems good enough to put 101C ERRTOL = (D1MACH(4)/16.0)**(1.0/6.0) 102C 103c ---------------------------------------------------------------------- 104C WARNING. CHANGES IN THE PROGRAM MAY IMPROVE SPEED AT THE 105C EXPENSE OF ROBUSTNESS. 106c ---------------------------------------------------------------------- 107C 108 IF (X .LT. 0.0d0 .OR. Y .EQ. 0.0d0) THEN 109 RC = 0.0d0 110 CALL DERM1 ('DRCVAL',1,0,'X < 0 or Y = 0','X',X,',') 111 CALL DERV1 ('Y',Y,'.') 112 IERR = 1 113 RETURN 114 END IF 115 IF (LOLIM .LT. 0.0d0) THEN 116 LOLIM = 5.0d0 * D1MACH(1) 117 UPLIM = D1MACH(2) / 5.0d0 118 XLU225 = (LOLIM*UPLIM)**2 / 25.0d0 119 Y2236L = -2.236d0/SQRT(LOLIM) 120 ERRTOL = (D1MACH(4)/16.0d0) ** (1.0d0/6.0d0) 121 END IF 122 YN = ABS(Y) 123 IF (X .GT. UPLIM .OR. YN .GT. UPLIM) THEN 124 RC = 0.0d0 125 CALL DERM1 ('DRCVAL',3,0,'X > UPLIM or ABS(Y) > UPLIM','X',X,',') 126 CALL DERV1 ('Y',Y,',') 127 CALL DERV1 ('UPLIM',UPLIM,'.') 128 IERR = 3 129 RETURN 130 END IF 131 IF (X + YN .LT. LOLIM) THEN 132 RC = 0.0d0 133 CALL DERM1 ('DRCVAL',2,0,'X + ABS(Y) < LOLIM','X',X,',') 134 CALL DERV1 ('Y',Y,',') 135 CALL DERV1 ('LOLIM',LOLIM,'.') 136 IERR = 2 137 RETURN 138 END IF 139 IF (X + YN .GT. UPLIM) THEN 140 RC = 0.0d0 141 CALL DERM1 ('DRCVAL',3,0,'X + ABS(Y) > UPLIM','X',X,',') 142 CALL DERV1 ('Y',Y,',') 143 CALL DERV1 ('UPLIM',UPLIM,'.') 144 IERR = 3 145 RETURN 146 END IF 147 IF (Y.LT.Y2236L .AND. X.GT.0.0d0 .AND. X.LT.XLU225) THEN 148 RC = 0.0d0 149 CALL DERM1 ('DRCVAL',4,0, 150 1 'Y < -2.236/SQRT(LOLIM) AND 0 < X < (LOLIM*UPLIM)**2/25', 151 2 'X',X,',') 152 CALL DERV1 ('Y',Y,',') 153 CALL DERV1 ('LOLIM',LOLIM,',') 154 CALL DERV1 ('UPLIM',UPLIM,'.') 155 IERR = 4 156 RETURN 157 END IF 158C 159 IERR = 0 160 IF (Y.GT.0.0d0) THEN 161 XN = X 162 W = 1.0d0 163 ELSE 164C TRANSFORM TO POSITIVE Y 165 XN = X - Y 166 W = SQRT(X) / SQRT(XN) 167 END IF 168C 169 20 CONTINUE 170 MU = (XN + YN + YN) / 3.0d0 171 SN = (YN + MU) / MU - 2.0d0 172 IF (ABS(SN) .LT. ERRTOL) THEN 173 RC = W * (1.0d0+SN*SN*(C1+SN*(C2+SN*(C3+SN*C4)))) / SQRT(MU) 174 RETURN 175 END IF 176 LAMDA = 2.0 * SQRT(XN) * SQRT(YN) + YN 177 XN = (XN + LAMDA) * 0.25d0 178 YN = (YN + LAMDA) * 0.25d0 179 GO TO 20 180 END 181