dbesj0.f
2.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
double precision function dbesj0 (x)
c july 1977 edition. w. fullerton, c3, los alamos scientific lab.
double precision x, bj0cs(19), ampl, theta, xsml, y,
1 d1mach, dcsevl, dcos, dsqrt
external d1mach, dcsevl, initds
c
c series for bj0 on the interval 0. to 1.60000e+01
c with weighted error 4.39e-32
c log weighted error 31.36
c significant figures required 31.21
c decimal places required 32.00
c
data bj0 cs( 1) / +.1002541619 6893913701 0731272640 74 d+0 /
data bj0 cs( 2) / -.6652230077 6440513177 6787578311 24 d+0 /
data bj0 cs( 3) / +.2489837034 9828131370 4604687266 80 d+0 /
data bj0 cs( 4) / -.3325272317 0035769653 8843415038 54 d-1 /
data bj0 cs( 5) / +.2311417930 4694015462 9049241177 29 d-2 /
data bj0 cs( 6) / -.9911277419 9508092339 0485193365 49 d-4 /
data bj0 cs( 7) / +.2891670864 3998808884 7339037470 78 d-5 /
data bj0 cs( 8) / -.6121085866 3032635057 8184074815 16 d-7 /
data bj0 cs( 9) / +.9838650793 8567841324 7687486364 15 d-9 /
data bj0 cs( 10) / -.1242355159 7301765145 5158970068 36 d-10 /
data bj0 cs( 11) / +.1265433630 2559045797 9158272103 63 d-12 /
data bj0 cs( 12) / -.1061945649 5287244546 9148175129 59 d-14 /
data bj0 cs( 13) / +.7470621075 8024567437 0989155840 00 d-17 /
data bj0 cs( 14) / -.4469703227 4412780547 6270079999 99 d-19 /
data bj0 cs( 15) / +.2302428158 4337436200 5230933333 33 d-21 /
data bj0 cs( 16) / -.1031914479 4166698148 5226666666 66 d-23 /
data bj0 cs( 17) / +.4060817827 4873322700 8000000000 00 d-26 /
data bj0 cs( 18) / -.1414383600 5240913919 9999999999 99 d-28 /
data bj0 cs( 19) / +.4391090549 6698880000 0000000000 00 d-31 /
c
data ntj0, xsml / 0, 0.d0 /
c
if (ntj0.ne.0) go to 10
ntj0 = initds (bj0cs, 19, 0.1*sngl(d1mach(3)))
xsml = dsqrt (4.0d0*d1mach(3))
c
10 y = dabs(x)
if (y.gt.4.0d0) go to 20
c
dbesj0 = 1.0d0
if (y.gt.xsml) dbesj0 = dcsevl (.125d0*y*y-1.d0, bj0cs, ntj0)
return
c
20 call d9b0mp (y, ampl, theta)
dbesj0 = ampl * dcos(theta)
c
return
end