bide.f 8.86 KB
function bide (x)
c july 1980 edition.  w. fullerton, bell labs.
c
c evaluate the dervative of the airy function bi(x) for x.lt.0
c and evaluate bid(x) * exp(-zeta) for x .ge. 0 where
c zeta = 2/3 * x**(3/2)
c
      dimension bifcs(8), bigcs(9), bif2cs(10), big2cs(10), bip1cs(24),
     1  bip2cs(29)
      external 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
c series for bip2 on the interval  0.00000e+00 to  1.25000e-01
c                                        with weighted error   2.07e-18
c                                         log weighted error  17.69
c                               significant figures required  16.51
c                                    decimal places required  18.42
c
      data bip2cs(  1) /    -0.1326970544 3526630495 e0/
      data bip2cs(  2) /    -0.0056844362 6045977481 e0/
      data bip2cs(  3) /    -0.0001564360 1119611610 e0/
      data bip2cs(  4) /    -0.0000113673 7203679562 e0/
      data bip2cs(  5) /    -0.0000014346 4350991284 e0/
      data bip2cs(  6) /    -0.0000001809 8531185164 e0/
      data bip2cs(  7) /     0.0000000092 6177343611 e0/
      data bip2cs(  8) /     0.0000000171 0005490721 e0/
      data bip2cs(  9) /     0.0000000047 6698163504 e0/
      data bip2cs( 10) /    -0.0000000003 5195022023 e0/
      data bip2cs( 11) /    -0.0000000005 8890614316 e0/
      data bip2cs( 12) /    -0.0000000000 6678499608 e0/
      data bip2cs( 13) /     0.0000000000 6395565102 e0/
      data bip2cs( 14) /     0.0000000000 1554529427 e0/
      data bip2cs( 15) /    -0.0000000000 0792397000 e0/
      data bip2cs( 16) /    -0.0000000000 0258326243 e0/
      data bip2cs( 17) /     0.0000000000 0121655048 e0/
      data bip2cs( 18) /     0.0000000000 0038707207 e0/
      data bip2cs( 19) /    -0.0000000000 0022487045 e0/
      data bip2cs( 20) /    -0.0000000000 0004953477 e0/
      data bip2cs( 21) /     0.0000000000 0004563782 e0/
      data bip2cs( 22) /     0.0000000000 0000332998 e0/
      data bip2cs( 23) /    -0.0000000000 0000921750 e0/
      data bip2cs( 24) /     0.0000000000 0000094157 e0/
      data bip2cs( 25) /     0.0000000000 0000167154 e0/
      data bip2cs( 26) /    -0.0000000000 0000055134 e0/
      data bip2cs( 27) /    -0.0000000000 0000022369 e0/
      data bip2cs( 28) /     0.0000000000 0000017487 e0/
      data bip2cs( 29) /     0.0000000000 0000000207 e0/
c
c series for bip1 on the interval  1.25000e-01 to  3.53553e-01
c                                        with weighted error   1.86e-17
c                                         log weighted error  16.73
c                               significant figures required  15.67
c                                    decimal places required  17.42
c
      data bip1cs(  1) /    -0.1729187351 079553719 e0/
      data bip1cs(  2) /    -0.0149358492 984694364 e0/
      data bip1cs(  3) /    -0.0005471104 951678566 e0/
      data bip1cs(  4) /     0.0001537966 292958408 e0/
      data bip1cs(  5) /     0.0000154353 476192179 e0/
      data bip1cs(  6) /    -0.0000065434 113851906 e0/
      data bip1cs(  7) /     0.0000003728 082407879 e0/
      data bip1cs(  8) /     0.0000002072 078388189 e0/
      data bip1cs(  9) /    -0.0000000658 173336470 e0/
      data bip1cs( 10) /     0.0000000074 926746354 e0/
      data bip1cs( 11) /     0.0000000011 101336884 e0/
      data bip1cs( 12) /    -0.0000000007 265140553 e0/
      data bip1cs( 13) /     0.0000000001 782723560 e0/
      data bip1cs( 14) /    -0.0000000000 217346352 e0/
      data bip1cs( 15) /    -0.0000000000 020302035 e0/
      data bip1cs( 16) /     0.0000000000 019311827 e0/
      data bip1cs( 17) /    -0.0000000000 006044953 e0/
      data bip1cs( 18) /     0.0000000000 001209450 e0/
      data bip1cs( 19) /    -0.0000000000 000125109 e0/
      data bip1cs( 20) /    -0.0000000000 000019917 e0/
      data bip1cs( 21) /     0.0000000000 000015154 e0/
      data bip1cs( 22) /    -0.0000000000 000004977 e0/
      data bip1cs( 23) /     0.0000000000 000001155 e0/
      data bip1cs( 24) /    -0.0000000000 000000186 e0/
c
      data atr / 8.750690570 8484345 e0 /
c atr = 16./(sqrt(8)-1)  and  btr = -(sqrt(8)+1)/(sqrt(8)-1)
      data btr /-2.093836321 3560543 e0 /
c
      data nbif, nbig, nbif2, nbig2, nbip1, nbip2 / 6*0 /
      data x2sml, x3sml, x32sml, xbig / 4*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)
      nbip1 = inits (bip1cs, 24, eta)
      nbip2 = inits (bip2cs, 29, eta)
c
      x2sml = sqrt (eta)
      x3sml = eta**0.3333
      x32sml = 1.3104*x3sml**2
      xbig = r1mach(2)**0.6666
c
 10   if (x.ge.(-1.0)) go to 20
      call r9admp (x, xn, phi)
      bide = 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
      bide = x2 * (csevl(x3, bifcs, nbif) + 0.25) +
     1  csevl (x3, bigcs, nbig) + 0.5
      if (x.gt.x32sml) bide = bide * exp (-2.0*x*sqrt(x)/3.0)
      return
c
 30   if (x.gt.2.0) go to 40
      z = (2.0*x**3 - 9.0) / 7.0
      bide = exp (-2.0*x*sqrt(x)/3.0) * (x*x * (0.25 +
     1  csevl (z, bif2cs, nbif2)) + 0.5 + csevl (z, big2cs, nbig2))
      return
c
 40   if (x.gt.4.0) go to 50
      sqrtx = sqrt(x)
      z = atr/(x*sqrtx) + btr
      bide = (0.625 + csevl (z, bip1cs, nbip1)) * sqrt(sqrtx)
      return
c
 50   sqrtx = sqrt(x)
      z = -1.0
      if (x.lt.xbig) z = 16.0/(x*sqrtx) - 1.0
      bide = (0.625 + csevl (z, bip2cs, nbip2)) * sqrt(sqrtx)
      return
c
      end