subroutine r9knus (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
        dimension alpha(15), beta(15), a(15), c0kcs(16), znu1cs(12)
        external csevl,gamma, inits, r1mach
  c
  c series for c0k        on the interval  0.          to  2.50000d-01
  c                                        with weighted error   1.60e-17
  c                                         log weighted error  16.79
  c                               significant figures required  15.99
  c                                    decimal places required  17.40
  c
        data c0k cs( 1) /    .0601830572 42626108e0 /
        data c0k cs( 2) /   -.1536487143 3017286e0 /
        data c0k cs( 3) /   -.0117511760 08210492e0 /
        data c0k cs( 4) /   -.0008524878 88919795e0 /
        data c0k cs( 5) /   -.0000613298 38767496e0 /
        data c0k cs( 6) /   -.0000044052 28124551e0 /
        data c0k cs( 7) /   -.0000003163 12467283e0 /
        data c0k cs( 8) /   -.0000000227 10719382e0 /
        data c0k cs( 9) /   -.0000000016 30564460e0 /
        data c0k cs(10) /   -.0000000001 17069392e0 /
        data c0k cs(11) /   -.0000000000 08405206e0 /
        data c0k cs(12) /   -.0000000000 00603466e0 /
        data c0k cs(13) /   -.0000000000 00043326e0 /
        data c0k cs(14) /   -.0000000000 00003110e0 /
        data c0k cs(15) /   -.0000000000 00000223e0 /
        data c0k cs(16) /   -.0000000000 00000016e0 /
  c
  c series for znu1       on the interval -7.00000d-01 to  0.
  c                                        with weighted error   1.43e-17
  c                                         log weighted error  16.85
  c                               significant figures required  16.08
  c                                    decimal places required  17.38
  c
        data znu1cs( 1) /    .2033067569 9419173e0 /
        data znu1cs( 2) /    .1400779334 1321977e0 /
        data znu1cs( 3) /    .0079167969 61001613e0 /
        data znu1cs( 4) /    .0003398011 82532104e0 /
        data znu1cs( 5) /    .0000117419 75688989e0 /
        data znu1cs( 6) /    .0000003393 57570612e0 /
        data znu1cs( 7) /    .0000000084 25941769e0 /
        data znu1cs( 8) /    .0000000001 83336677e0 /
        data znu1cs( 9) /    .0000000000 03549698e0 /
        data znu1cs(10) /    .0000000000 00061903e0 /
        data znu1cs(11) /    .0000000000 00000981e0 /
        data znu1cs(12) /    .0000000000 00000014e0 /
  c
        data euler / 0.5772156649 0153286e0 /
        data sqpi2 / 1.253314137 3155003e0 /
        data aln2 / 0.693147180 55994531e0 /
        data ntc0k, ntznu1 / 2*0 /
        data xnusml, xsml, alnsml, alnbig, alneps / 5*0.0 /
  c
        if (ntc0k.ne.0) go to 10
        ntc0k = inits (c0kcs, 16, 0.1*r1mach(3))
        ntznu1 = inits (znu1cs, 12, 0.1*r1mach(3))
  c
        xnusml = sqrt (r1mach(3)/8.0)
        xsml = 0.1*r1mach(3)
        alnsml = alog (r1mach(1))
        alnbig = alog (r1mach(2))
        alneps = alog (0.1*r1mach(3))
  c
   10   if (xnu.lt.0. .or. xnu.ge.1.0) call seteru (
       1  33hr9knus  xnu must be ge 0 and lt 1, 33, 1, 2)
        if (x.le.0.) call seteru (22hr9knus  x must be gt 0, 22, 2, 2)
  c
        iswtch = 0
        if (x.gt.2.0) 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.5) v = 1.0 - xnu
  c
  c carefully find (x/2)**xnu and z**xnu where z = x*x/4.
        alnz = 2.0 * (alog(x) - aln2)
  c
        if (x.gt.xnu) go to 20
        if (-0.5*xnu*alnz-aln2-alog(xnu).gt.alnbig) call seteru (
       1  45hr9knus  x so small bessel k-sub-xnu overflows, 45, 3, 2)
  c
   20   vlnz = v*alnz
        x2tov = exp (0.5*vlnz)
        ztov = 0.0
        if (vlnz.gt.alnsml) ztov = x2tov**2
  c
        a0 = 0.5*gamma(1.0+v)
        b0 = 0.5*gamma(1.0-v)
        c0 = -euler
        if (ztov.gt.0.5 .and. v.gt.xnusml) c0 = -0.75 +
       1  csevl ((8.0*v)*v-1., c0kcs, ntc0k)
  c
        if (ztov.le.0.5) alpha(1) = (a0-ztov*b0)/v
        if (ztov.gt.0.5) alpha(1) = c0 - alnz*(0.75 +
       1  csevl (vlnz/0.35+1.0, znu1cs, ntznu1))*b0
        beta(1) = -0.5*(a0+ztov*b0)
  c
        z = 0.0
        if (x.gt.xsml) z = 0.25*x*x
        nterms = max1 (2.0, 11.0+(8.*alnz-25.19-alneps)/(4.28-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.0*xi*a0)/(xi*(xi+v))
          beta(i) = (xi-0.5*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 = exp(x)
        bknu = expx*bknu/x2tov
  c
        if (-0.5*(xnu+1.)*alnz-2.0*aln2.gt.alnbig) iswtch = 1
        if (iswtch.eq.1) return
        bknud = expx*bknud*2.0/(x2tov*x)
  c
        if (xnu.le.0.5) bknu1 = v*bknu/x - bknud
        if (xnu.le.0.5) return
  c
        bknu0 = bknu
        bknu = -v*bknu/x - bknud
        bknu1 = 2.0*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 = sqrt(x)
        if (x.gt.1.0/xsml) go to 90
        an = -1.56 + 4.0/x
        bn = -0.29 - 0.22/x
        nterms = min0 (15, max1 (3., an+bn*alneps))
  c
        do 80 inu=1,2
          xmu = 0.
          if (inu.eq.1 .and. xnu.gt.xnusml) xmu = (4.0*xnu)*xnu
          if (inu.eq.2) xmu = 4.0*(abs(xnu)+1.)**2
  c
          a(1) = 1.0 - xmu
          a(2) = 9.0 - xmu
          a(3) = 25.0 - xmu
          if (a(2).eq.0.) result = sqpi2*(16.*x+xmu+7.)/(16.*x*sqrtx)
          if (a(2).eq.0.) go to 70
  c
          alpha(1) = 1.0
          alpha(2) = (16.*x+a(2))/a(2)
          alpha(3) = ((768.*x+48.*a(3))*x + a(2)*a(3))/(a(2)*a(3))
  c
          beta(1) = 1.0
          beta(2) = (16.*x+(xmu+7.))/a(2)
          beta(3) = ((768.*x+48.*(xmu+23.))*x + ((xmu+62.)*xmu+129.))
       1    / (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.)**2 - xmu
            qq = 16.*x2n/a(i)
            p1 = -x2n*(float(12*n*n-20*n)-a(1))/((x2n-2.)*a(i)) - qq*x
            p2 = (float(12*n*n-28*n+8)-a(1))/a(i) - qq*x
            p3 = -x2n*a(i-3)/((x2n-2.)*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