      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 alog, csevl, exp, gamma, inits, r1mach, sqrt
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