dbid.f 6.72 KB
double precision function dbid (x)
c july 1980 edition.  w. fullerton, bell labs.
c
c evaluate the derivative of the airy function bi(x).
c
      double precision x, bifcs(13), bigcs(13), bif2cs(15), big2cs(16),
     1  x2sml, x3sml, xmax, x3, x2, z, xn, phi, d1mach, dbide, dcsevl,
     2  dexp, dlog, dsin, dsqrt
      external d1mach, dbide, dcsevl, initds
c     1  sqrt
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
      data nbif, nbig, nbif2, nbig2 / 4*0 /
      data x2sml, x3sml, xmax / 3*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)
c
      x2sml = sqrt (eta)
      x3sml = eta**0.3333
      xmax = (1.5d0*dlog(d1mach(2)))**0.6666d0
c
 10   if (x.ge.(-1.0d0)) go to 20
      call d9admp (x, xn, phi)
      dbid = 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
      dbid = x2*(dcsevl (x3, bifcs, nbif) + 0.25d0) +
     1  dcsevl (x3, bigcs, nbig) + 0.5d0
      return
c
 30   if (x.gt.2.0d0) go to 40
      z = (2.0d0*x**3 - 9.0d0) / 7.0d0
      dbid = x*x*(dcsevl (z, bif2cs, nbif2) + 0.25d0) +
     1  dcsevl (z, big2cs, nbig2) + 0.5d0
      return
c
 40   if (x.gt.xmax) call seteru (
     1  36hdbid    x so big that dbid overflows, 36, 1, 2)
c
      dbid = dbide(x) * dexp (2.0d0*x*dsqrt(x)/3.0d0)
      return
c
      end