cpdsa Subroutine

subroutine cpdsa ( n, z, cdn )

  !*****************************************************************************80
  !
  !! CPDSA computes complex parabolic cylinder function Dn(z) for small 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 ) CDN, the value of DN(z).
  !
  implicit none

  complex ( kind = 8 ) ca0
  complex ( kind = 8 ) cb0
  complex ( kind = 8 ) cdn
  complex ( kind = 8 ) cdw
  complex ( kind = 8 ) cr
  real ( kind = 8 ) eps
  real ( kind = 8 ) g0
  real ( kind = 8 ) g1
  real ( kind = 8 ) ga0
  real ( kind = 8 ) gm
  integer ( kind = 4 ) m
  integer ( kind = 4 ) n
  real ( kind = 8 ) pd
  real ( kind = 8 ) pi
  real ( kind = 8 ) sq2
  real ( kind = 8 ) va0
  real ( kind = 8 ) vm
  real ( kind = 8 ) vt
  real ( kind = 8 ) xn
  complex ( kind = 8 ) z

  eps = 1.0D-15
  pi = 3.141592653589793D+00
  sq2 = sqrt ( 2.0D+00 )
  ca0 = exp ( - 0.25D+00 * z * z )
  va0 = 0.5D+00 * ( 1.0D+00 - n )

  if ( n == 0 ) then

     cdn = ca0

  else

     if ( abs ( z ) == 0.0D+00 ) then

        if ( va0 <= 0.0D+00 .and. va0 == int ( va0 ) ) then
           cdn = 0.0D+00
        else
           call gaih ( va0, ga0 )
           pd = sqrt ( pi ) / ( 2.0D+00 ** ( -0.5D+00 * n ) * ga0 )
           cdn = cmplx ( pd, 0.0D+00, kind = 8 )
        end if

     else

        xn = - n
        call gaih ( xn, g1 )
        cb0 = 2.0D+00 ** ( -0.5D+00 * n - 1.0D+00 ) * ca0 / g1
        vt = -0.5D+00 * n
        call gaih ( vt, g0 )
        cdn = cmplx ( g0, 0.0D+00, kind = 8 )
        cr = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )

        do m = 1, 250
           vm = 0.5D+00 * ( m - n )
           call gaih ( vm, gm )
           cr = - cr * sq2 * z / m
           cdw = gm * cr
           cdn = cdn + cdw
           if ( abs ( cdw ) < abs ( cdn ) * eps ) then
              exit
           end if
        end do

        cdn = cb0 * cdn

     end if

  end if

  return
end subroutine cpdsa