pbvv Subroutine

subroutine pbvv ( v, x, vv, vp, pvf, pvd )

  !*****************************************************************************80
  !
  !! PBVV computes parabolic cylinder functions Vv(x) and their derivatives.
  !
  !  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:
  !
  !    29 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.
  !
  !    Input, real ( kind = 8 ) X, the argument.
  !
  !    Output, real ( kind = 8 ) VV(0:*), VP(0:*), the values of Vv(x), Vv'(x).
  !
  !    Output, real ( kind = 8 ) PVF, PVD, the values of Vv(x) and Vv'(x).
  !
  implicit none

  real ( kind = 8 ) f
  real ( kind = 8 ) f0
  real ( kind = 8 ) f1
  integer ( kind = 4 ) ja
  integer ( kind = 4 ) k
  integer ( kind = 4 ) kv
  integer ( kind = 4 ) l
  integer ( kind = 4 ) m
  integer ( kind = 4 ) na
  integer ( kind = 4 ) nv
  real ( kind = 8 ) pi
  real ( kind = 8 ) pv0
  real ( kind = 8 ) pvd
  real ( kind = 8 ) pvf
  real ( kind = 8 ) q2p
  real ( kind = 8 ) qe
  real ( kind = 8 ) s0
  real ( kind = 8 ) v
  real ( kind = 8 ) v0
  real ( kind = 8 ) v1
  real ( kind = 8 ) v2
  real ( kind = 8 ) vh
  real ( kind = 8 ) vp(0:*)
  real ( kind = 8 ) vv(0:*)
  real ( kind = 8 ) x
  real ( kind = 8 ) xa

  pi = 3.141592653589793D+00
  xa = abs ( x )
  vh = v
  v = v + sign ( 1.0D+00, v )
  nv = int ( v )
  v0 = v - nv
  na = abs ( nv )
  qe = exp ( 0.25D+00 * x * x )
  q2p = sqrt ( 2.0D+00 / pi )

  if ( 1 <= na ) then
     ja = 1
  end if

  if ( v <= 0.0D+00 ) then

     if ( v0 == 0.0D+00 ) then

        if ( xa <= 7.5D+00 ) then 
           call vvsa ( v0, x, pv0 )
        else
           call vvla ( v0, x, pv0 )
        end if

        f0 = q2p * qe
        f1 = x * f0
        vv(0) = pv0
        vv(1) = f0
        vv(2) = f1

     else

        do l = 0, ja
           v1 = v0 - l
           if ( xa <= 7.5D+00 ) then
              call vvsa ( v1, x, f1 )
           else
              call vvla ( v1, x, f1 )
           end if
           if ( l == 0 ) then
              f0 = f1
           end if
        end do

        vv(0) = f0
        vv(1) = f1

     end if

     if ( v0 == 0.0D+00 ) then
        kv = 3
     else
        kv = 2
     end if

     do k = kv, na
        f = x * f1 + ( k - v0 - 2.0D+00 ) * f0
        vv(k) = f
        f0 = f1
        f1 = f
     end do

  else

     if ( 0.0D+00 <= x .and. x <= 7.5D+00 ) then

        v2 = v
        if ( v2 < 1.0D+00 ) then
           v2 = v2 + 1.0D+00
        end if

        call vvsa ( v2, x, f1 )
        v1 = v2 - 1.0D+00
        kv = int ( v2 )
        call vvsa ( v1, x, f0 )
        vv(kv) = f1
        vv(kv-1) = f0
        do k = kv - 2, 0, - 1
           f = x * f0 - ( k + v0 + 2.0D+00 ) * f1
           if ( k <= na ) then
              vv(k) = f
           end if
           f1 = f0
           f0 = f
        end do

     else if ( 7.5D+00 < x ) then

        call vvla ( v0, x, pv0 )
        m = 100 + abs ( na )
        vv(1) = pv0
        f1 = 0.0D+00
        f0 = 1.0D-40
        do k = m, 0, -1
           f = x * f0 - ( k + v0 + 2.0D+00 ) * f1
           if ( k <= na ) then
              vv(k) = f
           end if
           f1 = f0
           f0 = f
        end do
        s0 = pv0 / f
        do k = 0, na
           vv(k) = s0 * vv(k)
        end do

     else

        if ( xa <= 7.5D+00 ) then
           call vvsa ( v0, x, f0 )
           v1 = v0 + 1.0D+00
           call vvsa ( v1, x, f1 )
        else
           call vvla ( v0, x, f0 )
           v1 = v0 + 1.0D+00
           call vvla ( v1, x, f1 )
        end if

        vv(0) = f0
        vv(1) = f1
        do k = 2, na
           f = ( x * f1 - f0 ) / ( k + v0 )
           vv(k) = f
           f0 = f1
           f1 = f
        end do

     end if

  end if

  do k = 0, na - 1
     v1 = v0 + k
     if ( 0.0D+00 <= v ) then
        vp(k) = 0.5D+00 * x * vv(k) - ( v1 + 1.0D+00 ) * vv(k+1)
     else
        vp(k) = - 0.5D+00 * x * vv(k) + vv(k+1)
     end if
  end do

  pvf = vv(na-1)
  pvd = vp(na-1)
  v = vh

  return
end subroutine pbvv