spenc.f
3.33 KB
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
function spenc (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 - spenc (a*(c*z+d)/(a*d-b*c)) / c.
c
c ref -- k. mitchell, philosophical magazine, 40, p.351 (1949).
c stegun and abromowitz, ams 55, p.1004.
c
dimension spencs(19)
external csevl, inits, r1mach
c
c series for spen on the interval 0. to 5.00000d-01
c with weighted error 6.82e-17
c log weighted error 16.17
c significant figures required 15.22
c decimal places required 16.81
c
data spencs( 1) / .1527365598 892406e0 /
data spencs( 2) / .0816965805 8051014e0 /
data spencs( 3) / .0058141571 4077873e0 /
data spencs( 4) / .0005371619 8145415e0 /
data spencs( 5) / .0000572470 4675185e0 /
data spencs( 6) / .0000066745 4612164e0 /
data spencs( 7) / .0000008276 4673397e0 /
data spencs( 8) / .0000001073 3156730e0 /
data spencs( 9) / .0000000144 0077294e0 /
data spencs(10) / .0000000019 8444202e0 /
data spencs(11) / .0000000002 7940058e0 /
data spencs(12) / .0000000000 4003991e0 /
data spencs(13) / .0000000000 0582346e0 /
data spencs(14) / .0000000000 0085767e0 /
data spencs(15) / .0000000000 0012768e0 /
data spencs(16) / .0000000000 0001918e0 /
data spencs(17) / .0000000000 0000290e0 /
data spencs(18) / .0000000000 0000044e0 /
data spencs(19) / .0000000000 0000006e0 /
c
c pi26 = pi**2/6.0
data pi26 / 1.644934066 848226e0 /
data nspenc, xbig / 0, 0.0 /
c
if (nspenc.ne.0) go to 10
nspenc = inits (spencs, 19, 0.1*r1mach(3))
xbig = 1.0/r1mach(3)
c
10 if (x.gt.2.0) go to 60
if (x.gt.1.0) go to 50
if (x.gt.0.5) go to 40
if (x.ge.0.0) go to 30
if (x.gt.(-1.)) go to 20
c
c here if x .le. -1.0
c
aln = alog(1.0-x)
spenc = -pi26 - 0.5*aln*(2.0*alog(-x)-aln)
if (x.gt.(-xbig)) spenc = spenc
1 + (1.0 + csevl (4.0/(1.0-x)-1.0, spencs, nspenc)) / (1.0-x)
return
c
c -1.0 .lt. x .lt. 0.0
c
20 spenc = -0.5*alog(1.0-x)**2
1 - x*(1.0 + csevl (4.0*x/(x-1.0)-1.0, spencs, nspenc)) / (x-1.0)
return
c
c 0.0 .le. x .le. 0.5
c
30 spenc = x*(1.0 + csevl (4.0*x-1.0, spencs, nspenc))
return
c
c 0.5 .lt. x .le. 1.0
c
40 spenc = pi26
if (x.ne.1.0) spenc = pi26 - alog(x)*alog(1.0-x)
1 - (1.0-x)*(1.0 + csevl (4.0*(1.0-x)-1.0, spencs, nspenc))
return
c
c 1.0 .lt. x .le. 2.0
c
50 spenc = pi26 - 0.5*alog(x)*alog((x-1.0)**2/x)
1 + (x-1.)*(1.0 + csevl (4.0*(x-1.)/x-1.0, spencs, nspenc))/x
return
c
c x .gt. 2.0
c
60 spenc = 2.0*pi26 - 0.5*alog(x)**2
if (x.lt.xbig) spenc = spenc
1 - (1.0 + csevl (4.0/x-1.0, spencs, nspenc))/x
return
c
end