vvla Subroutine

subroutine vvla ( va, x, pv )

  !*****************************************************************************80
  !
  !! VVLA computes parabolic cylinder function Vv(x) for large 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 ) dsl
  real ( kind = 8 ) eps
  real ( kind = 8 ) gl
  integer ( kind = 4 ) k
  real ( kind = 8 ) pdl
  real ( kind = 8 ) pi
  real ( kind = 8 ) pv
  real ( kind = 8 ) qe
  real ( kind = 8 ) r
  real ( kind = 8 ) va
  real ( kind = 8 ) x
  real ( kind = 8 ) x1

  pi = 3.141592653589793D+00
  eps = 1.0D-12
  qe = exp ( 0.25D+00 * x * x )
  a0 = abs ( x ) ** ( -va - 1.0D+00 ) * sqrt ( 2.0D+00 / pi ) * qe

  r = 1.0D+00
  pv = 1.0D+00
  do k = 1, 18
     r = 0.5D+00 * r * ( 2.0D+00 * k + va - 1.0D+00 ) &
          * ( 2.0D+00 * k + va ) / ( k * x * x )
     pv = pv + r
     if ( abs ( r / pv ) < eps ) then
        exit
     end if
  end do

  pv = a0 * pv

  if ( x < 0.0D+00 ) then
     x1 = -x
     call dvla ( va, x1, pdl )
     call gammaf ( -va, gl )
     dsl = sin ( pi * va ) * sin ( pi * va )
     pv = dsl * gl / pi * pdl - cos ( pi * va ) * pv
  end if

  return
end subroutine vvla