subroutine d9knus (xnu, x, bknu, bknu1, iswtch)
c june 1977 edition.   w. fullerton, c3, los alamos scientific lab.
c
c compute bessel functions exp(x) * k-sub-xnu (x)  and
c exp(x) * k-sub-xnu+1 (x) for 0.0 .le. xnu .lt. 1.0 .
c
      double precision xnu, x, bknu, bknu1, alpha(32), beta(32), a(32),
     1  c0kcs(29), znu1cs(20), alnz, aln2, a0, bknud, bknu0,
     2  b0, c0, euler, expx, p1, p2, p3, qq, result, sqpi2, sqrtx, v,
     3  vlnz, xi, xmu, xnusml, xsml, x2n, x2tov, z, ztov, alnsml,
     4  alnbig
      real alneps, an, bn
      double precision d1mach, dcsevl, dexp, dgamma, dlog, dsqrt
      external d1mach, dcsevl, dgamma, initds
c
c series for c0k        on the interval  0.          to  2.50000e-01
c                                        with weighted error   2.16e-32
c                                         log weighted error  31.67
c                               significant figures required  30.86
c                                    decimal places required  32.40
c
      data c0k cs(  1) / +.6018305724 2626108387 5774451803 29 d-1     /
      data c0k cs(  2) / -.1536487143 3017286092 9597559431 24 d+0     /
      data c0k cs(  3) / -.1175117600 8210492040 0682292262 13 d-1     /
      data c0k cs(  4) / -.8524878889 1979509827 0484015509 87 d-3     /
      data c0k cs(  5) / -.6132983876 7496791874 0981769221 11 d-4     /
      data c0k cs(  6) / -.4405228124 5510444562 6798895485 05 d-5     /
      data c0k cs(  7) / -.3163124672 8384488192 9154458921 99 d-6     /
      data c0k cs(  8) / -.2271071938 2899588330 6737717933 96 d-7     /
      data c0k cs(  9) / -.1630564460 8077609552 2746205153 60 d-8     /
      data c0k cs( 10) / -.1170693929 9414776568 7560440431 30 d-9     /
      data c0k cs( 11) / -.8405206378 6464437174 5465934137 92 d-11    /
      data c0k cs( 12) / -.6034667011 8979991487 0960507371 98 d-12    /
      data c0k cs( 13) / -.4332696033 5681371952 0459973669 03 d-13    /
      data c0k cs( 14) / -.3110735803 0203546214 6346977722 37 d-14    /
      data c0k cs( 15) / -.2233407822 6736982254 4861334098 40 d-15    /
      data c0k cs( 16) / -.1603514671 6864226300 6357915286 10 d-16    /
      data c0k cs( 17) / -.1151271736 3666556196 0356977053 05 d-17    /
      data c0k cs( 18) / -.8265759174 6836959105 1694790892 58 d-19    /
      data c0k cs( 19) / -.5934548080 6383948172 3334366959 84 d-20    /
      data c0k cs( 20) / -.4260813819 6467143926 4996130239 76 d-21    /
      data c0k cs( 21) / -.3059126686 4812876299 2636983705 42 d-22    /
      data c0k cs( 22) / -.2196354142 6734575224 9755018155 16 d-23    /
      data c0k cs( 23) / -.1576911326 1495836071 1057506847 60 d-24    /
      data c0k cs( 24) / -.1132171393 5950320948 7577310480 56 d-25    /
      data c0k cs( 25) / -.8128624883 4598404082 7923497144 33 d-27    /
      data c0k cs( 26) / -.5836090089 3453226552 8293493159 49 d-28    /
      data c0k cs( 27) / -.4190124162 3610922519 4523377809 05 d-29    /
      data c0k cs( 28) / -.3008373796 0206435069 5305042128 62 d-30    /
      data c0k cs( 29) / -.2159915206 7808647728 3421680898 32 d-31    /
c
c series for znu1       on the interval -7.00000e-01 to  0.
c                                        with weighted error   2.45e-33
c                                         log weighted error  32.61
c                               significant figures required  31.85
c                                    decimal places required  33.26
c
      data znu1cs(  1) / +.2033067569 9419172967 4444001216 911 d+0    /
      data znu1cs(  2) / +.1400779334 1321977106 2943670790 563 d+0    /
      data znu1cs(  3) / +.7916796961 0016135284 0972241972 320 d-2    /
      data znu1cs(  4) / +.3398011825 3210404535 2930092205 750 d-3    /
      data znu1cs(  5) / +.1174197568 8989336666 4507228352 690 d-4    /
      data znu1cs(  6) / +.3393575706 1226168033 3825865475 121 d-6    /
      data znu1cs(  7) / +.8425941769 7621991019 4629891264 803 d-8    /
      data znu1cs(  8) / +.1833366770 2485008918 4748150900 090 d-9    /
      data znu1cs(  9) / +.3549698447 0441631086 3007064469 557 d-11   /
      data znu1cs( 10) / +.6190324964 6988733220 5244342078 407 d-13   /
      data znu1cs( 11) / +.9819645356 8043942496 0346115456 527 d-15   /
      data znu1cs( 12) / +.1428513143 9649047421 1473563005 985 d-16   /
      data znu1cs( 13) / +.1918949218 8782529896 6162467488 436 d-18   /
      data znu1cs( 14) / +.2394309797 3949891416 2313140597 128 d-20   /
      data znu1cs( 15) / +.2788902468 1534735483 5870465474 995 d-22   /
      data znu1cs( 16) / +.3046066506 3303344258 2845214092 865 d-24   /
      data znu1cs( 17) / +.3131732370 4219181577 1564260932 089 d-26   /
      data znu1cs( 18) / +.3041330989 8785495164 5174908005 034 d-28   /
      data znu1cs( 19) / +.2798403846 3683308434 3185097659 733 d-30   /
      data znu1cs( 20) / +.2446371862 7449759648 5238794922 666 d-32   /
c
      data euler / 0.5772156649 0153286060 6512090082 40d0 /
      data sqpi2 / +1.253314137 3155002512 0788264240 55 d0      /
      data aln2 / 0.6931471805 5994530941 7232121458 18d0 /
      data ntc0k, ntznu1 / 2*0 /
      data xnusml, xsml, alnsml, alnbig, alneps / 4*0.d0, 0.0 /
c
      if (ntc0k.ne.0) go to 10
      eta = 0.1d0*d1mach(3)
      ntc0k = initds (c0kcs, 29, eta)
      ntznu1 = initds (znu1cs, 20, eta)
c
      xnusml = dsqrt (d1mach(3)/8.d0)
      xsml = 0.1d0*d1mach(3)
      alnsml = dlog (d1mach(1))
      alnbig = dlog (d1mach(2))
      alneps = dlog (0.1d0*d1mach(3))
c
 10   if (xnu.lt.0.d0 .or. xnu.ge.1.d0) call seteru (
     1  33hd9knus  xnu must be ge 0 and lt 1, 33, 1, 2)
      if (x.le.0.) call seteru (22hd9knus  x must be gt 0, 22, 2, 2)
c
      iswtch = 0
      if (x.gt.2.0d0) go to 50
c
c x is small.  compute k-sub-xnu (x) and the derivative of k-sub-xnu (x)
c then find k-sub-xnu+1 (x).  xnu is reduced to the interval (-.5,+.5)
c then to (0., .5), because k of negative order (-nu) = k of positive
c order (+nu).
c
      v = xnu
      if (xnu.gt.0.5d0) v = 1.0d0 - xnu
c
c carefully find (x/2)**xnu and z**xnu where z = x*x/4.
      alnz = 2.d0 * (dlog(x) - aln2)
c
      if (x.gt.xnu) go to 20
      if (-0.5d0*xnu*alnz-aln2-dlog(xnu).gt.alnbig) call seteru (
     1  45hd9knus  x so small bessel k-sub-xnu overflows, 45, 3, 2)
c
 20   vlnz = v*alnz
      x2tov = dexp (0.5d0*vlnz)
      ztov = 0.0d0
      if (vlnz.gt.alnsml) ztov = x2tov**2
c
      a0 = 0.5d0*dgamma(1.0d0+v)
      b0 = 0.5d0*dgamma(1.0d0-v)
      c0 = -euler
      if (ztov.gt.0.5d0 .and. v.gt.xnusml) c0 = -0.75d0 +
     1  dcsevl ((8.0d0*v)*v-1.0d0, c0kcs, ntc0k)
c
      if (ztov.le.0.5d0) alpha(1) = (a0-ztov*b0)/v
      if (ztov.gt.0.5d0) alpha(1) = c0 - alnz*(0.75d0 +
     1  dcsevl (vlnz/0.35d0+1.0d0, znu1cs, ntznu1))*b0
      beta(1) = -0.5d0*(a0+ztov*b0)
c
      z = 0.0d0
      if (x.gt.xsml) z = 0.25d0*x*x
      nterms = max1 (2.0, 11.0+(8.*sngl(alnz)-25.19-alneps)
     1  /(4.28-sngl(alnz)))
      do 30 i=2,nterms
        xi = i - 1
        a0 = a0/(xi*(xi-v))
        b0 = b0/(xi*(xi+v))
        alpha(i) = (alpha(i-1)+2.0d0*xi*a0)/(xi*(xi+v))
        beta(i) = (xi-0.5d0*v)*alpha(i) - ztov*b0
 30   continue
c
      bknu = alpha(nterms)
      bknud = beta(nterms)
      do 40 ii=2,nterms
        i = nterms + 1 - ii
        bknu = alpha(i) + bknu*z
        bknud = beta(i) + bknud*z
 40   continue
c
      expx = dexp(x)
      bknu = expx*bknu/x2tov
c
      if (-0.5d0*(xnu+1.d0)*alnz-2.0d0*aln2.gt.alnbig) iswtch = 1
      if (iswtch.eq.1) return
      bknud = expx*bknud*2.0d0/(x2tov*x)
c
      if (xnu.le.0.5d0) bknu1 = v*bknu/x - bknud
      if (xnu.le.0.5d0) return
c
      bknu0 = bknu
      bknu = -v*bknu/x - bknud
      bknu1 = 2.0d0*xnu*bknu/x + bknu0
      return
c
c x is large.  find k-sub-xnu (x) and k-sub-xnu+1 (x) with y. l. luke-s
c rational expansion.
c
 50   sqrtx = dsqrt(x)
      if (x.gt.1.0d0/xsml) go to 90
      an = -0.60 - 1.02/sngl(x)
      bn = -0.27 - 0.53/sngl(x)
      nterms = min0 (32, max1 (3.0, an+bn*alneps))
c
      do 80 inu=1,2
        xmu = 0.d0
        if (inu.eq.1 .and. xnu.gt.xnusml) xmu = (4.0d0*xnu)*xnu
        if (inu.eq.2) xmu = 4.0d0*(dabs(xnu)+1.d0)**2
c
        a(1) = 1.0d0 - xmu
        a(2) = 9.0d0 - xmu
        a(3) = 25.0d0 - xmu
        if (a(2).eq.0.d0) result = sqpi2*(16.d0*x+xmu+7.d0) /
     1    (16.d0*x*sqrtx)
        if (a(2).eq.0.d0) go to 70
c
        alpha(1) = 1.0d0
        alpha(2) = (16.d0*x+a(2))/a(2)
        alpha(3) = ((768.d0*x+48.d0*a(3))*x + a(2)*a(3))/(a(2)*a(3))
c
        beta(1) = 1.0d0
        beta(2) = (16.d0*x+(xmu+7.d0))/a(2)
        beta(3) = ((768.d0*x+48.d0*(xmu+23.d0))*x +
     1    ((xmu+62.d0)*xmu+129.d0))/(a(2)*a(3))
c
        if (nterms.lt.4) go to 65
        do 60 i=4,nterms
          n = i - 1
          x2n = 2*n - 1
c
          a(i) = (x2n+2.d0)**2 - xmu
          qq = 16.d0*x2n/a(i)
          p1 = -x2n*(dble(float(12*n*n-20*n))-a(1))/((x2n-2.d0)*a(i))
     1      - qq*x
          p2 = (dble(float(12*n*n-28*n+8))-a(1))/a(i) - qq*x
          p3 = -x2n*a(i-3)/((x2n-2.d0)*a(i))
c
          alpha(i) = -p1*alpha(i-1) - p2*alpha(i-2) - p3*alpha(i-3)
          beta(i) = -p1*beta(i-1) - p2*beta(i-2) - p3*beta(i-3)
 60     continue
c
 65     result = sqpi2*beta(nterms)/(sqrtx*alpha(nterms))
c
 70     if (inu.eq.1) bknu = result
        if (inu.eq.2) bknu1 = result
 80   continue
      return
c
 90   bknu = sqpi2/sqrtx
      bknu1 = bknu
      return
c
      end