double precision function dbesk0 (x)
c july 1977 edition.  w. fullerton, c3, los alamos scientific lab.
      double precision x, bk0cs(16), xmax, xsml, y,
     1  d1mach, dcsevl, dbesi0, dbsk0e, dexp, dlog, dsqrt
      external d1mach, dbesi0, dbsk0e, dcsevl, dexp, dlog, dsqrt, initds
c
c series for bk0        on the interval  0.          to  4.00000e+00
c                                        with weighted error   3.08e-33
c                                         log weighted error  32.51
c                               significant figures required  32.05
c                                    decimal places required  33.11
c
      data bk0 cs(  1) / -.3532739323 3902768720 1140060063 153 d-1    /
      data bk0 cs(  2) / +.3442898999 2462848688 6344927529 213 d+0    /
      data bk0 cs(  3) / +.3597993651 5361501626 5721303687 231 d-1    /
      data bk0 cs(  4) / +.1264615411 4469259233 8479508673 447 d-2    /
      data bk0 cs(  5) / +.2286212103 1194517860 8269830297 585 d-4    /
      data bk0 cs(  6) / +.2534791079 0261494573 0790013428 354 d-6    /
      data bk0 cs(  7) / +.1904516377 2202088589 7214059381 366 d-8    /
      data bk0 cs(  8) / +.1034969525 7633624585 1008317853 089 d-10   /
      data bk0 cs(  9) / +.4259816142 7910825765 2445327170 133 d-13   /
      data bk0 cs( 10) / +.1374465435 8807508969 4238325440 000 d-15   /
      data bk0 cs( 11) / +.3570896528 5083735909 9688597333 333 d-18   /
      data bk0 cs( 12) / +.7631643660 1164373766 7498666666 666 d-21   /
      data bk0 cs( 13) / +.1365424988 4407818590 8053333333 333 d-23   /
      data bk0 cs( 14) / +.2075275266 9066680831 9999999999 999 d-26   /
      data bk0 cs( 15) / +.2712814218 0729856000 0000000000 000 d-29   /
      data bk0 cs( 16) / +.3082593887 9146666666 6666666666 666 d-32   /
c
      data ntk0, xsml, xmax / 0, 2*0.d0 /
c
      if (ntk0.ne.0) go to 10
      ntk0 = initds (bk0cs, 16, 0.1*sngl(d1mach(3)))
      xsml = dsqrt (4.0d0*d1mach(3))
      xmax = -dlog(d1mach(1))
      xmax = xmax - 0.5d0*xmax*dlog(xmax)/(xmax+0.5d0)
c
 10   if (x.le.0.d0) call seteru (29hdbesk0  x is zero or negative,
     1  29, 2, 2)
      if (x.gt.2.0d0) go to 20
c
      y = 0.d0
      if (x.gt.xsml) y = x*x
      dbesk0 = -dlog(0.5d0*x)*dbesi0(x) - 0.25d0 + dcsevl (.5d0*y-1.d0,
     1  bk0cs, ntk0)
      return
c
 20   dbesk0 = 0.d0
      if (x.gt.xmax) call seteru (30hdbesk0  x so big k0 underflows,
     1  30, 1, 0)
      if (x.gt.xmax) return
c
      dbesk0 = dexp(-x) * dbsk0e(x)
c
      return
      end