function besi1 (x)
c oct 1983 version.  w. fullerton, c3, los alamos scientific lab.
      dimension bi1cs(11)
      external alog, besi1e, csevl, exp, inits, r1mach, sqrt
c
c series for bi1        on the interval  0.          to  9.00000d+00
c                                        with weighted error   2.40e-17
c                                         log weighted error  16.62
c                               significant figures required  16.23
c                                    decimal places required  17.14
c
      data bi1 cs( 1) /   -.0019717132 61099859e0 /
      data bi1 cs( 2) /    .4073488766 7546481e0 /
      data bi1 cs( 3) /    .0348389942 99959456e0 /
      data bi1 cs( 4) /    .0015453945 56300123e0 /
      data bi1 cs( 5) /    .0000418885 21098377e0 /
      data bi1 cs( 6) /    .0000007649 02676483e0 /
      data bi1 cs( 7) /    .0000000100 42493924e0 /
      data bi1 cs( 8) /    .0000000000 99322077e0 /
      data bi1 cs( 9) /    .0000000000 00766380e0 /
      data bi1 cs(10) /    .0000000000 00004741e0 /
      data bi1 cs(11) /    .0000000000 00000024e0 /
c
c
      data nti1, xmin, xsml, xmax / 0, 3*0. /
c
      if (nti1.ne.0) go to 10
      nti1 = inits (bi1cs, 11, 0.1*r1mach(3))
      xmin = 2.0*r1mach(1)
      xsml = sqrt (8.0*r1mach(3))
      xmax = alog (r1mach(2))
c
 10   y = abs(x)
      if (y.gt.3.0) go to 20
c
      besi1 = 0.0
      if (y.eq.0.0) return
c
      if (y.le.xmin) call seteru (
     1  37hbesi1   abs(x) so small i1 underflows, 37, 1, 0)
      if (y.gt.xmin) besi1 = 0.5*x
      if (y.gt.xsml) besi1 = x * (.875 + csevl(y*y/4.5-1., bi1cs, nti1))
      return
c
 20   if (y.gt.xmax) call seteru (
     1  34hbesi1   abs(x) so big i1 overflows, 34, 2, 2)
c
      besi1 = exp(y) * besi1e(x)
c
      return
      end