sphk Subroutine

subroutine sphk ( n, x, nm, sk, dk )

  !*****************************************************************************80
  !
  !! SPHK computes modified spherical Bessel functions kn(x) and derivatives.
  !
  !  Discussion:
  !
  !    This procedure computes modified spherical Bessel functions
  !    of the second kind, kn(x) and kn'(x).
  !
  !  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:
  !
  !    15 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.
  !
  !    Input, real ( kind = 8 ) X, the argument.
  !
  !    Output, integer ( kind = 4 ) NM, the highest order computed.
  !
  !    Output, real ( kind = 8 ) SK(0:N), DK(0:N), the values of kn(x) and kn'(x).
  !
  implicit none

  integer ( kind = 4 ) n

  real ( kind = 8 ) dk(0:n)
  real ( kind = 8 ) f
  real ( kind = 8 ) f0
  real ( kind = 8 ) f1
  integer ( kind = 4 ) k
  integer ( kind = 4 ) nm
  real ( kind = 8 ) sk(0:n)
  real ( kind = 8 ) pi
  real ( kind = 8 ) x

  pi = 3.141592653589793D+00
  nm = n
  if ( x < 1.0D-60 ) then
     do k = 0,n
        sk(k) = 1.0D+300
        dk(k) = -1.0D+300
     end do
     return
  end if

  sk(0) = 0.5D+00 * pi / x * exp ( - x )
  sk(1) = sk(0) * ( 1.0D+00 + 1.0D+00 / x )
  f0 = sk(0)
  f1 = sk(1)
  do k = 2, n
     f = ( 2.0D+00 * k - 1.0D+00 ) * f1 / x + f0
     sk(k) = f
     if ( 1.0D+300 < abs ( f ) ) then
        exit
     end if
     f0 = f1
     f1 = f
  end do

  nm = k - 1

  dk(0) = -sk(1)
  do k = 1, nm
     dk(k) = -sk(k-1) - ( k + 1.0D+00 ) / x * sk(k)
  end do

  return
end subroutine sphk