FUNCTION DAW(XX) C---------------------------------------------------------------------- C C This function program evaluates Dawson's integral, C C 2 / x 2 C -x | t C F(x) = e | e dt C | C / 0 C C for a real argument x. C C The calling sequence for this function is C C Y=DAW(X) C C The main computation uses rational Chebyshev approximations C published in Math. Comp. 24, 171-178 (1970) by Cody, Paciorek C and Thacher. This transportable program is patterned after the C machine-dependent FUNPACK program DDAW(X), but cannot match that C version for efficiency or accuracy. This version uses rational C approximations that are theoretically accurate to about 19 C significant decimal digits. The accuracy achieved depends on the C arithmetic system, the compiler, the intrinsic functions, and C proper selection of the machine-dependent constants. C C******************************************************************* C******************************************************************* C C Explanation of machine-dependent constants C C XINF = largest positive machine number C XMIN = the smallest positive machine number. C EPS = smallest positive number such that 1+eps > 1. C Approximately beta**(-p), where beta is the machine C radix and p is the number of significant base-beta C digits in a floating-point number. C XMAX = absolute argument beyond which DAW(X) underflows. C XMAX = min(0.5/xmin, xinf). C XSMALL = absolute argument below DAW(X) may be represented C by X. We recommend XSMALL = sqrt(eps). C XLARGE = argument beyond which DAW(X) may be represented by C 1/(2x). We recommend XLARGE = 1/sqrt(eps). C C Approximate values for some important machines are C C beta p eps xmin xinf C C CDC 7600 (S.P.) 2 48 7.11E-15 3.14E-294 1.26E+322 C CRAY-1 (S.P.) 2 48 7.11E-15 4.58E-2467 5.45E+2465 C IEEE (IBM/XT, C SUN, etc.) (S.P.) 2 24 1.19E-07 1.18E-38 3.40E+38 C IEEE (IBM/XT, C SUN, etc.) (D.P.) 2 53 1.11D-16 2.23E-308 1.79D+308 C IBM 3033 (D.P.) 16 14 1.11D-16 5.40D-79 7.23D+75 C VAX 11/780 (S.P.) 2 24 5.96E-08 2.94E-39 1.70E+38 C (D.P.) 2 56 1.39D-17 2.94D-39 1.70D+38 C (G Format) (D.P.) 2 53 1.11D-16 5.57D-309 8.98D+307 C C XSMALL XLARGE XMAX C C CDC 7600 (S.P.) 5.96E-08 1.68E+07 1.59E+293 C CRAY-1 (S.P.) 5.96E-08 1.68E+07 5.65E+2465 C IEEE (IBM/XT, C SUN, etc.) (S.P.) 2.44E-04 4.10E+03 4.25E+37 C IEEE (IBM/XT, C SUN, etc.) (D.P.) 1.05E-08 9.49E+07 2.24E+307 C IBM 3033 (D.P.) 3.73D-09 2.68E+08 7.23E+75 C VAX 11/780 (S.P.) 2.44E-04 4.10E+03 1.70E+38 C (D.P.) 3.73E-09 2.68E+08 1.70E+38 C (G Format) (D.P.) 1.05E-08 9.49E+07 8.98E+307 C C******************************************************************* C******************************************************************* C C Error Returns C C The program returns 0.0 for |X| > XMAX. C C Intrinsic functions required are: C C ABS C C C Author: W. J. Cody C Mathematics and Computer Science Division C Argonne National Laboratory C Argonne, IL 60439 C C Latest modification: June 15, 1988 C C---------------------------------------------------------------------- INTEGER I CS REAL CD DOUBLE PRECISION 1 DAW,FRAC,HALF,ONE,ONE225,P1,P2,P3,P4,Q1,Q2,Q3,Q4,SIX25, 2 SUMP,SUMQ,TWO5,W2,X,XX,Y,XLARGE,XMAX,XSMALL,ZERO DIMENSION P1(10),P2(10),P3(10),P4(10),Q1(10),Q2(9),Q3(9),Q4(9) C---------------------------------------------------------------------- C Mathematical constants. C---------------------------------------------------------------------- CS DATA ZERO,HALF,ONE/0.0E0,0.5E0,1.0E0/, CS 1 SIX25,ONE225,TWO5/6.25E0,12.25E0,25.0E0/ CD DATA ZERO,HALF,ONE/0.0D0,0.5D0,1.0D0/, CD 1 SIX25,ONE225,TWO5/6.25D0,12.25D0,25.0D0/ C---------------------------------------------------------------------- C Machine-dependent constants C---------------------------------------------------------------------- CS DATA XSMALL/2.44E-04/, XLARGE/4.10E+03/, XMAX/4.25E+37/ CD DATA XSMALL/1.05D-08/, XLARGE/9.49D+07/, XMAX/2.24D+307/ C---------------------------------------------------------------------- C Coefficients for R(9,9) approximation for |x| < 2.5 C---------------------------------------------------------------------- CS DATA P1/-2.69020398788704782410E-12, 4.18572065374337710778E-10, CS 1 -1.34848304455939419963E-08, 9.28264872583444852976E-07, CS 2 -1.23877783329049120592E-05, 4.07205792429155826266E-04, CS 3 -2.84388121441008500446E-03, 4.70139022887204722217E-02, CS 4 -1.38868086253931995101E-01, 1.00000000000000000004E+00/ CS DATA Q1/ 1.71257170854690554214E-10, 1.19266846372297253797E-08, CS 1 4.32287827678631772231E-07, 1.03867633767414421898E-05, CS 2 1.78910965284246249340E-04, 2.26061077235076703171E-03, CS 3 2.07422774641447644725E-02, 1.32212955897210128811E-01, CS 4 5.27798580412734677256E-01, 1.00000000000000000000E+00/ CD DATA P1/-2.69020398788704782410D-12, 4.18572065374337710778D-10, CD 1 -1.34848304455939419963D-08, 9.28264872583444852976D-07, CD 2 -1.23877783329049120592D-05, 4.07205792429155826266D-04, CD 3 -2.84388121441008500446D-03, 4.70139022887204722217D-02, CD 4 -1.38868086253931995101D-01, 1.00000000000000000004D+00/ CD DATA Q1/ 1.71257170854690554214D-10, 1.19266846372297253797D-08, CD 1 4.32287827678631772231D-07, 1.03867633767414421898D-05, CD 2 1.78910965284246249340D-04, 2.26061077235076703171D-03, CD 3 2.07422774641447644725D-02, 1.32212955897210128811D-01, CD 4 5.27798580412734677256D-01, 1.00000000000000000000D+00/ C---------------------------------------------------------------------- C Coefficients for R(9,9) approximation in J-fraction form C for x in [2.5, 3.5) C---------------------------------------------------------------------- CS DATA P2/-1.70953804700855494930E+00,-3.79258977271042880786E+01, CS 1 2.61935631268825992835E+01, 1.25808703738951251885E+01, CS 2 -2.27571829525075891337E+01, 4.56604250725163310122E+00, CS 3 -7.33080089896402870750E+00, 4.65842087940015295573E+01, CS 4 -1.73717177843672791149E+01, 5.00260183622027967838E-01/ CS DATA Q2/ 1.82180093313514478378E+00, 1.10067081034515532891E+03, CS 1 -7.08465686676573000364E+00, 4.53642111102577727153E+02, CS 2 4.06209742218935689922E+01, 3.02890110610122663923E+02, CS 3 1.70641269745236227356E+02, 9.51190923960381458747E+02, CS 4 2.06522691539642105009E-01/ CD DATA P2/-1.70953804700855494930D+00,-3.79258977271042880786D+01, CD 1 2.61935631268825992835D+01, 1.25808703738951251885D+01, CD 2 -2.27571829525075891337D+01, 4.56604250725163310122D+00, CD 3 -7.33080089896402870750D+00, 4.65842087940015295573D+01, CD 4 -1.73717177843672791149D+01, 5.00260183622027967838D-01/ CD DATA Q2/ 1.82180093313514478378D+00, 1.10067081034515532891D+03, CD 1 -7.08465686676573000364D+00, 4.53642111102577727153D+02, CD 2 4.06209742218935689922D+01, 3.02890110610122663923D+02, CD 3 1.70641269745236227356D+02, 9.51190923960381458747D+02, CD 4 2.06522691539642105009D-01/ C---------------------------------------------------------------------- C Coefficients for R(9,9) approximation in J-fraction form C for x in [3.5, 5.0] C---------------------------------------------------------------------- CS DATA P3/-4.55169503255094815112E+00,-1.86647123338493852582E+01, CS 1 -7.36315669126830526754E+00,-6.68407240337696756838E+01, CS 2 4.84507265081491452130E+01, 2.69790586735467649969E+01, CS 3 -3.35044149820592449072E+01, 7.50964459838919612289E+00, CS 4 -1.48432341823343965307E+00, 4.99999810924858824981E-01/ CS DATA Q3/ 4.47820908025971749852E+01, 9.98607198039452081913E+01, CS 1 1.40238373126149385228E+01, 3.48817758822286353588E+03, CS 2 -9.18871385293215873406E+00, 1.24018500009917163023E+03, CS 3 -6.88024952504512254535E+01,-2.31251575385145143070E+00, CS 4 2.50041492369922381761E-01/ CD DATA P3/-4.55169503255094815112D+00,-1.86647123338493852582D+01, CD 1 -7.36315669126830526754D+00,-6.68407240337696756838D+01, CD 2 4.84507265081491452130D+01, 2.69790586735467649969D+01, CD 3 -3.35044149820592449072D+01, 7.50964459838919612289D+00, CD 4 -1.48432341823343965307D+00, 4.99999810924858824981D-01/ CD DATA Q3/ 4.47820908025971749852D+01, 9.98607198039452081913D+01, CD 1 1.40238373126149385228D+01, 3.48817758822286353588D+03, CD 2 -9.18871385293215873406D+00, 1.24018500009917163023D+03, CD 3 -6.88024952504512254535D+01,-2.31251575385145143070D+00, CD 4 2.50041492369922381761D-01/ C---------------------------------------------------------------------- C Coefficients for R(9,9) approximation in J-fraction form C for |x| > 5.0 C---------------------------------------------------------------------- CS DATA P4/-8.11753647558432685797E+00,-3.84043882477454453430E+01, CS 1 -2.23787669028751886675E+01,-2.88301992467056105854E+01, CS 2 -5.99085540418222002197E+00,-1.13867365736066102577E+01, CS 3 -6.52828727526980741590E+00,-4.50002293000355585708E+00, CS 4 -2.50000000088955834952E+00, 5.00000000000000488400E-01/ CS DATA Q4/ 2.69382300417238816428E+02, 5.04198958742465752861E+01, CS 1 6.11539671480115846173E+01, 2.08210246935564547889E+02, CS 2 1.97325365692316183531E+01,-1.22097010558934838708E+01, CS 3 -6.99732735041547247161E+00,-2.49999970104184464568E+00, CS 4 7.49999999999027092188E-01/ CD DATA P4/-8.11753647558432685797D+00,-3.84043882477454453430D+01, CD 1 -2.23787669028751886675D+01,-2.88301992467056105854D+01, CD 2 -5.99085540418222002197D+00,-1.13867365736066102577D+01, CD 3 -6.52828727526980741590D+00,-4.50002293000355585708D+00, CD 4 -2.50000000088955834952D+00, 5.00000000000000488400D-01/ CD DATA Q4/ 2.69382300417238816428D+02, 5.04198958742465752861D+01, CD 1 6.11539671480115846173D+01, 2.08210246935564547889D+02, CD 2 1.97325365692316183531D+01,-1.22097010558934838708D+01, CD 3 -6.99732735041547247161D+00,-2.49999970104184464568D+00, CD 4 7.49999999999027092188D-01/ C---------------------------------------------------------------------- X = XX IF (ABS(X) .GT. XLARGE) THEN IF (ABS(X) .LE. XMAX) THEN DAW = HALF / X ELSE DAW = ZERO END IF ELSE IF (ABS(X) .LT. XSMALL) THEN DAW = X ELSE Y = X * X IF (Y .LT. SIX25) THEN C---------------------------------------------------------------------- C ABS(X) .LT. 2.5 C---------------------------------------------------------------------- SUMP = P1(1) SUMQ = Q1(1) DO 100 I = 2, 10 SUMP = SUMP * Y + P1(I) SUMQ = SUMQ * Y + Q1(I) 100 CONTINUE DAW = X * SUMP / SUMQ ELSE IF (Y .LT. ONE225) THEN C---------------------------------------------------------------------- C 2.5 .LE. ABS(X) .LT. 3.5 C---------------------------------------------------------------------- FRAC = ZERO DO 200 I = 1, 9 200 FRAC = Q2(I) / (P2(I) + Y + FRAC) DAW = (P2(10) + FRAC) / X ELSE IF (Y .LT. TWO5) THEN C---------------------------------------------------------------------- C 3.5 .LE. ABS(X) .LT. 5.0 C--------------------------------------------------------------------- FRAC = ZERO DO 300 I = 1, 9 300 FRAC = Q3(I) / (P3(I) + Y + FRAC) DAW = (P3(10) + FRAC) / X ELSE C---------------------------------------------------------------------- C 5.0 .LE. ABS(X) .LE. XLARGE C------------------------------------------------------------------ W2 = ONE / X / X FRAC = ZERO DO 400 I = 1, 9 400 FRAC = Q4(I) / (P4(I) + Y + FRAC) FRAC = P4(10) + FRAC DAW = (HALF + HALF * W2 * FRAC) / X END IF END IF RETURN C---------- Last line of DAW ---------- END