Blame view
fvn_fnlib/dbide.f
17.3 KB
38581db0c 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 |