cpbdn Subroutine

subroutine cpbdn ( n, z, cpb, cpd )

  !*****************************************************************************80
  !
  !! CPBDN: parabolic cylinder function Dn(z) and Dn'(z) for complex 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:
  !
  !    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, integer ( kind = 4 ) N, the order.
  !
  !    Input, complex ( kind = 8 ) Z, the argument.
  !
  !    Output, complex ( kind = 8 ) CPB(0:N), CPD(0:N), the values of Dn(z) 
  !    and Dn'(z).
  !
  implicit none

  integer ( kind = 4 ) n

  real ( kind = 8 ) a0
  complex ( kind = 8 ) c0
  complex ( kind = 8 ) ca0
  complex ( kind = 8 ) cf
  complex ( kind = 8 ) cf0
  complex ( kind = 8 ) cf1
  complex ( kind = 8 ) cfa
  complex ( kind = 8 ) cfb
  complex ( kind = 8 ) cpb(0:n)
  complex ( kind = 8 ) cpd(0:n)
  complex ( kind = 8 ) cs0
  integer ( kind = 4 ) k
  integer ( kind = 4 ) m
  integer ( kind = 4 ) n0
  integer ( kind = 4 ) n1
  integer ( kind = 4 ) nm1
  real ( kind = 8 ) pi
  real ( kind = 8 ) x
  complex ( kind = 8 ) z
  complex ( kind = 8 ) z1

  pi = 3.141592653589793D+00
  x = real ( z, kind = 8 )
  a0 = abs ( z )
  c0 = cmplx ( 0.0D+00, 0.0D+00, kind = 8 )
  ca0 = exp ( -0.25D+00 * z * z )

  if ( 0 <= n ) then

     cf0 = ca0
     cf1 = z * ca0
     cpb(0) = cf0
     cpb(1) = cf1
     do k = 2, n
        cf = z * cf1 - ( k - 1.0D+00 ) * cf0
        cpb(k) = cf
        cf0 = cf1
        cf1 = cf
     end do

  else

     n0 = -n

     if ( x <= 0.0D+00 .or. abs ( z ) .eq. 0.0D+00 ) then

        cf0 = ca0
        cpb(0) = cf0
        z1 = - z
        if ( a0 <= 7.0D+00 ) then
           call cpdsa ( -1, z1, cf1 )
        else
           call cpdla ( -1, z1, cf1 )
        end if
        cf1 = sqrt ( 2.0D+00 * pi ) / ca0 - cf1
        cpb(1) = cf1
        do k = 2, n0
           cf = ( - z * cf1 + cf0 ) / ( k - 1.0D+00 )
           cpb(k) = cf
           cf0 = cf1
           cf1 = cf
        end do

     else

        if ( a0 <= 3.0D+00 ) then

           call cpdsa ( -n0, z, cfa )
           cpb(n0) = cfa
           n1 = n0 + 1
           call cpdsa ( -n1, z, cfb )
           cpb(n1) = cfb
           nm1 = n0 - 1
           do k = nm1, 0, -1
              cf = z * cfa + ( k + 1.0D+00 ) * cfb
              cpb(k) = cf
              cfb = cfa
              cfa = cf
           end do

        else

           m = 100 + abs ( n )
           cfa = c0
           cfb = cmplx ( 1.0D-30, 0.0D+00, kind = 8 )
           do k = m, 0, -1
              cf = z * cfb + ( k + 1.0D+00 ) * cfa
              if ( k <= n0 ) then
                 cpb(k) = cf
              end if
              cfa = cfb
              cfb = cf
           end do
           cs0 = ca0 / cf
           do k = 0, n0
              cpb(k) = cs0 * cpb(k)
           end do

        end if

     end if

  end if

  cpd(0) = -0.5D+00 * z * cpb(0)

  if ( 0 <= n ) then
     do k = 1, n
        cpd(k) = -0.5D+00 * z * cpb(k) + k * cpb(k-1)
     end do
  else
     do k = 1, n0
        cpd(k) = 0.5D+00 * z * cpb(k) - cpb(k-1)
     end do
  end if

  return
end subroutine cpbdn