itsh0 Subroutine

subroutine itsh0 ( x, th0 )

  !*****************************************************************************80
  !
  !! ITSH0 integrates the Struve function H0(t) from 0 to x.
  !
  !  Discussion:
  !
  !    This procedure evaluates the integral of Struve function
  !    H0(t) with respect to t from 0 and 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:
  !
  !    25 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 upper limit of the integral.
  !
  !    Output, real ( kind = 8 ) TH0, the integral of H0(t) from 0 to x.
  !
  implicit none

  real ( kind = 8 ) a(25)
  real ( kind = 8 ) a0
  real ( kind = 8 ) a1
  real ( kind = 8 ) af
  real ( kind = 8 ) bf
  real ( kind = 8 ) bg
  real ( kind = 8 ) el
  integer ( kind = 4 ) k
  real ( kind = 8 ) pi
  real ( kind = 8 ) r
  real ( kind = 8 ) rd
  real ( kind = 8 ) s
  real ( kind = 8 ) s0
  real ( kind = 8 ) th0
  real ( kind = 8 ) ty
  real ( kind = 8 ) x
  real ( kind = 8 ) xp

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

  if ( x <= 30.0D+00 ) then

     s = 0.5D+00

     do k = 1, 100

        if ( k == 1 ) then
           rd = 0.5D+00
        else
           rd = 1.0D+00
        end if

        r = - r * rd * k / ( k + 1.0D+00 ) &
             * ( x / ( 2.0D+00 * k + 1.0D+00 ) ) ** 2
        s = s + r

        if ( abs ( r ) < abs ( s ) * 1.0D-12 ) then
           exit
        end if

     end do

     th0 = 2.0D+00 / pi * x * x * s

  else

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

     el = 0.57721566490153D+00
     s0 = s / ( pi * x * x ) + 2.0D+00 / pi &
          * ( log ( 2.0D+00 * x ) + el )
     a0 = 1.0D+00
     a1 = 5.0D+00 / 8.0D+00
     a(1) = a1
     do k = 1, 20
        af = ( ( 1.5D+00 * ( k + 0.5D+00 ) &
             * ( k + 5.0D+00 / 6.0D+00 ) * a1 - 0.5D+00 &
             * ( k + 0.5D+00 ) * ( k + 0.5D+00 ) &
             * ( k - 0.5D+00 ) * a0 ) ) / ( k + 1.0D+00 )
        a(k+1) = af
        a0 = a1
        a1 = af
     end do

     bf = 1.0D+00
     r = 1.0D+00
     do k = 1, 10
        r = - r / ( x * x )
        bf = bf + a(2*k) * r
     end do
     bg = a(1) / x
     r = 1.0D+00 / x
     do k = 1, 10
        r = - r / ( x * x ) 
        bg = bg + a(2*k+1) * r
     end do
     xp = x + 0.25D+00 * pi
     ty = sqrt ( 2.0D+00 / ( pi * x ) ) &
          * ( bg * cos ( xp ) - bf * sin ( xp ) )
     th0 = ty + s0

  end if

  return
end subroutine itsh0