1*DECK DRF 2 DOUBLE PRECISION FUNCTION DRF (X, Y, Z, IER) 3C***BEGIN PROLOGUE DRF 4C***PURPOSE Compute the incomplete or complete elliptic integral of the 5C 1st kind. For X, Y, and Z non-negative and at most one of 6C them zero, RF(X,Y,Z) = Integral from zero to infinity of 7C -1/2 -1/2 -1/2 8C (1/2)(t+X) (t+Y) (t+Z) dt. 9C If X, Y or Z is zero, the integral is complete. 10C***LIBRARY SLATEC 11C***CATEGORY C14 12C***TYPE DOUBLE PRECISION (RF-S, DRF-D) 13C***KEYWORDS COMPLETE ELLIPTIC INTEGRAL, DUPLICATION THEOREM, 14C INCOMPLETE ELLIPTIC INTEGRAL, INTEGRAL OF THE FIRST KIND, 15C TAYLOR SERIES 16C***AUTHOR Carlson, B. C. 17C Ames Laboratory-DOE 18C Iowa State University 19C Ames, IA 50011 20C Notis, E. M. 21C Ames Laboratory-DOE 22C Iowa State University 23C Ames, IA 50011 24C Pexton, R. L. 25C Lawrence Livermore National Laboratory 26C Livermore, CA 94550 27C***DESCRIPTION 28C 29C 1. DRF 30C Evaluate an INCOMPLETE (or COMPLETE) ELLIPTIC INTEGRAL 31C of the first kind 32C Standard FORTRAN function routine 33C Double precision version 34C The routine calculates an approximation result to 35C DRF(X,Y,Z) = Integral from zero to infinity of 36C 37C -1/2 -1/2 -1/2 38C (1/2)(t+X) (t+Y) (t+Z) dt, 39C 40C where X, Y, and Z are nonnegative and at most one of them 41C is zero. If one of them is zero, the integral is COMPLETE. 42C The duplication theorem is iterated until the variables are 43C nearly equal, and the function is then expanded in Taylor 44C series to fifth order. 45C 46C 2. Calling sequence 47C DRF( X, Y, Z, IER ) 48C 49C Parameters On entry 50C Values assigned by the calling routine 51C 52C X - Double precision, nonnegative variable 53C 54C Y - Double precision, nonnegative variable 55C 56C Z - Double precision, nonnegative variable 57C 58C 59C 60C On Return (values assigned by the DRF routine) 61C 62C DRF - Double precision approximation to the integral 63C 64C IER - Integer 65C 66C IER = 0 Normal and reliable termination of the 67C routine. It is assumed that the requested 68C accuracy has been achieved. 69C 70C IER > 0 Abnormal termination of the routine 71C 72C X, Y, Z are unaltered. 73C 74C 75C 3. Error Messages 76C 77C 78C Value of IER assigned by the DRF routine 79C 80C Value assigned Error Message Printed 81C IER = 1 MIN(X,Y,Z) .LT. 0.0D0 82C = 2 MIN(X+Y,X+Z,Y+Z) .LT. LOLIM 83C = 3 MAX(X,Y,Z) .GT. UPLIM 84C 85C 86C 87C 4. Control Parameters 88C 89C Values of LOLIM, UPLIM, and ERRTOL are set by the 90C routine. 91C 92C LOLIM and UPLIM determine the valid range of X, Y and Z 93C 94C LOLIM - Lower limit of valid arguments 95C 96C Not less than 5 * (machine minimum). 97C 98C UPLIM - Upper limit of valid arguments 99C 100C Not greater than (machine maximum) / 5. 101C 102C 103C Acceptable values for: LOLIM UPLIM 104C IBM 360/370 SERIES : 3.0D-78 1.0D+75 105C CDC 6000/7000 SERIES : 1.0D-292 1.0D+321 106C UNIVAC 1100 SERIES : 1.0D-307 1.0D+307 107C CRAY : 2.3D-2466 1.09D+2465 108C VAX 11 SERIES : 1.5D-38 3.0D+37 109C 110C 111C 112C ERRTOL determines the accuracy of the answer 113C 114C The value assigned by the routine will result 115C in solution precision within 1-2 decimals of 116C "machine precision". 117C 118C 119C 120C ERRTOL - Relative error due to truncation is less than 121C ERRTOL ** 6 / (4 * (1-ERRTOL) . 122C 123C 124C 125C The accuracy of the computed approximation to the integral 126C can be controlled by choosing the value of ERRTOL. 127C Truncation of a Taylor series after terms of fifth order 128C introduces an error less than the amount shown in the 129C second column of the following table for each value of 130C ERRTOL in the first column. In addition to the truncation 131C error there will be round-off error, but in practice the 132C total error from both sources is usually less than the 133C amount given in the table. 134C 135C 136C 137C 138C 139C Sample choices: ERRTOL Relative Truncation 140C error less than 141C 1.0D-3 3.0D-19 142C 3.0D-3 2.0D-16 143C 1.0D-2 3.0D-13 144C 3.0D-2 2.0D-10 145C 1.0D-1 3.0D-7 146C 147C 148C Decreasing ERRTOL by a factor of 10 yields six more 149C decimal digits of accuracy at the expense of one or 150C two more iterations of the duplication theorem. 151C 152C *Long Description: 153C 154C DRF Special Comments 155C 156C 157C 158C Check by addition theorem: DRF(X,X+Z,X+W) + DRF(Y,Y+Z,Y+W) 159C = DRF(0,Z,W), where X,Y,Z,W are positive and X * Y = Z * W. 160C 161C 162C On Input: 163C 164C X, Y, and Z are the variables in the integral DRF(X,Y,Z). 165C 166C 167C On Output: 168C 169C 170C X, Y, Z are unaltered. 171C 172C 173C 174C ******************************************************** 175C 176C WARNING: Changes in the program may improve speed at the 177C expense of robustness. 178C 179C 180C 181C Special double precision functions via DRF 182C 183C 184C 185C 186C Legendre form of ELLIPTIC INTEGRAL of 1st kind 187C 188C ----------------------------------------- 189C 190C 191C 192C 2 2 2 193C F(PHI,K) = SIN(PHI) DRF(COS (PHI),1-K SIN (PHI),1) 194C 195C 196C 2 197C K(K) = DRF(0,1-K ,1) 198C 199C 200C PI/2 2 2 -1/2 201C = INT (1-K SIN (PHI) ) D PHI 202C 0 203C 204C 205C 206C Bulirsch form of ELLIPTIC INTEGRAL of 1st kind 207C 208C ----------------------------------------- 209C 210C 211C 2 2 2 212C EL1(X,KC) = X DRF(1,1+KC X ,1+X ) 213C 214C 215C Lemniscate constant A 216C 217C ----------------------------------------- 218C 219C 220C 1 4 -1/2 221C A = INT (1-S ) DS = DRF(0,1,2) = DRF(0,2,1) 222C 0 223C 224C 225C 226C ------------------------------------------------------------------- 227C 228C***REFERENCES B. C. Carlson and E. M. Notis, Algorithms for incomplete 229C elliptic integrals, ACM Transactions on Mathematical 230C Software 7, 3 (September 1981), pp. 398-403. 231C B. C. Carlson, Computing elliptic integrals by 232C duplication, Numerische Mathematik 33, (1979), 233C pp. 1-16. 234C B. C. Carlson, Elliptic integrals of the first kind, 235C SIAM Journal of Mathematical Analysis 8, (1977), 236C pp. 231-242. 237C***ROUTINES CALLED D1MACH, XERMSG 238C***REVISION HISTORY (YYMMDD) 239C 790801 DATE WRITTEN 240C 890531 Changed all specific intrinsics to generic. (WRB) 241C 891009 Removed unreferenced statement labels. (WRB) 242C 891009 REVISION DATE from Version 3.2 243C 891214 Prologue converted to Version 4.0 format. (BAB) 244C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) 245C 900326 Removed duplicate information from DESCRIPTION section. 246C (WRB) 247C 900510 Changed calls to XERMSG to standard form, and some 248C editorial changes. (RWC)) 249C 920501 Reformatted the REFERENCES section. (WRB) 250C***END PROLOGUE DRF 251 CHARACTER*16 XERN3, XERN4, XERN5, XERN6 252 INTEGER IER 253 DOUBLE PRECISION LOLIM, UPLIM, EPSLON, ERRTOL, D1MACH 254 DOUBLE PRECISION C1, C2, C3, E2, E3, LAMDA 255 DOUBLE PRECISION MU, S, X, XN, XNDEV 256 DOUBLE PRECISION XNROOT, Y, YN, YNDEV, YNROOT, Z, ZN, ZNDEV, 257 * ZNROOT 258 LOGICAL FIRST 259 SAVE ERRTOL,LOLIM,UPLIM,C1,C2,C3,FIRST 260 DATA FIRST /.TRUE./ 261C 262C***FIRST EXECUTABLE STATEMENT DRF 263C 264 IF (FIRST) THEN 265 ERRTOL = (4.0D0*D1MACH(3))**(1.0D0/6.0D0) 266 LOLIM = 5.0D0 * D1MACH(1) 267 UPLIM = D1MACH(2)/5.0D0 268C 269 C1 = 1.0D0/24.0D0 270 C2 = 3.0D0/44.0D0 271 C3 = 1.0D0/14.0D0 272 ENDIF 273 FIRST = .FALSE. 274C 275C CALL ERROR HANDLER IF NECESSARY. 276C 277 DRF = 0.0D0 278 IF (MIN(X,Y,Z).LT.0.0D0) THEN 279 IER = 1 280 WRITE (XERN3, '(1PE15.6)') X 281 WRITE (XERN4, '(1PE15.6)') Y 282 WRITE (XERN5, '(1PE15.6)') Z 283 CALL XERMSG ('SLATEC', 'DRF', 284 * 'MIN(X,Y,Z).LT.0 WHERE X = ' // XERN3 // ' Y = ' // XERN4 // 285 * ' AND Z = ' // XERN5, 1, 1) 286 RETURN 287 ENDIF 288C 289 IF (MAX(X,Y,Z).GT.UPLIM) THEN 290 IER = 3 291 WRITE (XERN3, '(1PE15.6)') X 292 WRITE (XERN4, '(1PE15.6)') Y 293 WRITE (XERN5, '(1PE15.6)') Z 294 WRITE (XERN6, '(1PE15.6)') UPLIM 295 CALL XERMSG ('SLATEC', 'DRF', 296 * 'MAX(X,Y,Z).GT.UPLIM WHERE X = ' // XERN3 // ' Y = ' // 297 * XERN4 // ' Z = ' // XERN5 // ' AND UPLIM = ' // XERN6, 3, 1) 298 RETURN 299 ENDIF 300C 301 IF (MIN(X+Y,X+Z,Y+Z).LT.LOLIM) THEN 302 IER = 2 303 WRITE (XERN3, '(1PE15.6)') X 304 WRITE (XERN4, '(1PE15.6)') Y 305 WRITE (XERN5, '(1PE15.6)') Z 306 WRITE (XERN6, '(1PE15.6)') LOLIM 307 CALL XERMSG ('SLATEC', 'DRF', 308 * 'MIN(X+Y,X+Z,Y+Z).LT.LOLIM WHERE X = ' // XERN3 // 309 * ' Y = ' // XERN4 // ' Z = ' // XERN5 // ' AND LOLIM = ' // 310 * XERN6, 2, 1) 311 RETURN 312 ENDIF 313C 314 IER = 0 315 XN = X 316 YN = Y 317 ZN = Z 318C 319 30 MU = (XN+YN+ZN)/3.0D0 320 XNDEV = 2.0D0 - (MU+XN)/MU 321 YNDEV = 2.0D0 - (MU+YN)/MU 322 ZNDEV = 2.0D0 - (MU+ZN)/MU 323 EPSLON = MAX(ABS(XNDEV),ABS(YNDEV),ABS(ZNDEV)) 324 IF (EPSLON.LT.ERRTOL) GO TO 40 325 XNROOT = SQRT(XN) 326 YNROOT = SQRT(YN) 327 ZNROOT = SQRT(ZN) 328 LAMDA = XNROOT*(YNROOT+ZNROOT) + YNROOT*ZNROOT 329 XN = (XN+LAMDA)*0.250D0 330 YN = (YN+LAMDA)*0.250D0 331 ZN = (ZN+LAMDA)*0.250D0 332 GO TO 30 333C 334 40 E2 = XNDEV*YNDEV - ZNDEV*ZNDEV 335 E3 = XNDEV*YNDEV*ZNDEV 336 S = 1.0D0 + (C1*E2-0.10D0-C2*E3)*E2 + C3*E3 337 DRF = S/SQRT(MU) 338C 339 RETURN 340 END 341