chgul Subroutine

subroutine chgul ( a, b, x, hu, id )

  !*****************************************************************************80
  !
  !! CHGUL: confluent hypergeometric function U(a,b,x) for large argument 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:
  !
  !    22 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, real ( kind = 8 ) A, B, parameters.
  !
  !    Input, real ( kind = 8 ) X, the argument.
  !
  !    Output, real ( kind = 8 ) HU, the value of U(a,b,x).
  !
  !    Output, integer ( kind = 4 ) ID, the estimated number of 
  !    significant digits.
  !
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) aa
  real ( kind = 8 ) b
  real ( kind = 8 ) hu
  integer ( kind = 4 ) id
  logical il1
  logical il2
  integer ( kind = 4 ) k
  integer ( kind = 4 ) nm
  real ( kind = 8 ) r
  real ( kind = 8 ) ra
  real ( kind = 8 ) r0
  real ( kind = 8 ) x

  id = -100
  aa = a - b + 1.0D+00
  il1 = ( a == int ( a ) ) .and. ( a <= 0.0D+00 )
  il2 = ( aa == int ( aa ) ) .and. ( aa <= 0.0D+00 )

  if ( il1 .or. il2 ) then

     if ( il1 ) then
        nm = abs ( a )
     end if

     if ( il2 ) then
        nm = abs ( aa )
     end if

     hu = 1.0D+00
     r = 1.0D+00
     do k = 1, nm
        r = - r * ( a + k - 1.0D+00 ) * ( a - b + k ) / ( k * x )
        hu = hu + r
     end do
     hu = x ** ( - a ) * hu
     id = 10

  else

     hu = 1.0D+00
     r = 1.0D+00
     do k = 1, 25
        r = - r * ( a + k - 1.0D+00 ) * ( a - b + k ) / ( k * x )
        ra = abs ( r )
        if ( ( 5 < k .and. r0 <= ra ) .or. ra < 1.0D-15 ) then
           exit
        end if
        r0 = ra
        hu = hu + r
     end do

     id = abs ( log10 ( ra ) )
     hu = x ** ( - a ) * hu

  end if

  return
end subroutine chgul