cpdsa Subroutine

      subroutine cpdsa(n,z,Cdn)

      integer,intent(in) :: n !! Order of D(z) (n = 0,-1,-2,...)
      complex(wp),intent(in) :: z !! complex argument of D(z)
      complex(wp),intent(out) :: Cdn !! Dn(z)

      complex(wp) :: ca0 , cb0 , cdw , cr
      real(wp) :: g0 , g1 , ga0 , gm , pd , va0 , vm , vt , xn
      integer :: m

      real(wp),parameter :: eps = 1.0e-15_wp

      ca0 = exp(-0.25_wp*z*z)
      va0 = 0.5_wp*(1.0_wp-n)
      if ( n==0 ) then
         Cdn = ca0
      elseif ( abs(z)==0.0_wp ) then
         if ( va0<=0.0_wp .and. va0==int(va0) ) then
            Cdn = 0.0_wp
         else
            call gaih(va0,ga0)
            pd = sqrtpi/(2.0_wp**(-0.5_wp*n)*ga0)
            Cdn = cmplx(pd,0.0_wp,wp)
         endif
      else
         xn = -n
         call gaih(xn,g1)
         cb0 = 2.0_wp**(-0.5_wp*n-1.0_wp)*ca0/g1
         vt = -0.5_wp*n
         call gaih(vt,g0)
         Cdn = cmplx(g0,0.0_wp,wp)
         cr = (1.0_wp,0.0_wp)
         do m = 1 , 250
            vm = 0.5_wp*(m-n)
            call gaih(vm,gm)
            cr = -cr*sq2*z/m
            cdw = gm*cr
            Cdn = Cdn + cdw
            if ( abs(cdw)<abs(Cdn)*eps ) exit
         enddo
         Cdn = cb0*Cdn
      endif
   end subroutine cpdsa