subroutine dbesks (xnu, x, nin, bk) c july 1980 edition. w. fullerton, c3, los alamos scientific lab. double precision xnu, x, bk(1), expxi, xmax, d1mach, dexp, dlog external d1mach, dexp, dlog data xmax / 0.d0 / c if (xmax.ne.0.0d0) go to 10 xmax = -dlog(d1mach(1)) xmax = xmax + 0.5d0*dlog(3.14d0*0.5d0/xmax) c 10 if (x.gt.xmax) call seteru ( 1 36hdbesks x so big bessel k underflows, 36, 1, 2) c call dbskes (xnu, x, nin, bk) c expxi = dexp (-x) n = iabs (nin) do 20 i=1,n bk(i) = expxi * bk(i) 20 continue c return end