Blame view

fvn_fnlib/r9sifg.f 8.16 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
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
        subroutine r9sifg (x, f, g)
  c december 1980 edition.  w. fullerton, bell labs
        dimension f1cs(20), f2cs(29), g1cs(21), g2cs(34)
        external csevl, inits, r1mach
  c
  c series for f1   on the interval  2.00000e-02 to  6.25000e-02
  c                                        with weighted error   2.82e-17
  c                                         log weighted error  16.55
  c                               significant figures required  15.36
  c                                    decimal places required  17.20
  c
        data f1  cs(  1) /    -0.1191081969 051363610 e0/
        data f1  cs(  2) /    -0.0247823144 996236248 e0/
        data f1  cs(  3) /     0.0011910281 453357821 e0/
        data f1  cs(  4) /    -0.0000927027 714388562 e0/
        data f1  cs(  5) /     0.0000093373 141568271 e0/
        data f1  cs(  6) /    -0.0000011058 287820557 e0/
        data f1  cs(  7) /     0.0000001464 772071460 e0/
        data f1  cs(  8) /    -0.0000000210 694496288 e0/
        data f1  cs(  9) /     0.0000000032 293492367 e0/
        data f1  cs( 10) /    -0.0000000005 206529618 e0/
        data f1  cs( 11) /     0.0000000000 874878885 e0/
        data f1  cs( 12) /    -0.0000000000 152176187 e0/
        data f1  cs( 13) /     0.0000000000 027257192 e0/
        data f1  cs( 14) /    -0.0000000000 005007053 e0/
        data f1  cs( 15) /     0.0000000000 000940241 e0/
        data f1  cs( 16) /    -0.0000000000 000180014 e0/
        data f1  cs( 17) /     0.0000000000 000035063 e0/
        data f1  cs( 18) /    -0.0000000000 000006935 e0/
        data f1  cs( 19) /     0.0000000000 000001391 e0/
        data f1  cs( 20) /    -0.0000000000 000000282 e0/
  c
  c series for f2   on the interval  0.00000e+00 to  2.00000e-02
  c                                        with weighted error   4.32e-17
  c                                         log weighted error  16.36
  c                               significant figures required  14.75
  c                                    decimal places required  17.10
  c
        data f2  cs(  1) /    -0.0348409253 897013234 e0/
        data f2  cs(  2) /    -0.0166842205 677959686 e0/
        data f2  cs(  3) /     0.0006752901 241237738 e0/
        data f2  cs(  4) /    -0.0000535066 622544701 e0/
        data f2  cs(  5) /     0.0000062693 421779007 e0/
        data f2  cs(  6) /    -0.0000009526 638801991 e0/
        data f2  cs(  7) /     0.0000001745 629224251 e0/
        data f2  cs(  8) /    -0.0000000368 795403065 e0/
        data f2  cs(  9) /     0.0000000087 202677705 e0/
        data f2  cs( 10) /    -0.0000000022 601970392 e0/
        data f2  cs( 11) /     0.0000000006 324624977 e0/
        data f2  cs( 12) /    -0.0000000001 888911889 e0/
        data f2  cs( 13) /     0.0000000000 596774674 e0/
        data f2  cs( 14) /    -0.0000000000 198044313 e0/
        data f2  cs( 15) /     0.0000000000 068641396 e0/
        data f2  cs( 16) /    -0.0000000000 024731020 e0/
        data f2  cs( 17) /     0.0000000000 009226360 e0/
        data f2  cs( 18) /    -0.0000000000 003552364 e0/
        data f2  cs( 19) /     0.0000000000 001407606 e0/
        data f2  cs( 20) /    -0.0000000000 000572623 e0/
        data f2  cs( 21) /     0.0000000000 000238654 e0/
        data f2  cs( 22) /    -0.0000000000 000101714 e0/
        data f2  cs( 23) /     0.0000000000 000044259 e0/
        data f2  cs( 24) /    -0.0000000000 000019634 e0/
        data f2  cs( 25) /     0.0000000000 000008868 e0/
        data f2  cs( 26) /    -0.0000000000 000004074 e0/
        data f2  cs( 27) /     0.0000000000 000001901 e0/
        data f2  cs( 28) /    -0.0000000000 000000900 e0/
        data f2  cs( 29) /     0.0000000000 000000432 e0/
  c
  c series for g1   on the interval  2.00000e-02 to  6.25000e-02
  c                                        with weighted error   5.48e-17
  c                                         log weighted error  16.26
  c                               significant figures required  15.47
  c                                    decimal places required  16.92
  c
        data g1  cs(  1) /    -0.3040578798 253495954 e0/
        data g1  cs(  2) /    -0.0566890984 597120588 e0/
        data g1  cs(  3) /     0.0039046158 173275644 e0/
        data g1  cs(  4) /    -0.0003746075 959202261 e0/
        data g1  cs(  5) /     0.0000435431 556559844 e0/
        data g1  cs(  6) /    -0.0000057417 294453025 e0/
        data g1  cs(  7) /     0.0000008282 552104503 e0/
        data g1  cs(  8) /    -0.0000001278 245892595 e0/
        data g1  cs(  9) /     0.0000000207 978352949 e0/
        data g1  cs( 10) /    -0.0000000035 313205922 e0/
        data g1  cs( 11) /     0.0000000006 210824236 e0/
        data g1  cs( 12) /    -0.0000000001 125215474 e0/
        data g1  cs( 13) /     0.0000000000 209088918 e0/
        data g1  cs( 14) /    -0.0000000000 039715832 e0/
        data g1  cs( 15) /     0.0000000000 007690431 e0/
        data g1  cs( 16) /    -0.0000000000 001514697 e0/
        data g1  cs( 17) /     0.0000000000 000302892 e0/
        data g1  cs( 18) /    -0.0000000000 000061400 e0/
        data g1  cs( 19) /     0.0000000000 000012601 e0/
        data g1  cs( 20) /    -0.0000000000 000002615 e0/
        data g1  cs( 21) /     0.0000000000 000000548 e0/
  c
  c series for g2   on the interval  0.00000e+00 to  2.00000e-02
  c                                        with weighted error   5.01e-17
  c                                         log weighted error  16.30
  c                               significant figures required  15.12
  c                                    decimal places required  17.07
  c
        data g2  cs(  1) /    -0.0967329367 532432218 e0/
        data g2  cs(  2) /    -0.0452077907 957459871 e0/
        data g2  cs(  3) /     0.0028190005 352706523 e0/
        data g2  cs(  4) /    -0.0002899167 740759160 e0/
        data g2  cs(  5) /     0.0000407444 664601121 e0/
        data g2  cs(  6) /    -0.0000071056 382192354 e0/
        data g2  cs(  7) /     0.0000014534 723163019 e0/
        data g2  cs(  8) /    -0.0000003364 116512503 e0/
        data g2  cs(  9) /     0.0000000859 774367886 e0/
        data g2  cs( 10) /    -0.0000000238 437656302 e0/
        data g2  cs( 11) /     0.0000000070 831906340 e0/
        data g2  cs( 12) /    -0.0000000022 318068154 e0/
        data g2  cs( 13) /     0.0000000007 401087359 e0/
        data g2  cs( 14) /    -0.0000000002 567171162 e0/
        data g2  cs( 15) /     0.0000000000 926707021 e0/
        data g2  cs( 16) /    -0.0000000000 346693311 e0/
        data g2  cs( 17) /     0.0000000000 133950573 e0/
        data g2  cs( 18) /    -0.0000000000 053290754 e0/
        data g2  cs( 19) /     0.0000000000 021775312 e0/
        data g2  cs( 20) /    -0.0000000000 009118621 e0/
        data g2  cs( 21) /     0.0000000000 003905864 e0/
        data g2  cs( 22) /    -0.0000000000 001708459 e0/
        data g2  cs( 23) /     0.0000000000 000762015 e0/
        data g2  cs( 24) /    -0.0000000000 000346151 e0/
        data g2  cs( 25) /     0.0000000000 000159996 e0/
        data g2  cs( 26) /    -0.0000000000 000075213 e0/
        data g2  cs( 27) /     0.0000000000 000035970 e0/
        data g2  cs( 28) /    -0.0000000000 000017530 e0/
        data g2  cs( 29) /     0.0000000000 000008738 e0/
        data g2  cs( 30) /    -0.0000000000 000004487 e0/
        data g2  cs( 31) /     0.0000000000 000002397 e0/
        data g2  cs( 32) /    -0.0000000000 000001347 e0/
        data g2  cs( 33) /     0.0000000000 000000801 e0/
        data g2  cs( 34) /    -0.0000000000 000000501 e0/
  c
        data nf1, nf2, ng1, ng2 / 4*0 /
        data xbnd, xbig, xmaxf, xmaxg / 4*0.0 /
  c
        if (nf1.ne.0) go to 10
        tol = 0.1*r1mach(3)
        nf1 = inits (f1cs, 20, tol)
        nf2 = inits (f2cs, 29, tol)
        ng1 = inits (g1cs, 21, tol)
        ng2 = inits (g2cs, 34, tol)
  c
        xbig = sqrt (1.0/r1mach(3))
        xmaxf = exp (amin1(-alog(r1mach(1)), alog(r1mach(2))) - 0.01)
        xmaxg = 1.0/sqrt(r1mach(1))
        xbnd = sqrt(50.0)
  c
   10   if (x.lt.4.0) call seteru (
       1  34hr9sifg  approxs invalid for x lt 4, 34, 1, 2)
  c
        if (x.gt.xbnd) go to 20
        f = (1.0 + csevl ((1.0/x**2-0.04125)/0.02125, f1cs, nf1))/x
        g = (1.0 + csevl ((1.0/x**2-0.04125)/0.02125, g1cs, ng1))/x**2
        return
  c
   20   if (x.gt.xbig) go to 30
        f = (1.0 + csevl (100./x**2-1., f2cs, nf2))/x
        g = (1.0 + csevl (100./x**2-1., g2cs, ng2))/x**2
        return
  c
   30   f = 0.0
        if (x.lt.xmaxf) f = 1.0/x
        g = 0.0
        if (x.lt.xmaxg) g = 1.0/x**2
        return
  c
        end