Blame view

fvn_fnlib/dbide.f 17.3 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
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
        double precision function dbide (x)
  c july 1980 edition.  w. fullerton, bell labs.
  c
  c evaluate the derivative of the airy function bi(x) for x .le. 0
  c and evaluate  bid(x) * exp(-zeta)  for x .ge. 0  where
  c      zeta = 2/3 * x**(3/2)
  c
        double precision x, bifcs(13), bigcs(13), bif2cs(15), big2cs(16),
       1  bip1cs(47), bip2cs(88), atr, btr, x2sml, x3sml, xbig, xn, phi,
       2  x3, x32sml, x2, z, sqrtx,   d1mach, dcsevl, dexp, dsin, dsqrt
        external d1mach, dcsevl,  initds
  c
  c series for bif    on the interval -1.00000e+00 to  1.00000e+00
  c                                        with weighted error   3.82e-34
  c                                         log weighted error  33.42
  c                               significant figures required  32.21
  c                                    decimal places required  33.98
  c
        data bif cs(  1) /  0.1153536790 8285702426 7474446284 908879d0/
        data bif cs(  2) /  0.0205007894 0491928753 0357789445 940252d0/
        data bif cs(  3) /  0.0002135290 2789028758 1892679619 451158d0/
        data bif cs(  4) /  0.0000010783 9606146768 3042209155 523569d0/
        data bif cs(  5) /  0.0000000032 0947088332 0666783353 670420d0/
        data bif cs(  6) /  0.0000000000 0629304067 1833540390 213316d0/
        data bif cs(  7) /  0.0000000000 0000874030 4300063083 340121d0/
        data bif cs(  8) /  0.0000000000 0000000904 7915683496 049529d0/
        data bif cs(  9) /  0.0000000000 0000000000 7249923164 709251d0/
        data bif cs( 10) /  0.0000000000 0000000000 0004629576 649604d0/
        data bif cs( 11) /  0.0000000000 0000000000 0000002411 236436d0/
        data bif cs( 12) /  0.0000000000 0000000000 0000000001 043825d0/
        data bif cs( 13) /  0.0000000000 0000000000 0000000000 000382d0/
  c
  c series for big    on the interval -1.00000e+00 to  1.00000e+00
  c                                        with weighted error   4.79e-32
  c                                         log weighted error  31.32
  c                               significant figures required  30.52
  c                                    decimal places required  31.88
  c
        data big cs(  1) / -0.0971964404 1644353738 9779097460 6802d0/
        data big cs(  2) /  0.1495035768 4316706657 1084344532 6264d0/
        data big cs(  3) /  0.0031135253 8712132604 1941917683 9631d0/
        data big cs(  4) /  0.0000247085 7057982129 6777702192 0569d0/
        data big cs(  5) /  0.0000001029 4962773137 8608198732 4295d0/
        data big cs(  6) /  0.0000000002 6397037398 6943289267 6778d0/
        data big cs(  7) /  0.0000000000 0045827927 0780320660 8181d0/
        data big cs(  8) /  0.0000000000 0000057428 2974089344 7321d0/
        data big cs(  9) /  0.0000000000 0000000054 3827538523 8549d0/
        data big cs( 10) /  0.0000000000 0000000000 0402834726 7083d0/
        data big cs( 11) /  0.0000000000 0000000000 0000239782 3826d0/
        data big cs( 12) /  0.0000000000 0000000000 0000000117 1956d0/
        data big cs( 13) /  0.0000000000 0000000000 0000000000 0479d0/
  c
  c series for bif2 on the interval  1.00000e+00 to  8.00000e+00
  c                                        with weighted error   1.38e-33
  c                                         log weighted error  32.86
  c                               significant figures required  32.12
  c                                    decimal places required  33.45
  c
        data bif2cs(  1) /  0.3234939876 0352203352 1191935962 66015d0/
        data bif2cs(  2) /  0.0862978715 3556355913 8888353238 11100d0/
        data bif2cs(  3) /  0.0029940255 5265539742 6138210507 27155d0/
        data bif2cs(  4) /  0.0000514305 2836466163 7204643169 50821d0/
        data bif2cs(  5) /  0.0000005258 4025003681 1460260330 98613d0/
        data bif2cs(  6) /  0.0000000035 6175137395 7700281027 30600d0/
        data bif2cs(  7) /  0.0000000000 1714686400 7145848305 18308d0/
        data bif2cs(  8) /  0.0000000000 0006166351 9692325554 06693d0/
        data bif2cs(  9) /  0.0000000000 0000017191 0821543159 85806d0/
        data bif2cs( 10) /  0.0000000000 0000000038 2368895188 03943d0/
        data bif2cs( 11) /  0.0000000000 0000000000 0694241736 24884d0/
        data bif2cs( 12) /  0.0000000000 0000000000 0001048339 32510d0/
        data bif2cs( 13) /  0.0000000000 0000000000 0000001337 21972d0/
        data bif2cs( 14) /  0.0000000000 0000000000 0000000001 45986d0/
        data bif2cs( 15) /  0.0000000000 0000000000 0000000000 00138d0/
  c
  c series for big2 on the interval  1.00000e+00 to  8.00000e+00
  c                                        with weighted error   1.91e-34
  c                                         log weighted error  33.72
  c                               significant figures required  33.76
  c                                    decimal places required  34.32
  c
        data big2cs(  1) /  1.6062999463 6212945775 9284537862 622883d0/
        data big2cs(  2) /  0.7449088819 8760886520 1476685194 753972d0/
        data big2cs(  3) /  0.0470138738 6102773796 4095177635 353019d0/
        data big2cs(  4) /  0.0012284422 0625482390 7016188785 848091d0/
        data big2cs(  5) /  0.0000173222 4122566236 2670987355 613727d0/
        data big2cs(  6) /  0.0000001521 9016523680 1893711508 366563d0/
        data big2cs(  7) /  0.0000000009 1135602491 1957704145 528786d0/
        data big2cs(  8) /  0.0000000000 0395479184 2356644201 722554d0/
        data big2cs(  9) /  0.0000000000 0001300173 7033862320 007309d0/
        data big2cs( 10) /  0.0000000000 0000003349 3506858269 079763d0/
        data big2cs( 11) /  0.0000000000 0000000006 9419094403 694057d0/
        data big2cs( 12) /  0.0000000000 0000000000 0118248256 604581d0/
        data big2cs( 13) /  0.0000000000 0000000000 0000168462 493472d0/
        data big2cs( 14) /  0.0000000000 0000000000 0000000203 684674d0/
        data big2cs( 15) /  0.0000000000 0000000000 0000000000 211619d0/
        data big2cs( 16) /  0.0000000000 0000000000 0000000000 000191d0/
  c
  c series for bip2 on the interval  0.00000e+00 to  1.25000e-01
  c                                        with weighted error   4.37e-33
  c                                         log weighted error  32.36
  c                               significant figures required  31.18
  c                                    decimal places required  33.33
  c
        data bip2cs(  1) / -0.1326970544 3526630494 9370312102 17135d0/
        data bip2cs(  2) / -0.0056844362 6045977481 3060463390 37428d0/
        data bip2cs(  3) / -0.0001564360 1119611609 6236984712 16660d0/
        data bip2cs(  4) / -0.0000113673 7203679562 2673360532 07940d0/
        data bip2cs(  5) / -0.0000014346 4350991283 6696431369 51338d0/
        data bip2cs(  6) / -0.0000001809 8531185164 1318687464 81700d0/
        data bip2cs(  7) /  0.0000000092 6177343610 8655462295 11422d0/
        data bip2cs(  8) /  0.0000000171 0005490720 5921818872 96162d0/
        data bip2cs(  9) /  0.0000000047 6698163503 7817082526 86849d0/
        data bip2cs( 10) / -0.0000000003 5195022023 1631419453 97159d0/
        data bip2cs( 11) / -0.0000000005 8890614315 8868715741 47635d0/
        data bip2cs( 12) / -0.0000000000 6678499607 7955375976 12089d0/
        data bip2cs( 13) /  0.0000000000 6395565101 7203911906 97713d0/
        data bip2cs( 14) /  0.0000000000 1554529427 0643941064 03245d0/
        data bip2cs( 15) / -0.0000000000 0792396999 7446129716 84001d0/
        data bip2cs( 16) / -0.0000000000 0258326242 6897177989 47525d0/
        data bip2cs( 17) /  0.0000000000 0121655047 7878491179 78773d0/
        data bip2cs( 18) /  0.0000000000 0038707207 1728999859 42258d0/
        data bip2cs( 19) / -0.0000000000 0022487045 4796182291 30656d0/
        data bip2cs( 20) / -0.0000000000 0004953476 5156840462 93493d0/
        data bip2cs( 21) /  0.0000000000 0004563781 6015269127 56017d0/
        data bip2cs( 22) /  0.0000000000 0000332998 3143450141 18494d0/
        data bip2cs( 23) / -0.0000000000 0000921750 1858328742 02719d0/
        data bip2cs( 24) /  0.0000000000 0000094156 6706589582 05765d0/
        data bip2cs( 25) /  0.0000000000 0000167153 9526407161 57721d0/
        data bip2cs( 26) / -0.0000000000 0000055134 2687821824 10852d0/
        data bip2cs( 27) / -0.0000000000 0000022368 6515720066 17795d0/
        data bip2cs( 28) /  0.0000000000 0000017486 9489765200 89209d0/
        data bip2cs( 29) /  0.0000000000 0000000206 5186663523 29750d0/
        data bip2cs( 30) / -0.0000000000 0000003973 0600181307 12479d0/
        data bip2cs( 31) /  0.0000000000 0000001154 8369357248 92335d0/
        data bip2cs( 32) /  0.0000000000 0000000553 9060536782 76421d0/
        data bip2cs( 33) / -0.0000000000 0000000457 1744273964 78267d0/
        data bip2cs( 34) /  0.0000000000 0000000026 5671118582 84432d0/
        data bip2cs( 35) /  0.0000000000 0000000101 5991481541 67823d0/
        data bip2cs( 36) / -0.0000000000 0000000044 8212312721 96246d0/
        data bip2cs( 37) / -0.0000000000 0000000007 9591496616 17295d0/
        data bip2cs( 38) /  0.0000000000 0000000014 5836156161 65794d0/
        data bip2cs( 39) / -0.0000000000 0000000004 0151278930 61405d0/
        data bip2cs( 40) / -0.0000000000 0000000002 0791529637 43616d0/
        data bip2cs( 41) /  0.0000000000 0000000001 9726304496 34388d0/
        data bip2cs( 42) / -0.0000000000 0000000000 3360334040 01683d0/
        data bip2cs( 43) / -0.0000000000 0000000000 3765048326 85507d0/
        data bip2cs( 44) /  0.0000000000 0000000000 2699355088 25595d0/
        data bip2cs( 45) / -0.0000000000 0000000000 0269859460 69808d0/
        data bip2cs( 46) / -0.0000000000 0000000000 0617940117 88222d0/
        data bip2cs( 47) /  0.0000000000 0000000000 0387826933 11711d0/
        data bip2cs( 48) / -0.0000000000 0000000000 0024200940 05071d0/
        data bip2cs( 49) / -0.0000000000 0000000000 0098440510 58925d0/
        data bip2cs( 50) /  0.0000000000 0000000000 0059543533 58494d0/
        data bip2cs( 51) / -0.0000000000 0000000000 0003612744 46366d0/
        data bip2cs( 52) / -0.0000000000 0000000000 0015526345 78088d0/
        data bip2cs( 53) /  0.0000000000 0000000000 0009778193 80304d0/
        data bip2cs( 54) / -0.0000000000 0000000000 0000922394 47509d0/
        data bip2cs( 55) / -0.0000000000 0000000000 0002415459 03934d0/
        data bip2cs( 56) /  0.0000000000 0000000000 0001695586 52255d0/
        data bip2cs( 57) / -0.0000000000 0000000000 0000267624 08641d0/
        data bip2cs( 58) / -0.0000000000 0000000000 0000361881 16265d0/
        data bip2cs( 59) /  0.0000000000 0000000000 0000303724 04951d0/
        data bip2cs( 60) / -0.0000000000 0000000000 0000074228 76903d0/
        data bip2cs( 61) / -0.0000000000 0000000000 0000049306 78544d0/
        data bip2cs( 62) /  0.0000000000 0000000000 0000054687 90028d0/
        data bip2cs( 63) / -0.0000000000 0000000000 0000019203 15188d0/
        data bip2cs( 64) / -0.0000000000 0000000000 0000005163 35154d0/
        data bip2cs( 65) /  0.0000000000 0000000000 0000009577 23167d0/
        data bip2cs( 66) / -0.0000000000 0000000000 0000004636 59079d0/
        data bip2cs( 67) / -0.0000000000 0000000000 0000000045 09226d0/
        data bip2cs( 68) /  0.0000000000 0000000000 0000001556 17519d0/
        data bip2cs( 69) / -0.0000000000 0000000000 0000001041 56509d0/
        data bip2cs( 70) /  0.0000000000 0000000000 0000000195 65323d0/
        data bip2cs( 71) /  0.0000000000 0000000000 0000000213 35380d0/
        data bip2cs( 72) / -0.0000000000 0000000000 0000000214 61958d0/
        data bip2cs( 73) /  0.0000000000 0000000000 0000000078 75791d0/
        data bip2cs( 74) /  0.0000000000 0000000000 0000000017 13768d0/
        data bip2cs( 75) / -0.0000000000 0000000000 0000000039 17137d0/
        data bip2cs( 76) /  0.0000000000 0000000000 0000000022 33559d0/
        data bip2cs( 77) / -0.0000000000 0000000000 0000000002 69383d0/
        data bip2cs( 78) / -0.0000000000 0000000000 0000000005 77764d0/
        data bip2cs( 79) /  0.0000000000 0000000000 0000000005 19650d0/
        data bip2cs( 80) / -0.0000000000 0000000000 0000000001 83361d0/
        data bip2cs( 81) / -0.0000000000 0000000000 0000000000 45763d0/
        data bip2cs( 82) /  0.0000000000 0000000000 0000000000 99235d0/
        data bip2cs( 83) / -0.0000000000 0000000000 0000000000 58938d0/
        data bip2cs( 84) /  0.0000000000 0000000000 0000000000 09568d0/
        data bip2cs( 85) /  0.0000000000 0000000000 0000000000 13758d0/
        data bip2cs( 86) / -0.0000000000 0000000000 0000000000 14066d0/
        data bip2cs( 87) /  0.0000000000 0000000000 0000000000 05964d0/
        data bip2cs( 88) /  0.0000000000 0000000000 0000000000 00437d0/
  c
  c series for bip1 on the interval  1.25000e-01 to  3.53553e-01
  c                                        with weighted error   1.75e-32
  c                                         log weighted error  31.76
  c                               significant figures required  30.70
  c                                    decimal places required  32.59
  c
        data bip1cs(  1) / -0.1729187351 0795537186 1246798237 41003d0/
        data bip1cs(  2) / -0.0149358492 9846943639 4862310218 18675d0/
        data bip1cs(  3) / -0.0005471104 9516785663 9906586978 74460d0/
        data bip1cs(  4) /  0.0001537966 2929584083 4495737278 56666d0/
        data bip1cs(  5) /  0.0000154353 4761921794 1310289480 22869d0/
        data bip1cs(  6) / -0.0000065434 1138519060 1292260871 06765d0/
        data bip1cs(  7) /  0.0000003728 0824078787 0322321522 75240d0/
        data bip1cs(  8) /  0.0000002072 0783881887 4800808107 10514d0/
        data bip1cs(  9) / -0.0000000658 1733364696 1916894958 83922d0/
        data bip1cs( 10) /  0.0000000074 9267463539 2882129860 48985d0/
        data bip1cs( 11) /  0.0000000011 1013368840 7071476988 90101d0/
        data bip1cs( 12) / -0.0000000007 2651405529 1595123238 80794d0/
        data bip1cs( 13) /  0.0000000001 7827235598 4701539621 65668d0/
        data bip1cs( 14) / -0.0000000000 2173463524 8095062696 56807d0/
        data bip1cs( 15) / -0.0000000000 0203020349 6538825940 17049d0/
        data bip1cs( 16) /  0.0000000000 0193118272 2940775193 19859d0/
        data bip1cs( 17) / -0.0000000000 0060449525 0482902960 23117d0/
        data bip1cs( 18) /  0.0000000000 0012094496 2489336642 77802d0/
        data bip1cs( 19) / -0.0000000000 0001251088 3600744797 84619d0/
        data bip1cs( 20) / -0.0000000000 0000199173 8324248813 44036d0/
        data bip1cs( 21) /  0.0000000000 0000151540 8163428643 03038d0/
        data bip1cs( 22) / -0.0000000000 0000049768 9270598162 40250d0/
        data bip1cs( 23) /  0.0000000000 0000011545 9597318105 01403d0/
        data bip1cs( 24) / -0.0000000000 0000001863 2868629079 83871d0/
        data bip1cs( 25) /  0.0000000000 0000000099 3303923447 59104d0/
        data bip1cs( 26) /  0.0000000000 0000000068 1820836674 12417d0/
        data bip1cs( 27) / -0.0000000000 0000000034 8544564796 50551d0/
        data bip1cs( 28) /  0.0000000000 0000000010 8603821342 35961d0/
        data bip1cs( 29) / -0.0000000000 0000000002 5992901852 40166d0/
        data bip1cs( 30) /  0.0000000000 0000000000 4768953704 59000d0/
        data bip1cs( 31) / -0.0000000000 0000000000 0519469407 77177d0/
        data bip1cs( 32) / -0.0000000000 0000000000 0059255750 44912d0/
        data bip1cs( 33) /  0.0000000000 0000000000 0057460089 70972d0/
        data bip1cs( 34) / -0.0000000000 0000000000 0021861198 06494d0/
        data bip1cs( 35) /  0.0000000000 0000000000 0006241242 94738d0/
        data bip1cs( 36) / -0.0000000000 0000000000 0001460034 21785d0/
        data bip1cs( 37) /  0.0000000000 0000000000 0000274938 93904d0/
        data bip1cs( 38) / -0.0000000000 0000000000 0000034746 78018d0/
        data bip1cs( 39) / -0.0000000000 0000000000 0000001093 03694d0/
        data bip1cs( 40) /  0.0000000000 0000000000 0000002619 72744d0/
        data bip1cs( 41) / -0.0000000000 0000000000 0000001123 65018d0/
        data bip1cs( 42) /  0.0000000000 0000000000 0000000351 52059d0/
        data bip1cs( 43) / -0.0000000000 0000000000 0000000091 67601d0/
        data bip1cs( 44) /  0.0000000000 0000000000 0000000020 40203d0/
        data bip1cs( 45) / -0.0000000000 0000000000 0000000003 73038d0/
        data bip1cs( 46) /  0.0000000000 0000000000 0000000000 46070d0/
        data bip1cs( 47) /  0.0000000000 0000000000 0000000000 01748d0/
  c
        data atr / 8.750690570 8484345088 0771988210 148d0 /
  c atr = 16./(sqrt(8)-1)  and  btr = -(sqrt(8)+1)/(sqrt(8)-1)
        data btr /-2.093836321 3560543136 0096498526 268d0 /
  c
        data nbif, nbig, nbif2, nbig2, nbip1, nbip2 / 6*0 /
        data x2sml, x3sml, x32sml, xbig / 4*0.0d0 /
  c
        if (nbif.ne.0) go to 10
        eta = 0.1*sngl(d1mach(3))
        nbif = initds (bifcs, 13, eta)
        nbig = initds (bigcs, 13, eta)
        nbif2 = initds (bif2cs, 15, eta)
        nbig2 = initds (big2cs, 16, eta)
        nbip1 = initds (bip1cs, 47, eta)
        nbip2 = initds (bip2cs, 88, eta)
  c
        x2sml = sqrt (eta)
        x3sml = eta**0.3333
        x32sml = 1.3104d0*x3sml**2
        xbig = d1mach(2)**0.6666d0
  c
   10   if (x.ge.(-1.0d0)) go to 20
        call d9admp (x, xn, phi)
        dbide = xn * dsin (phi)
        return
  c
   20   if (x.gt.1.0d0) go to 30
        x3 = 0.0d0
        if (dabs(x).gt.x3sml) x3 = x**3
        x2 = 0.0d0
        if (dabs(x).gt.x2sml) x2 = x*x
        dbide = x2*(dcsevl (x3, bifcs, nbif) + 0.25d0) +
       1  dcsevl (x3, bigcs, nbig) + 0.5d0
        if (x.gt.x32sml) dbide = dbide * dexp (-2.0d0*x*dsqrt(x)/3.0d0)
        return
  c
   30   if (x.gt.2.0d0) go to 40
        z = (2.0d0*x**3 - 9.0d0) / 7.0d0
        dbide = dexp (-2.0d0*x*dsqrt(x)/3.0d0) * (x*x*(0.25d0 +
       1  dcsevl (z, bif2cs, nbif2)) + 0.5d0 + dcsevl (z, big2cs, nbig2))
        return
  c
   40   if (x.gt.4.0d0) go to 50
        sqrtx = dsqrt (x)
        z = atr/(x*sqrtx) + btr
        dbide = (0.625d0 + dcsevl (z, bip1cs, nbip1)) * dsqrt(sqrtx)
        return
  c
   50   sqrtx = dsqrt (x)
        z = -1.0d0
        if (x.lt.xbig) z = 16.0d0/(x*sqrtx) - 1.0d0
        dbide = (0.625d0 + dcsevl (z, bip2cs, nbip2)) * dsqrt(sqrtx)
        return
  c
        end