dbesy0.f
2.62 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
double precision function dbesy0 (x)
c august 1980 edition. w. fullerton, c3, los alamos scientific lab.
double precision x, by0cs(19), ampl, theta, twodpi, xsml,
1 y, alnhaf, d1mach, dcsevl, dbesj0, dlog, dsin, dsqrt
external d1mach, dbesj0, dcsevl, initds
c
c series for by0 on the interval 0. to 1.60000e+01
c with weighted error 8.14e-32
c log weighted error 31.09
c significant figures required 30.31
c decimal places required 31.73
c
data by0 cs( 1) / -.1127783939 2865573217 9398054602 8 d-1 /
data by0 cs( 2) / -.1283452375 6042034604 8088453183 8 d+0 /
data by0 cs( 3) / -.1043788479 9794249365 8176227661 8 d+0 /
data by0 cs( 4) / +.2366274918 3969695409 2415926461 3 d-1 /
data by0 cs( 5) / -.2090391647 7004862391 9622395034 2 d-2 /
data by0 cs( 6) / +.1039754539 3905725209 9924657638 1 d-3 /
data by0 cs( 7) / -.3369747162 4239720967 1877534503 7 d-5 /
data by0 cs( 8) / +.7729384267 6706671585 2136721637 1 d-7 /
data by0 cs( 9) / -.1324976772 6642595914 4347606896 4 d-8 /
data by0 cs( 10) / +.1764823261 5404527921 0038936315 8 d-10 /
data by0 cs( 11) / -.1881055071 5801962006 0282301206 9 d-12 /
data by0 cs( 12) / +.1641865485 3661495027 9223718574 9 d-14 /
data by0 cs( 13) / -.1195659438 6046060857 4599100672 0 d-16 /
data by0 cs( 14) / +.7377296297 4401858424 9411242666 6 d-19 /
data by0 cs( 15) / -.3906843476 7104373307 4090666666 6 d-21 /
data by0 cs( 16) / +.1795503664 4361579498 2912000000 0 d-23 /
data by0 cs( 17) / -.7229627125 4480104789 3333333333 3 d-26 /
data by0 cs( 18) / +.2571727931 6351685973 3333333333 3 d-28 /
data by0 cs( 19) / -.8141268814 1636949333 3333333333 3 d-31 /
c
data twodpi / 0.6366197723 6758134307 5535053490 057 d0 /
data alnhaf /-0.6931471805 5994530941 7232121458 18d0 /
data nty0, xsml / 0, 0.d0 /
c
if (nty0.ne.0) go to 10
nty0 = initds (by0cs, 19, 0.1*sngl(d1mach(3)))
xsml = dsqrt (4.0d0*d1mach(3))
c
10 if (x.le.0.d0) call seteru (29hdbesy0 x is zero or negative, 29,
1 1, 2)
if (x.gt.4.0d0) go to 20
c
y = 0.d0
if (x.gt.xsml) y = x*x
dbesy0 = twodpi*(alnhaf+dlog(x))*dbesj0(x) + .375d0
1 + dcsevl (.125d0*y-1.d0, by0cs, nty0)
return
c
20 call d9b0mp (x, ampl, theta)
dbesy0 = ampl * dsin(theta)
return
c
end