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