clpn Subroutine

subroutine clpn ( n, x, y, cpn, cpd )

  !*****************************************************************************80
  !
  !! CLPN computes Legendre functions and derivatives for complex argument.
  !
  !  Discussion:
  !
  !    Compute Legendre polynomials Pn(z) and their derivatives Pn'(z) for 
  !    a 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:
  !
  !    15 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 degree.
  !
  !    Input, real ( kind = 8 ) X, Y, the real and imaginary parts 
  !    of the argument.
  !
  !    Output, complex ( kind = 8 ) CPN(0:N), CPD(0:N), the values of Pn(z)
  !    and Pn'(z).
  !
  implicit none

  integer ( kind = 4 ) n

  complex ( kind = 8 ) cp0
  complex ( kind = 8 ) cp1
  complex ( kind = 8 ) cpd(0:n)
  complex ( kind = 8 ) cpf
  complex ( kind = 8 ) cpn(0:n)
  integer ( kind = 4 ) k
  real ( kind = 8 ) x
  real ( kind = 8 ) y
  complex ( kind = 8 ) z

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

  cpn(0) = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
  cpn(1) = z
  cpd(0) = cmplx ( 0.0D+00, 0.0D+00, kind = 8 )
  cpd(1) = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )

  cp0 = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
  cp1 = z
  do k = 2, n
     cpf = ( 2.0D+00 * k - 1.0D+00 ) / k * z * cp1 - ( k - 1.0D+00 ) / k * cp0
     cpn(k) = cpf
     if ( abs ( x ) == 1.0D+00 .and. y == 0.0D+00 ) then
        cpd(k) = 0.5D+00 * x ** ( k + 1 ) * k * ( k + 1.0D+00 )
     else
        cpd(k) = k * ( cp1 - z * cpf ) / ( 1.0D+00 - z * z )
     end if
     cp0 = cp1
     cp1 = cpf
  end do

  return
end subroutine clpn