stvhv Subroutine

subroutine stvhv ( v, x, hv )

  !*****************************************************************************80
  !
  !! STVHV computes the Struve function Hv(x) with arbitrary order v.
  !
  !  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:
  !
  !    24 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 ) V, the order of the function.
  !
  !    Input, real ( kind = 8 ) X, the argument.
  !
  !    Output, real ( kind = 8 ) HV, the value of Hv(x).
  !
  implicit none

  real ( kind = 8 ) bf
  real ( kind = 8 ) bf0
  real ( kind = 8 ) bf1
  real ( kind = 8 ) by0
  real ( kind = 8 ) by1
  real ( kind = 8 ) byv
  real ( kind = 8 ) ga
  real ( kind = 8 ) gb
  real ( kind = 8 ) hv
  integer ( kind = 4 ) k
  integer ( kind = 4 ) l
  integer ( kind = 4 ) n
  real ( kind = 8 ) pi
  real ( kind = 8 ) pu0
  real ( kind = 8 ) pu1
  real ( kind = 8 ) qu0
  real ( kind = 8 ) qu1
  real ( kind = 8 ) r1
  real ( kind = 8 ) r2
  real ( kind = 8 ) s
  real ( kind = 8 ) s0
  real ( kind = 8 ) sa
  real ( kind = 8 ) sr
  real ( kind = 8 ) t0
  real ( kind = 8 ) t1
  real ( kind = 8 ) u
  real ( kind = 8 ) u0
  real ( kind = 8 ) v
  real ( kind = 8 ) v0
  real ( kind = 8 ) va
  real ( kind = 8 ) vb
  real ( kind = 8 ) vt
  real ( kind = 8 ) x

  pi = 3.141592653589793D+00

  if ( x .eq. 0.0D+00 ) then
     if ( -1.0D+00 < v .or. int ( v ) - v == 0.5D+00 ) then
        hv = 0.0D+00
     else if ( v < -1.0D+00 ) then
        hv = ( -1.0D+00 ) ** ( int ( 0.5D+00 - v ) - 1 ) * 1.0D+300
     else if ( v == -1.0D+00 ) then
        hv = 2.0D+00 / pi
     end if
     return
  end if

  if ( x <= 20.0D+00 ) then

     v0 = v + 1.5D+00
     call gammaf ( v0, ga )
     s = 2.0D+00 / ( sqrt ( pi ) * ga )
     r1 = 1.0D+00

     do k = 1, 100
        va = k + 1.5D+00
        call gammaf ( va, ga )
        vb = v + k + 1.5D+00
        call gammaf ( vb, gb )
        r1 = -r1 * ( 0.5D+00 * x ) ** 2
        r2 = r1 / ( ga * gb )
        s = s + r2
        if ( abs ( r2 ) < abs ( s ) * 1.0D-12 ) then
           exit
        end if
     end do

     hv = ( 0.5D+00 * x ) ** ( v + 1.0D+00 ) * s

  else

     sa = ( 0.5D+00 * x ) ** ( v - 1.0D+00 ) / pi
     v0 = v + 0.5D+00
     call gammaf ( v0, ga )
     s = sqrt ( pi ) / ga
     r1 = 1.0D+00

     do k = 1, 12
        va = k + 0.5D+00
        call gammaf ( va, ga )
        vb = - k + v + 0.5D+00
        call gammaf ( vb, gb )
        r1 = r1 / ( 0.5D+00 * x ) ** 2
        s = s + r1 * ga / gb
     end do

     s0 = sa * s
     u = abs ( v )
     n = int ( u )
     u0 = u - n

     do l = 0, 1

        vt = 4.0D+00 * ( u0 + l ) ** 2
        r1 = 1.0D+00
        pu1 = 1.0D+00
        do k = 1, 12
           r1 = -0.0078125D+00 * r1 &
                * ( vt - ( 4.0D+00 * k - 3.0D+00 ) ** 2 ) &
                * ( vt - ( 4.0D+00 * k - 1.0D+00 ) ** 2 ) &
                / ( ( 2.0D+00 * k - 1.0D+00 ) * k * x * x )
           pu1 = pu1 + r1
        end do

        qu1 = 1.0D+00
        r2 = 1.0D+00
        do k = 1, 12
           r2 = -0.0078125D+00 * r2 &
                * ( vt - ( 4.0D+00 * k - 1.0D+00 ) ** 2 ) &
                * ( vt - ( 4.0D+00 * k + 1.0D+00 ) ** 2 ) &
                / ( ( 2.0D+00 * k + 1.0D+00 ) * k * x * x )
           qu1 = qu1 + r2
        end do
        qu1 = 0.125D+00 * ( vt - 1.0D+00 ) / x * qu1

        if ( l == 0 ) then
           pu0 = pu1
           qu0 = qu1
        end if

     end do

     t0 = x - ( 0.5D+00 * u0 + 0.25D+00 ) * pi
     t1 = x - ( 0.5D+00 * u0 + 0.75D+00 ) * pi
     sr = sqrt ( 2.0D+00 / ( pi * x ) )
     by0 = sr * ( pu0 * sin ( t0 ) + qu0 * cos ( t0 ) )
     by1 = sr * ( pu1 * sin ( t1 ) + qu1 * cos ( t1 ) )
     bf0 = by0
     bf1 = by1
     do k = 2, n
        bf = 2.0D+00 * ( k - 1.0D+00 + u0 ) / x * bf1 - bf0
        bf0 = bf1
        bf1 = bf
     end do

     if ( n == 0 ) then
        byv = by0
     else if ( n == 1 ) then
        byv = by1
     else
        byv = bf
     end if
     hv = byv + s0
  end if

  return
end subroutine stvhv