dbide.f 17.3 KB
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