clqn Subroutine

   subroutine clqn(n,x,y,Cqn,Cqd)

      integer,intent(in) :: n !! Degree of Qn(z), n = 0,1,2,...
      real(wp),intent(in) :: x !! Real part of z
      real(wp),intent(in) :: y !! Imaginary part of z
      complex(wp),intent(out) :: Cqn(0:n) !! `Cqn`
      complex(wp),intent(out) :: Cqd(0:n) !! `Cqd`

      complex(wp) :: cq0 , cq1 , cqf0 , cqf1 , cqf2 , z
      integer :: k , km , ls

      z = cmplx(x,y,kind=wp)
      if ( z==1.0_wp ) then
         do k = 0 , n
            Cqn(k) = (1.0e+300_wp,0.0_wp)
            Cqd(k) = (1.0e+300_wp,0.0_wp)
         enddo
         return
      endif
      ls = 1
      if ( abs(z)>1.0_wp ) ls = -1
      cq0 = 0.5_wp*log(ls*(1.0_wp+z)/(1.0_wp-z))
      cq1 = z*cq0 - 1.0_wp
      Cqn(0) = cq0
      Cqn(1) = cq1
      if ( abs(z)<1.0001_wp ) then
         cqf0 = cq0
         cqf1 = cq1
         do k = 2 , n
            cqf2 = ((2.0_wp*k-1.0_wp)*z*cqf1-(k-1.0_wp)*cqf0)/k
            Cqn(k) = cqf2
            cqf0 = cqf1
            cqf1 = cqf2
         enddo
      else
         if ( abs(z)>1.1_wp ) then
            km = 40 + n
         else
            km = (40+n)*int(-1.0_wp-1.8_wp*log(abs(z-1.0_wp)))
         endif
         cqf2 = 0.0_wp
         cqf1 = 1.0_wp
         do k = km , 0 , -1
            cqf0 = ((2*k+3.0_wp)*z*cqf1-(k+2.0_wp)*cqf2)/(k+1.0_wp)
            if ( k<=n ) Cqn(k) = cqf0
            cqf2 = cqf1
            cqf1 = cqf0
         enddo
         do k = 0 , n
            Cqn(k) = Cqn(k)*cq0/cqf0
         enddo
      endif
      Cqd(0) = (Cqn(1)-z*Cqn(0))/(z*z-1.0_wp)
      do k = 1 , n
         Cqd(k) = (k*z*Cqn(k)-k*Cqn(k-1))/(z*z-1.0_wp)
      enddo

   end subroutine clqn