chgus Subroutine

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

  !*****************************************************************************80
  !
  !! CHGUS: confluent hypergeometric function U(a,b,x) for small 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:
  !
  !    27 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, U(a,b,x).
  !
  !    Output, integer ( kind = 4 ) ID, the estimated number of 
  !    significant digits.
  !
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) b
  real ( kind = 8 ) d1
  real ( kind = 8 ) d2
  real ( kind = 8 ) ga
  real ( kind = 8 ) gab
  real ( kind = 8 ) gb
  real ( kind = 8 ) gb2
  real ( kind = 8 ) h0
  real ( kind = 8 ) hmax
  real ( kind = 8 ) hmin
  real ( kind = 8 ) hu
  real ( kind = 8 ) hu0
  real ( kind = 8 ) hua
  integer ( kind = 4 ) id
  integer ( kind = 4 ) j
  real ( kind = 8 ) pi
  real ( kind = 8 ) r1
  real ( kind = 8 ) r2
  real ( kind = 8 ) x
  real ( kind = 8 ) xg1
  real ( kind = 8 ) xg2

  id = -100
  pi = 3.141592653589793D+00
  call gammaf ( a, ga )
  call gammaf ( b, gb )
  xg1 = 1.0D+00 + a - b
  call gammaf ( xg1, gab )
  xg2 = 2.0D+00 - b
  call gammaf ( xg2, gb2 )
  hu0 = pi / sin ( pi * b )
  r1 = hu0 / ( gab * gb )
  r2 = hu0 * x ** ( 1.0D+00 - b ) / ( ga * gb2 )
  hu = r1 - r2
  hmax = 0.0D+00
  hmin = 1.0D+300
  do j = 1, 150
     r1 = r1 * ( a + j - 1.0D+00 ) / ( j * ( b + j - 1.0D+00 ) ) * x
     r2 = r2 * ( a - b + j ) / ( j * ( 1.0D+00 - b + j ) ) * x
     hu = hu + r1 - r2
     hua = abs ( hu )
     hmax = max ( hmax, hua )
     hmin = min ( hmin, hua )
     if ( abs ( hu - h0 ) < abs ( hu ) * 1.0D-15 ) then
        exit
     end if
     h0 = hu
  end do

  d1 = log10 ( hmax )
  if ( hmin /= 0.0D+00 ) then
     d2 = log10 ( hmin )
  end if
  id = 15 - abs ( d1 - d2 )

  return
end subroutine chgus