dvsa Subroutine

subroutine dvsa ( va, x, pd )

  !*****************************************************************************80
  !
  !! DVSA computes parabolic cylinder functions Dv(x) for small argument.
  !
  !  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:
  !
  !    07 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 ) VA, the order.
  !
  !    Input, real ( kind = 8 ) X, the argument.
  !
  !    Output, real ( kind = 8 ) PD, the function value.
  !
  implicit none

  real ( kind = 8 ) a0
  real ( kind = 8 ) ep
  real ( kind = 8 ) eps
  real ( kind = 8 ) g0
  real ( kind = 8 ) g1
  real ( kind = 8 ) ga0
  real ( kind = 8 ) gm
  integer ( kind = 4 ) m
  real ( kind = 8 ) pd
  real ( kind = 8 ) pi
  real ( kind = 8 ) r
  real ( kind = 8 ) r1
  real ( kind = 8 ) sq2
  real ( kind = 8 ) va
  real ( kind = 8 ) va0
  real ( kind = 8 ) vm
  real ( kind = 8 ) vt
  real ( kind = 8 ) x

  eps = 1.0D-15
  pi = 3.141592653589793D+00
  sq2 = sqrt ( 2.0D+00 )
  ep = exp ( -0.25D+00 * x * x )
  va0 = 0.5D+00 * ( 1.0D+00 - va )

  if ( va == 0.0D+00 ) then

     pd = ep

  else

     if ( x == 0.0D+00 ) then
        if ( va0 <= 0.0D+00 .and. va0 == int ( va0 ) ) then
           pd = 0.0D+00
        else
           call gammaf ( va0, ga0 )
           pd = sqrt ( pi ) / ( 2.0D+00 ** ( -0.5D+00 * va ) * ga0 )
        end if

     else

        call gammaf ( -va, g1 )
        a0 = 2.0D+00 ** ( -0.5D+00 * va - 1.0D+00 ) * ep / g1
        vt = -0.5D+00 * va
        call gammaf ( vt, g0 )
        pd = g0
        r = 1.0D+00
        do m = 1, 250
           vm = 0.5D+00 * ( m - va )
           call gammaf ( vm, gm )
           r = -r * sq2 * x / m
           r1 = gm * r
           pd = pd + r1
           if ( abs ( r1 ) < abs ( pd ) * eps ) then
              exit
           end if
        end do

        pd = a0 * pd

     end if

  end if

  return
end subroutine dvsa