pbdv Subroutine

subroutine pbdv ( v, x, dv, dp, pdf, pdd )

  !*****************************************************************************80
  !
  !! PBDV computes parabolic cylinder functions Dv(x) and 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 ) DV(0:*), DP(0:*), the values of
  !    Dn+v0(x), Dn+v0'(x).
  !
  !    Output, real ( kind = 8 ) PDF, PDD, the values of Dv(x) and Dv'(x).
  !
  implicit none

  real ( kind = 8 ) dp(0:*)
  real ( kind = 8 ) dv(0:*)
  real ( kind = 8 ) ep
  real ( kind = 8 ) f
  real ( kind = 8 ) f0
  real ( kind = 8 ) f1
  integer ( kind = 4 ) ja
  integer ( kind = 4 ) k
  integer ( kind = 4 ) l
  integer ( kind = 4 ) m
  integer ( kind = 4 ) na
  integer ( kind = 4 ) nk
  integer ( kind = 4 ) nv
  real ( kind = 8 ) pd
  real ( kind = 8 ) pd0
  real ( kind = 8 ) pd1
  real ( kind = 8 ) pdd
  real ( kind = 8 ) pdf
  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 ) x
  real ( kind = 8 ) xa

  xa = abs ( x )
  vh = v
  v = v + sign ( 1.0D+00, v )
  nv = int ( v )
  v0 = v - nv
  na = abs ( nv )
  ep = exp ( -0.25D+00 * x * x )

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

  if ( 0.0D+00 <= v ) then
     if ( v0 == 0.0D+00 ) then
        pd0 = ep
        pd1 = x * ep
     else
        do l = 0, ja
           v1 = v0 + l
           if ( xa <= 5.8D+00 ) then
              call dvsa ( v1, x, pd1 )
           else
              call dvla ( v1, x, pd1 )
           end if
           if ( l == 0 ) then
              pd0 = pd1
           end if
        end do
     end if

     dv(0) = pd0
     dv(1) = pd1
     do k = 2, na
        pdf = x * pd1 - ( k + v0 - 1.0D+00 ) * pd0
        dv(k) = pdf
        pd0 = pd1
        pd1 = pdf
     end do

  else

     if ( x <= 0.0D+00 ) then

        if ( xa <= 5.8D+00 )  then
           call dvsa ( v0, x, pd0 )
           v1 = v0 - 1.0D+00
           call dvsa ( v1, x, pd1 )
        else
           call dvla ( v0, x, pd0 )
           v1 = v0 - 1.0D+00
           call dvla ( v1, x, pd1 )
        end if

        dv(0) = pd0
        dv(1) = pd1
        do k = 2, na
           pd = ( - x * pd1 + pd0 ) / ( k - 1.0D+00 - v0 )
           dv(k) = pd
           pd0 = pd1
           pd1 = pd
        end do

     else if ( x <= 2.0D+00 ) then

        v2 = nv + v0
        if ( nv == 0 ) then
           v2 = v2 - 1.0D+00
        end if

        nk = int ( - v2 )
        call dvsa ( v2, x, f1 )
        v1 = v2 + 1.0D+00
        call dvsa ( v1, x, f0 )
        dv(nk) = f1
        dv(nk-1) = f0
        do k = nk - 2, 0, -1
           f = x * f0 + ( k - v0 + 1.0D+00 ) * f1
           dv(k) = f
           f1 = f0
           f0 = f
        end do

     else

        if ( xa <= 5.8D+00 ) then
           call dvsa ( v0, x, pd0 )
        else
           call dvla ( v0, x, pd0 )
        end if

        dv(0) = pd0
        m = 100 + na
        f1 = 0.0D+00
        f0 = 1.0D-30
        do k = m, 0, -1
           f = x * f0 + ( k - v0 + 1.0D+00 ) * f1
           if ( k <= na ) then
              dv(k) = f
           end if
           f1 = f0
           f0 = f
        end do
        s0 = pd0 / f
        do k = 0, na
           dv(k) = s0 * dv(k)
        end do

     end if

  end if

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

  pdf = dv(na-1)
  pdd = dp(na-1)
  v = vh

  return
end subroutine pbdv