clqn Subroutine

subroutine clqn ( n, x, y, cqn, cqd )

  !*****************************************************************************80
  !
  !! CLQN: Legendre function Qn(z) and derivative Wn'(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:
  !
  !    01 August 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 degree of Qn(z).
  !
  !    Input, real ( kind = 8 ) X, Y, the real and imaginary parts of the 
  !    argument Z.
  !
  !    Output, complex ( kind = 8 ) CQN(0:N), CQD(0:N), the values of Qn(z) 
  !    and Qn'(z.
  !
  implicit none

  integer ( kind = 4 ) n

  complex ( kind = 8 ) cq0
  complex ( kind = 8 ) cq1
  complex ( kind = 8 ) cqf0
  complex ( kind = 8 ) cqf1
  complex ( kind = 8 ) cqf2
  complex ( kind = 8 ) cqn(0:n)
  complex ( kind = 8 ) cqd(0:n)
  integer ( kind = 4 ) k
  integer ( kind = 4 ) km
  integer ( kind = 4 ) ls
  real ( kind = 8 ) x
  real ( kind = 8 ) y
  complex ( kind = 8 ) z

  z = cmplx ( x, y, kind = 8 )

  if ( z == 1.0D+00 ) then
     do k = 0, n
        cqn(k) = cmplx ( 1.0D+30, 0.0D+00, kind = 8 )
        cqd(k) = cmplx ( 1.0D+30, 0.0D+00, kind = 8 )
     end do
     return
  end if

  if ( 1.0D+00 < abs ( z ) ) then
     ls = -1
  else
     ls = +1
  end if

  cq0 = 0.5D+00 * log ( ls * ( 1.0D+00 + z ) / ( 1.0D+00 - z ) )
  cq1 = z * cq0 - 1.0D+00
  cqn(0) = cq0
  cqn(1) = cq1

  if ( abs ( z ) < 1.0001D+00 ) then

     cqf0 = cq0
     cqf1 = cq1
     do k = 2, n
        cqf2 = ( ( 2.0D+00 * k - 1.0D+00 ) * z * cqf1 &
             - ( k - 1.0D+00 ) * cqf0 ) / k
        cqn(k) = cqf2
        cqf0 = cqf1
        cqf1 = cqf2
     end do

  else

     if ( 1.1D+00 < abs ( z ) ) then
        km = 40 + n
     else
        km = ( 40 + n ) * int ( - 1.0D+00 &
             - 1.8D+00 * log ( abs ( z - 1.0D+00 ) ) )
     end if

     cqf2 = 0.0D+00
     cqf1 = 1.0D+00
     do k = km, 0, -1
        cqf0 = ( ( 2 * k + 3.0D+00 ) * z * cqf1 &
             - ( k + 2.0D+00 ) * cqf2 ) / ( k + 1.0D+00 )
        if ( k <= n ) then
           cqn(k) = cqf0
        end if
        cqf2 = cqf1
        cqf1 = cqf0
     end do
     do k = 0, n
        cqn(k) = cqn(k) * cq0 / cqf0
     end do
  end if

  cqd(0) = ( cqn(1) - z * cqn(0) ) / ( z * z - 1.0D+00 )
  do k = 1, n
     cqd(k) = ( k * z * cqn(k) - k * cqn(k-1) ) / ( z * z - 1.0D+00 )
  end do

  return
end subroutine clqn