bid.f 4.64 KB
function bid (x)
c july 1980 edition.  w. fullerton, bell labs.
c
c evaluate the derivative of the airy function bi(x).
c
      dimension bifcs(8), bigcs(9), bif2cs(10), big2cs(10)
      external  bide, csevl, inits, r1mach
c
c series for bif on the interval -1.00000e+00 to  1.00000e+00
c                                        with weighted error   9.05e-18
c                                         log weighted error  17.04
c                               significant figures required  15.83
c                                    decimal places required  17.49
c
      data bifcs(  1) /     0.1153536790 828570243 e0/
      data bifcs(  2) /     0.0205007894 049192875 e0/
      data bifcs(  3) /     0.0002135290 278902876 e0/
      data bifcs(  4) /     0.0000010783 960614677 e0/
      data bifcs(  5) /     0.0000000032 094708833 e0/
      data bifcs(  6) /     0.0000000000 062930407 e0/
      data bifcs(  7) /     0.0000000000 000087403 e0/
      data bifcs(  8) /     0.0000000000 000000090 e0/
c
c series for big on the interval -1.00000e+00 to  1.00000e+00
c                                        with weighted error   5.44e-19
c                                         log weighted error  18.26
c                               significant figures required  17.46
c                                    decimal places required  18.74
c
      data bigcs(  1) /    -0.0971964404 1644353739 0e0/
      data bigcs(  2) /     0.1495035768 4316706657 1e0/
      data bigcs(  3) /     0.0031135253 8712132604 2e0/
      data bigcs(  4) /     0.0000247085 7057982129 7e0/
      data bigcs(  5) /     0.0000001029 4962773137 9e0/
      data bigcs(  6) /     0.0000000002 6397037398 7e0/
      data bigcs(  7) /     0.0000000000 0045827927 1e0/
      data bigcs(  8) /     0.0000000000 0000057428 3e0/
      data bigcs(  9) /     0.0000000000 0000000054 4e0/
c
c series for bif2 on the interval  1.00000e+00 to  8.00000e+00
c                                        with weighted error   3.82e-19
c                                         log weighted error  18.42
c                               significant figures required  17.68
c                                    decimal places required  18.92
c
      data bif2cs(  1) /     0.3234939876 0352203352 1e0/
      data bif2cs(  2) /     0.0862978715 3556355913 9e0/
      data bif2cs(  3) /     0.0029940255 5265539742 6e0/
      data bif2cs(  4) /     0.0000514305 2836466163 7e0/
      data bif2cs(  5) /     0.0000005258 4025003681 1e0/
      data bif2cs(  6) /     0.0000000035 6175137395 8e0/
      data bif2cs(  7) /     0.0000000000 1714686400 7e0/
      data bif2cs(  8) /     0.0000000000 0006166352 0e0/
      data bif2cs(  9) /     0.0000000000 0000017191 1e0/
      data bif2cs( 10) /     0.0000000000 0000000038 2e0/
c
c series for big2 on the interval  1.00000e+00 to  8.00000e+00
c                                        with weighted error   3.35e-17
c                                         log weighted error  16.48
c                               significant figures required  16.52
c                                    decimal places required  16.98
c
      data big2cs(  1) /     1.6062999463 621294578 e0/
      data big2cs(  2) /     0.7449088819 876088652 e0/
      data big2cs(  3) /     0.0470138738 610277380 e0/
      data big2cs(  4) /     0.0012284422 062548239 e0/
      data big2cs(  5) /     0.0000173222 412256624 e0/
      data big2cs(  6) /     0.0000001521 901652368 e0/
      data big2cs(  7) /     0.0000000009 113560249 e0/
      data big2cs(  8) /     0.0000000000 039547918 e0/
      data big2cs(  9) /     0.0000000000 000130017 e0/
      data big2cs( 10) /     0.0000000000 000000335 e0/
c
      data nbif, nbig, nbif2, nbig2 / 4*0 /
      data x2sml, x3sml, xmax / 3*0.0 /
c
      if (nbif.ne.0) go to 10
      eta = 0.1*r1mach(3)
      nbif = inits (bifcs, 8, eta)
      nbig = inits (bigcs, 9, eta)
      nbif2 = inits (bif2cs, 10, eta)
      nbig2 = inits (big2cs, 10, eta)
c
      x2sml = sqrt (eta)
      x3sml = eta**0.3333
      xmax = (1.5*alog(r1mach(2)))**0.6666
c
 10   if (x.ge.(-1.0)) go to 20
      call r9admp (x, xn, phi)
      bid = xn * sin (phi)
      return
c
 20   if (x.gt.1.0) go to 30
      x3 = 0.0
      if (abs(x).gt.x3sml) x3 = x**3
      x2 = 0.0
      if (abs(x).gt.x2sml) x2 = x*x
      bid = x2*(csevl (x3, bifcs, nbif) + 0.25) +
     1  csevl (x3, bigcs, nbig) + 0.5
      return
c
 30   if (x.gt.2.0) go to 40
      z = (2.0*x**3 - 9.0) / 7.0
      bid = x*x*(csevl (z, bif2cs, nbif2) + 0.25) +
     1  csevl (z, big2cs, nbig2) + 0.5
      return
c
 40   if (x.gt.xmax) call seteru (
     1  35hbid     x so big that bid overflows, 35, 1, 2)
c
      bid = bide(x) * exp (2.0*x*sqrt(x)/3.0)
c
      return
      end