stvl0 Subroutine

subroutine stvl0(x, sl0)

************80

! STVL0 computes the modified Struve function L0(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:

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, real ( kind = 8 ) X, the argument.

Output, real ( kind = 8 ) SL0, the function value.

Arguments

Type IntentOptional Attributes Name
real(kind=8) :: x
real(kind=8) :: sl0

Source Code

subroutine stvl0 ( x, sl0 )

  !*****************************************************************************80
  !
  !! STVL0 computes the modified Struve function L0(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:
  !
  !    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, real ( kind = 8 ) X, the argument.
  !
  !    Output, real ( kind = 8 ) SL0, the function value.
  !
  implicit none

  real ( kind = 8 ) a0
  real ( kind = 8 ) a1
  real ( kind = 8 ) bi0
  integer ( kind = 4 ) k
  integer ( kind = 4 ) km
  real ( kind = 8 ) pi
  real ( kind = 8 ) r
  real ( kind = 8 ) s
  real ( kind = 8 ) sl0
  real ( kind = 8 ) x

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

  if ( x <= 20.0D+00 ) then

     a0 = 2.0D+00 * x / pi

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

     sl0 = a0 * s

  else

     if ( x < 50.0D+00 ) then
        km = int ( 0.5D+00 * ( x + 1.0D+00 ) )
     else
        km = 25
     end if

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

     a1 = exp ( x ) / sqrt ( 2.0D+00 * pi * x )
     r = 1.0D+00
     bi0 = 1.0D+00
     do k = 1, 16
        r = 0.125D+00 * r * ( 2.0D+00 * k - 1.0D+00 ) ** 2 / ( k * x )
        bi0 = bi0 + r
        if ( abs ( r / bi0 ) < 1.0D-12 ) then
           exit
        end if
     end do

     bi0 = a1 * bi0
     sl0 = - 2.0D+00 / ( pi * x ) * s + bi0

  end if

  return
end subroutine stvl0