vvsa Subroutine

subroutine vvsa ( va, x, pv )

  !*****************************************************************************80
  !
  !! VVSA computes parabolic cylinder function V(nu,x) for small arguments.
  !
  !  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 ) X, the argument.
  !
  !    Input, real ( kind = 8 ) VA, the order nu.
  !
  !    Output, real ( kind = 8 ) PV, the value of V(nu,x).
  !
  implicit none

  real ( kind = 8 ) a0
  real ( kind = 8 ) ep
  real ( kind = 8 ) eps
  real ( kind = 8 ) fac
  real ( kind = 8 ) g1
  real ( kind = 8 ) ga0
  real ( kind = 8 ) gm
  real ( kind = 8 ) gw
  integer ( kind = 4 ) m
  real ( kind = 8 ) pi
  real ( kind = 8 ) pv
  real ( kind = 8 ) r
  real ( kind = 8 ) r1
  real ( kind = 8 ) sq2
  real ( kind = 8 ) sv
  real ( kind = 8 ) sv0
  real ( kind = 8 ) v1
  real ( kind = 8 ) va
  real ( kind = 8 ) va0
  real ( kind = 8 ) vb0
  real ( kind = 8 ) vm
  real ( kind = 8 ) x

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

  if ( x == 0.0D+00 ) then

     if ( ( va0 <= 0.0D+00 .and. va0 == int ( va0 ) ) .or. &
          va == 0.0D+00 ) then
        pv = 0.0D+00
     else
        vb0 = -0.5D+00 * va
        sv0 = sin ( va0 * pi )
        call gammaf ( va0, ga0 )
        pv = 2.0D+00 ** vb0 * sv0 / ga0
     end if

  else

     sq2 = sqrt ( 2.0D+00 )
     a0 = 2.0D+00 ** ( -0.5D+00 * va ) * ep / ( 2.0D+00 * pi )
     sv = sin ( - ( va + 0.5D+00 ) * pi )
     v1 = -0.5D+00 * va
     call gammaf ( v1, g1 )
     pv = ( sv + 1.0D+00 ) * g1
     r = 1.0D+00
     fac = 1.0D+00

     do m = 1, 250
        vm = 0.5D+00 * ( m - va )
        call gammaf ( vm, gm )
        r = r * sq2 * x / m
        fac = - fac
        gw = fac * sv + 1.0D+00
        r1 = gw * r * gm
        pv = pv + r1
        if ( abs ( r1 / pv ) < eps .and. gw /= 0.0D+00 ) then
           exit
        end if
     end do

     pv = a0 * pv

  end if

  return
end subroutine vvsa