Blame view
fvn_fnlib/dspenc.f
5.27 KB
38581db0c
git-svn-id: https...
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 |
double precision function dspenc (x)
c feb 1978 edition. w. fullerton, c3, los alamos scientific lab.
c
c evaluate a form of spence-s function defined by
c integral from 0 to x of -alog(abs(1-y))/y dy.
c for abs(x).le.1, the uniformly convergent expansion
c spenc = sum k=1,infinity x**k / k**2 is valid.
c
c spence-s function can be used to evaluate much more general integral
c forms. for example,
c integral from 0 to z of alog(a*x+b)/(c*x+d) dx =
c alog(abs(b-a*d/c))*alog(abs(a*(c*z+d)/(a*d-b*c)))/c
c - dspenc (a*(c*z+d)/(a*d-b*c)) / c.
c
c ref -- k. mitchell, philosophical magazine, 40, p.351 (1949).
c stegun and abramowitz, ams 55, p.1004.
c
double precision x, spencs(38), aln, pi26, xbig, d1mach, dcsevl,
1 dlog
external d1mach, dcsevl, initds
c
c series for spen on the interval 0. to 5.00000e-01
c with weighted error 4.74e-32
c log weighted error 31.32
c significant figures required 30.37
c decimal places required 32.11
c
data spencs( 1) / +.1527365598 8924058729 4668491002 8 d+0 /
data spencs( 2) / +.8169658058 0510144035 0183818527 1 d-1 /
data spencs( 3) / +.5814157140 7787308729 7735064118 2 d-2 /
data spencs( 4) / +.5371619814 5415275422 4788900531 9 d-3 /
data spencs( 5) / +.5724704675 1858262332 1060305478 2 d-4 /
data spencs( 6) / +.6674546121 6493363436 0783543858 9 d-5 /
data spencs( 7) / +.8276467339 7156769815 8439168901 1 d-6 /
data spencs( 8) / +.1073315673 0306789512 7000587335 4 d-6 /
data spencs( 9) / +.1440077294 3032394023 3459033151 3 d-7 /
data spencs( 10) / +.1984442029 9659063678 9887713960 8 d-8 /
data spencs( 11) / +.2794005822 1636387202 0199482161 5 d-9 /
data spencs( 12) / +.4003991310 8833118230 7258044590 8 d-10 /
data spencs( 13) / +.5823462892 0446384713 6813583575 7 d-11 /
data spencs( 14) / +.8576708692 6386892780 9791477122 4 d-12 /
data spencs( 15) / +.1276862586 2801930459 8948303343 3 d-12 /
data spencs( 16) / +.1918826209 0425170811 6238041606 2 d-13 /
data spencs( 17) / +.2907319206 9771381777 9579971967 3 d-14 /
data spencs( 18) / +.4437112685 2767804625 5747364174 5 d-15 /
data spencs( 19) / +.6815727787 4145995278 6735913560 7 d-16 /
data spencs( 20) / +.1053017386 0155744295 4701941664 4 d-16 /
data spencs( 21) / +.1635389806 7523771000 5182173457 0 d-17 /
data spencs( 22) / +.2551852874 9404639323 1090164258 1 d-18 /
data spencs( 23) / +.3999020621 9993601127 7047037951 9 d-19 /
data spencs( 24) / +.6291501645 2168118765 1414917119 9 d-20 /
data spencs( 25) / +.9933827435 6756776438 0388775253 3 d-21 /
data spencs( 26) / +.1573679570 7499648167 2176380586 6 d-21 /
data spencs( 27) / +.2500595316 8494761293 6927095466 6 d-22 /
data spencs( 28) / +.3984740918 3838111392 1066325333 3 d-23 /
data spencs( 29) / +.6366473210 0828438926 9132629333 3 d-24 /
data spencs( 30) / +.1019674287 2396783670 7706197333 3 d-24 /
data spencs( 31) / +.1636881058 9135188411 1107413333 3 d-25 /
data spencs( 32) / +.2633310439 4176501173 4527999999 9 d-26 /
data spencs( 33) / +.4244811560 1239768172 2436266666 6 d-27 /
data spencs( 34) / +.6855411983 6800529168 2474666666 6 d-28 /
data spencs( 35) / +.1109122433 4380564340 1898666666 6 d-28 /
data spencs( 36) / +.1797431304 9998914573 6533333333 3 d-29 /
data spencs( 37) / +.2917505845 9760951732 9066666666 6 d-30 /
data spencs( 38) / +.4742646808 9286710613 3333333333 3 d-31 /
c
c pi26 = pi**2/6.0
data pi26 / +1.644934066 8482264364 7241516664 6025189219 d0 /
data nspenc, xbig / 0, 0.0d0 /
c
if (nspenc.ne.0) go to 10
nspenc = initds (spencs, 38, 0.1*sngl(d1mach(3)))
xbig = 1.0d0/d1mach(3)
c
10 if (x.gt.2.0d0) go to 60
if (x.gt.1.0d0) go to 50
if (x.gt.0.5d0) go to 40
if (x.ge.0.0d0) go to 30
if (x.gt.(-1.d0)) go to 20
c
c here if x .le. -1.0
c
aln = dlog(1.0d0-x)
dspenc = -pi26 - 0.5d0*aln*(2.0d0*dlog(-x)-aln)
if (x.gt.(-xbig)) dspenc = dspenc
1 + (1.d0 + dcsevl (4.d0/(1.d0-x)-1.d0, spencs, nspenc))/(1.d0-x)
return
c
c -1.0 .lt. x .lt. 0.0
c
20 dspenc = -0.5d0*dlog(1.0d0-x)**2
1 - x*(1.d0+dcsevl(4.d0*x/(x-1.d0)-1.d0, spencs, nspenc))/(x-1.d0)
return
c
c 0.0 .le. x .le. 0.5
c
30 dspenc = x*(1.d0 + dcsevl (4.d0*x-1.d0, spencs, nspenc))
return
c
c 0.5 .lt. x .le. 1.0
c
40 dspenc = pi26
if (x.ne.1.d0) dspenc = pi26 - dlog(x)*dlog(1.0d0-x)
1 - (1.d0-x)*(1.d0+dcsevl(4.d0*(1.d0-x)-1.d0, spencs, nspenc))
return
c
c 1.0 .lt. x .le. 2.0
c
50 dspenc = pi26 - 0.5d0*dlog(x)*dlog((x-1.d0)**2/x)
1 + (x-1.d0)*(1.d0+dcsevl(4.d0*(x-1.d0)/x-1.d0, spencs, nspenc))/x
return
c
c x .gt. 2.0
c
60 dspenc = 2.0d0*pi26 - 0.5d0*dlog(x)**2
if (x.lt.xbig) dspenc = dspenc
1 - (1.d0 + dcsevl (4.d0/x-1.d0, spencs, nspenc))/x
return
c
end
|