stvlv Subroutine

subroutine stvlv ( v, x, slv )

  !*****************************************************************************80
  !
  !! STVLV computes the modified Struve function Lv(x) with arbitary order.
  !
  !  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:
  !
  !    04 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 Lv(x).
  !
  !    Input, real ( kind = 8 ) X, the argument of Lv(x).
  !
  !    Output, real ( kind = 8 ) SLV, the value of Lv(x).
  !
  implicit none

  real ( kind = 8 ) bf
  real ( kind = 8 ) bf0
  real ( kind = 8 ) bf1
  real ( kind = 8 ) biv
  real ( kind = 8 ) biv0
  real ( kind = 8 ) ga
  real ( kind = 8 ) gb
  integer ( kind = 4 ) k
  integer ( kind = 4 ) l
  integer ( kind = 4 ) n
  real ( kind = 8 ) pi
  real ( kind = 8 ) r
  real ( kind = 8 ) r1
  real ( kind = 8 ) r2
  real ( kind = 8 ) s
  real ( kind = 8 ) s0
  real ( kind = 8 ) sa
  real ( kind = 8 ) slv
  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 == 0.0D+00 ) then

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

  else if ( x <= 40.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 / s ) < 1.0D-12 ) then
           exit
        end if
     end do

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

  else

     sa = -1.0D+00 / pi * ( 0.5D+00 * x ) ** ( v - 1.0D+00 )
     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 = u0 + l
        r = 1.0D+00
        biv = 1.0D+00
        do k = 1, 16
           r = -0.125D+00 * r * ( 4.0D+00 * vt * vt - &
                ( 2.0D+00 * k - 1.0D+00 )**2 ) / ( k * x )
           biv = biv + r
           if ( abs ( r / biv ) < 1.0D-12 ) then
              exit
           end if
        end do

        if ( l == 0 ) then
           biv0 = biv
        end if

     end do

     bf0 = biv0
     bf1 = biv
     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
        biv = biv0
     else if ( 1 < n ) then
        biv = bf
     end if

     slv = exp ( x ) / sqrt ( 2.0D+00 * pi * x ) * biv + s0

  end if

  return
end subroutine stvlv