itth0 Subroutine

subroutine itth0 ( x, tth )

  !*****************************************************************************80
  !
  !! ITTH0 integrates H0(t)/t from x to oo.
  !
  !  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:
  !
  !    23 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 ) X, the lower limit of the integral.
  !
  !    Output, real ( kind = 8 ) TTH, the integral of H0(t)/t from x to oo.
  !
  implicit none

  real ( kind = 8 ) f0
  real ( kind = 8 ) g0
  integer ( kind = 4 ) k
  real ( kind = 8 ) pi
  real ( kind = 8 ) r
  real ( kind = 8 ) s
  real ( kind = 8 ) t
  real ( kind = 8 ) tth
  real ( kind = 8 ) tty
  real ( kind = 8 ) x
  real ( kind = 8 ) xt

  pi = 3.141592653589793D+00
  s = 1.0D+00
  r = 1.0D+00

  if ( x < 24.5D+00 ) then

     do k = 1, 60
        r = - r * x * x * ( 2.0D+00 * k - 1.0D+00 ) &
             / ( 2.0D+00 * k + 1.0D+00 ) ** 3
        s = s + r
        if ( abs ( r ) < abs ( s ) * 1.0D-12 ) then
           exit
        end if
     end do

     tth = pi / 2.0D+00 - 2.0D+00 / pi * x * s

  else

     do k = 1, 10
        r = - r * ( 2.0D+00 * k - 1.0D+00 ) ** 3 &
             / ( ( 2.0D+00 * k + 1.0D+00 ) * x * x )
        s = s + r
        if ( abs ( r ) < abs ( s ) * 1.0D-12 ) then
           exit
        end if
     end do

     tth = 2.0D+00 / ( pi * x ) * s
     t = 8.0D+00 / x
     xt = x + 0.25D+00 * pi
     f0 = ((((( &
          0.18118D-02 * t &
          - 0.91909D-02 ) * t &
          + 0.017033D+00 ) * t &
          - 0.9394D-03 ) * t &
          - 0.051445D+00 ) * t &
          - 0.11D-05 ) * t &
          + 0.7978846D+00
     g0 = ((((( &
          - 0.23731D-02 * t &
          + 0.59842D-02 ) * t &
          + 0.24437D-02 ) * t &
          - 0.0233178D+00 ) * t &
          + 0.595D-04 ) * t &
          + 0.1620695D+00 ) * t
     tty = ( f0 * sin ( xt ) - g0 * cos ( xt ) ) / ( sqrt ( x ) * x )
     tth = tth + tty

  end if

  return
end subroutine itth0