chguit Subroutine

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

  !*****************************************************************************80
  !
  !! CHGUIT computes the hypergeometric function using Gauss-Legendre integration.
  !
  !  Discussion:
  !
  !    This procedure computes the hypergeometric function U(a,b,x) by
  !    using Gaussian-Legendre integration (n = 60)
  !
  !  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, double precision A, B, parameters.
  !
  !    Input, double precision X, the argument.
  !
  !    Output, double precision HU, U(a,b,z).
  !
  !    Output, integer ID, the estimated number of significant digits.
  !
  implicit none

  double precision a
  double precision a1
  double precision b
  double precision b1
  double precision c
  double precision d
  double precision f1
  double precision f2
  double precision g
  double precision ga
  double precision hu
  double precision hu0
  double precision hu1
  double precision hu2
  integer id
  integer j
  integer k
  integer m
  double precision s
  double precision, save, dimension ( 30 ) :: t = (/ &
       0.259597723012478D-01, 0.778093339495366D-01, &
       0.129449135396945D+00, 0.180739964873425D+00, &
       0.231543551376029D+00, 0.281722937423262D+00, &
       0.331142848268448D+00, 0.379670056576798D+00, &
       0.427173741583078D+00, 0.473525841761707D+00, &
       0.518601400058570D+00, 0.562278900753945D+00, &
       0.604440597048510D+00, 0.644972828489477D+00, &
       0.683766327381356D+00, 0.720716513355730D+00, &
       0.755723775306586D+00, 0.788693739932264D+00, &
       0.819537526162146D+00, 0.848171984785930D+00, &
       0.874519922646898D+00, 0.898510310810046D+00, &
       0.920078476177628D+00, 0.939166276116423D+00, &
       0.955722255839996D+00, 0.969701788765053D+00, &
       0.981067201752598D+00, 0.989787895222222D+00, &
       0.995840525118838D+00, 0.999210123227436D+00 /)
  double precision t1
  double precision t2
  double precision t3
  double precision t4
  double precision, save, dimension ( 30 ) :: w = (/ &
       0.519078776312206D-01, 0.517679431749102D-01, &
       0.514884515009810D-01, 0.510701560698557D-01, &
       0.505141845325094D-01, 0.498220356905502D-01, &
       0.489955754557568D-01, 0.480370318199712D-01, &
       0.469489888489122D-01, 0.457343797161145D-01, &
       0.443964787957872D-01, 0.429388928359356D-01, &
       0.413655512355848D-01, 0.396806954523808D-01, &
       0.378888675692434D-01, 0.359948980510845D-01, &
       0.340038927249464D-01, 0.319212190192963D-01, &
       0.297524915007890D-01, 0.275035567499248D-01, &
       0.251804776215213D-01, 0.227895169439978D-01, &
       0.203371207294572D-01, 0.178299010142074D-01, &
       0.152746185967848D-01, 0.126781664768159D-01, &
       0.100475571822880D-01, 0.738993116334531D-02, &
       0.471272992695363D-02, 0.202681196887362D-02 /)
  double precision x

  id = 7
  a1 = a - 1.0D+00
  b1 = b - a - 1.0D+00
  c = 12.0D+00 / x

  do m = 10, 100, 5

     hu1 = 0.0D+00
     g = 0.5D+00 * c / m
     d = g
     do j = 1, m
        s = 0.0D+00
        do k = 1, 30
           t1 = d + g * t(k)
           t2 = d - g * t(k)
           f1 = exp ( - x * t1 ) * t1 ** a1 * ( 1.0D+00 + t1 ) ** b1
           f2 = exp ( - x * t2 ) * t2 ** a1 * ( 1.0D+00 + t2 ) ** b1
           s = s + w(k) * ( f1 + f2 )
        end do
        hu1 = hu1 + s * g
        d = d + 2.0D+00 * g
     end do

     if ( abs ( 1.0D+00 - hu0 / hu1 ) < 1.0D-07 ) then
        exit
     end if

     hu0 = hu1

  end do

  call gammaf ( a, ga )
  hu1 = hu1 / ga

  do m = 2, 10, 2
     hu2 = 0.0D+00
     g = 0.5D+00 / m
     d = g
     do j = 1, m
        s = 0.0D+00
        do k = 1, 30
           t1 = d + g * t(k)
           t2 = d - g * t(k)
           t3 = c / ( 1.0D+00 - t1 )
           t4 = c / ( 1.0D+00 - t2 ) 
           f1 = t3 * t3 / c * exp ( - x * t3 ) * t3 ** a1 * ( 1.0D+00 + t3 ) ** b1
           f2 = t4 * t4 / c * exp ( - x * t4 ) * t4 ** a1 * ( 1.0D+00 + t4 ) ** b1
           s = s + w(k) * ( f1 + f2 )
        end do
        hu2 = hu2 + s * g
        d = d + 2.0D+00 * g
     end do

     if ( abs ( 1.0D+00 - hu0 / hu2 ) < 1.0D-07 ) then
        exit
     end if

     hu0 = hu2

  end do

  call gammaf ( a, ga )
  hu2 = hu2 / ga
  hu = hu1 + hu2

  return
end subroutine chguit