rcty Subroutine

subroutine rcty ( n, x, nm, ry, dy )

  !*****************************************************************************80
  !
  !! RCTY computes Riccati-Bessel function of the second kind, and 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:
  !
  !    18 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 yn(x).
  !
  !    Input, real ( kind = 8 ) X, the argument.
  !
  !    Output, integer ( kind = 4 ) NM, the highest order computed.
  !
  !    Output, real ( kind = 8 ) RY(0:N), the values of x yn(x).
  !
  !    Output, real ( kind = 8 ) DY(0:N), the values of [x yn(x)]'.
  !
  implicit none

  integer ( kind = 4 ) n

  real ( kind = 8 ) dy(0:n)
  integer ( kind = 4 ) k
  integer ( kind = 4 ) nm
  real ( kind = 8 ) rf0
  real ( kind = 8 ) rf1
  real ( kind = 8 ) rf2
  real ( kind = 8 ) ry(0:n)
  real ( kind = 8 ) x

  nm = n

  if ( x < 1.0D-60 ) then
     do k = 0, n
        ry(k) = -1.0D+300
        dy(k) = 1.0D+300
     end do
     ry(0) = -1.0D+00
     dy(0) = 0.0D+00
     return
  end if

  ry(0) = - cos ( x )
  ry(1) = ry(0) / x - sin ( x )
  rf0 = ry(0)
  rf1 = ry(1)
  do k = 2, n
     rf2 = ( 2.0D+00 * k - 1.0D+00 ) * rf1 / x - rf0
     if ( 1.0D+300 < abs ( rf2 ) ) then
        exit
     end if
     ry(k) = rf2
     rf0 = rf1
     rf1 = rf2
  end do

  nm = k - 1
  dy(0) = sin ( x )
  do k = 1, nm
     dy(k) = - k * ry(k) / x + ry(k-1)
  end do

  return
end subroutine rcty