Blame view

fvn_fnlib/besk0e.f 4.04 KB
38581db0c   daniau   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
        function besk0e (x)
  c july 1980 version.  w. fullerton, c3, los alamos scientific lab.
        dimension bk0cs(11), ak0cs(17), ak02cs(14)
        external  besi0, csevl, inits, r1mach
  c
  c series for bk0        on the interval  0.          to  4.00000d+00
  c                                        with weighted error   3.57e-19
  c                                         log weighted error  18.45
  c                               significant figures required  17.99
  c                                    decimal places required  18.97
  c
        data bk0 cs( 1) /   -.0353273932 3390276872e0 /
        data bk0 cs( 2) /    .3442898999 246284869e0 /
        data bk0 cs( 3) /    .0359799365 1536150163e0 /
        data bk0 cs( 4) /    .0012646154 1144692592e0 /
        data bk0 cs( 5) /    .0000228621 2103119451e0 /
        data bk0 cs( 6) /    .0000002534 7910790261e0 /
        data bk0 cs( 7) /    .0000000019 0451637722e0 /
        data bk0 cs( 8) /    .0000000000 1034969525e0 /
        data bk0 cs( 9) /    .0000000000 0004259816e0 /
        data bk0 cs(10) /    .0000000000 0000013744e0 /
        data bk0 cs(11) /    .0000000000 0000000035e0 /
  c
  c series for ak0        on the interval  1.25000d-01 to  5.00000d-01
  c                                        with weighted error   5.34e-17
  c                                         log weighted error  16.27
  c                               significant figures required  14.92
  c                                    decimal places required  16.89
  c
        data ak0 cs( 1) /   -.0764394790 3327941e0 /
        data ak0 cs( 2) /   -.0223565260 5699819e0 /
        data ak0 cs( 3) /    .0007734181 1546938e0 /
        data ak0 cs( 4) /   -.0000428100 6688886e0 /
        data ak0 cs( 5) /    .0000030817 0017386e0 /
        data ak0 cs( 6) /   -.0000002639 3672220e0 /
        data ak0 cs( 7) /    .0000000256 3713036e0 /
        data ak0 cs( 8) /   -.0000000027 4270554e0 /
        data ak0 cs( 9) /    .0000000003 1694296e0 /
        data ak0 cs(10) /   -.0000000000 3902353e0 /
        data ak0 cs(11) /    .0000000000 0506804e0 /
        data ak0 cs(12) /   -.0000000000 0068895e0 /
        data ak0 cs(13) /    .0000000000 0009744e0 /
        data ak0 cs(14) /   -.0000000000 0001427e0 /
        data ak0 cs(15) /    .0000000000 0000215e0 /
        data ak0 cs(16) /   -.0000000000 0000033e0 /
        data ak0 cs(17) /    .0000000000 0000005e0 /
  c
  c series for ak02       on the interval  0.          to  1.25000d-01
  c                                        with weighted error   2.34e-17
  c                                         log weighted error  16.63
  c                               significant figures required  14.67
  c                                    decimal places required  17.20
  c
        data ak02cs( 1) /   -.0120186982 6307592e0 /
        data ak02cs( 2) /   -.0091748526 9102569e0 /
        data ak02cs( 3) /    .0001444550 9317750e0 /
        data ak02cs( 4) /   -.0000040136 1417543e0 /
        data ak02cs( 5) /    .0000001567 8318108e0 /
        data ak02cs( 6) /   -.0000000077 7011043e0 /
        data ak02cs( 7) /    .0000000004 6111825e0 /
        data ak02cs( 8) /   -.0000000000 3158592e0 /
        data ak02cs( 9) /    .0000000000 0243501e0 /
        data ak02cs(10) /   -.0000000000 0020743e0 /
        data ak02cs(11) /    .0000000000 0001925e0 /
        data ak02cs(12) /   -.0000000000 0000192e0 /
        data ak02cs(13) /    .0000000000 0000020e0 /
        data ak02cs(14) /   -.0000000000 0000002e0 /
  c
        data ntk0, ntak0, ntak02, xsml / 3*0, 0. /
  c
        if (ntk0.ne.0) go to 10
        ntk0 = inits (bk0cs, 11, 0.1*r1mach(3))
        ntak0 = inits (ak0cs, 17, 0.1*r1mach(3))
        ntak02 = inits (ak02cs, 14, 0.1*r1mach(3))
        xsml = sqrt (4.0*r1mach(3))
  c
   10   if (x.le.0.) call seteru (29hbesk0e  x is zero or negative, 29,
       1  1, 2)
        if (x.gt.2.) go to 20
  c
        y = 0.
        if (x.gt.xsml) y = x*x
        besk0e = exp(x) * (-alog(0.5*x)*besi0(x)
       1  - .25 + csevl (.5*y-1., bk0cs, ntk0) )
        return
  c
   20   if (x.le.8.) besk0e = (1.25 + csevl ((16./x-5.)/3., ak0cs, ntak0))
       1  / sqrt(x)
        if (x.gt.8.) besk0e = (1.25 + csevl (16./x-1., ak02cs, ntak02))
       1  / sqrt(x)
  c
        return
        end