iknb Subroutine

subroutine iknb ( n, x, nm, bi, di, bk, dk )

  !*****************************************************************************80
  !
  !! IKNB compute Bessel function In(x) and Kn(x).
  !
  !  Discussion:
  !
  !    Compute modified Bessel functions In(x) and Kn(x),
  !    and their derivatives.
  !
  !  Licensing:
  !
  !    This routine is copyrighted by Shanjie Zhang and Jianming Jin.  However, 
  !    they give permission to incorporate this routine into a user program 
  !    provided that the copyright is acknowledged.
  !
  !  Modified:
  !
  !    17 July 2012
  !
  !  Author:
  !
  !    Shanjie Zhang, Jianming Jin
  !
  !  Reference:
  !
  !    Shanjie Zhang, Jianming Jin,
  !    Computation of Special Functions,
  !    Wiley, 1996,
  !    ISBN: 0-471-11963-6,
  !    LC: QA351.C45.
  !
  !  Parameters:
  !
  !    Input, integer ( kind = 4 ) N, the order of In(x) and Kn(x).
  !
  !    Input, real ( kind = 8 ) X, the argument.
  !
  !    Output, integer ( kind = 4 ) NM, the highest order computed.
  !
  !    Output, real ( kind = 8 ) BI(0:N), DI(0:N), BK(0:N), DK(0:N),
  !    the values of In(x), In'(x), Kn(x), Kn'(x).
  !
  implicit none

  integer ( kind = 4 ) n

  real ( kind = 8 ) a0
  real ( kind = 8 ) bi(0:n)
  real ( kind = 8 ) bk(0:n)
  real ( kind = 8 ) bkl
  real ( kind = 8 ) bs
  real ( kind = 8 ) di(0:n)
  real ( kind = 8 ) dk(0:n)
  real ( kind = 8 ) el
  real ( kind = 8 ) f
  real ( kind = 8 ) f0
  real ( kind = 8 ) f1
  real ( kind = 8 ) g
  real ( kind = 8 ) g0
  real ( kind = 8 ) g1
  integer ( kind = 4 ) k
  integer ( kind = 4 ) k0
  integer ( kind = 4 ) l
  integer ( kind = 4 ) m
    ! integer ( kind = 4 ) msta1
  ! integer ( kind = 4 ) msta2
  integer ( kind = 4 ) nm
  real ( kind = 8 ) pi
  real ( kind = 8 ) r
  real ( kind = 8 ) s0
  real ( kind = 8 ) sk0
  real ( kind = 8 ) vt
  real ( kind = 8 ) x

  pi = 3.141592653589793D+00
  el = 0.5772156649015329d0
  nm = n

  if ( x <= 1.0D-100 ) then
     do k = 0, n
        bi(k) = 0.0D+00
        di(k) = 0.0D+00
        bk(k) = 1.0D+300
        dk(k) = -1.0D+300
     end do
     bi(0) = 1.0D+00
     di(1) = 0.5D+00
     return
  end if

  if ( n == 0 ) then
     nm = 1
  end if

  m = msta1 ( x, 200 )
  if ( m < nm ) then
     nm = m
  else
     m = msta2 ( x, nm, 15 )
  end if

  bs = 0.0D+00
  sk0 = 0.0D+00
  f0 = 0.0D+00
  f1 = 1.0D-100
  do k = m, 0, -1
     f = 2.0D+00 * ( k + 1.0D+00 ) / x * f1 + f0
     if ( k <= nm ) then
        bi(k) = f
     end if
     if ( k /= 0 .and. k == 2 * int ( k / 2 ) ) then
        sk0 = sk0 + 4.0D+00 * f / k
     end if
     bs = bs + 2.0D+00 * f
     f0 = f1
     f1 = f
  end do

  s0 = exp ( x ) / ( bs - f )
  do k = 0, nm
     bi(k) = s0 * bi(k)
  end do

  if ( x <= 8.0D+00 ) then
     bk(0) = - ( log ( 0.5D+00 * x ) + el ) * bi(0) + s0 * sk0
     bk(1) = ( 1.0D+00 / x - bi(1) * bk(0) ) / bi(0)
  else
     a0 = sqrt ( pi / ( 2.0D+00 * x ) ) * exp ( - x ) 

     if ( x < 25.0D+00 ) then
        k0 = 16
     else if ( x < 80.0D+00 ) then
        k0 = 10
     else if ( x < 200.0D+00 ) then
        k0 = 8
     else
        k0 = 6
     end if

     do l = 0, 1
        bkl = 1.0D+00
        vt = 4.0D+00 * l
        r = 1.0D+00
        do k = 1, k0
           r = 0.125D+00 * r * ( vt - ( 2.0D+00 * k - 1.0D+00 ) ** 2 ) / ( k * x )
           bkl = bkl + r
        end do
        bk(l) = a0 * bkl
     end do
  end if

  g0 = bk(0)
  g1 = bk(1)
  do k = 2, nm
     g = 2.0D+00 * ( k - 1.0D+00 ) / x * g1 + g0
     bk(k) = g
     g0 = g1
     g1 = g
  end do

  di(0) = bi(1)
  dk(0) = -bk(1)
  do k = 1, nm
     di(k) = bi(k-1) - k / x * bi(k)
     dk(k) = -bk(k-1) - k / x * bk(k)
  end do

  return
end subroutine iknb