clpmn Subroutine

subroutine clpmn ( mm, m, n, x, y, cpm, cpd )

  !*****************************************************************************80
  !
  !! CLPMN: associated Legendre functions and derivatives for complex argument.
  !
  !  Discussion:
  !
  !    Compute the associated Legendre functions Pmn(z)   
  !    and their derivatives Pmn'(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:
  !
  !    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 ) MM, the physical dimension of CPM and CPD.
  !
  !    Input, integer ( kind = 4 ) M, N, the order and degree of Pmn(z).
  !
  !    Input, real ( kind = 8 ) X, Y, the real and imaginary parts of 
  !    the argument Z.
  !
  !    Output, complex ( kind = 8 ) CPM(0:MM,0:N), CPD(0:MM,0:N), the values of
  !    Pmn(z) and Pmn'(z).
  !
  implicit none

  integer ( kind = 4 ) mm

  complex ( kind = 8 ) cpd(0:mm,0:n)
  complex ( kind = 8 ) cpm(0:mm,0:n)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) j
  integer ( kind = 4 ) ls
  integer ( kind = 4 ) m
  integer ( kind = 4 ) n
  real ( kind = 8 ) x
  real ( kind = 8 ) y
  complex ( kind = 8 ) z
  complex ( kind = 8 ) zq
  complex ( kind = 8 ) zs

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

  do i = 0, n
     do j = 0, m
        cpm(j,i) = cmplx ( 0.0D+00, 0.0D+00, kind = 8 )
        cpd(j,i) = cmplx ( 0.0D+00, 0.0D+00, kind = 8 )
     end do
  end do

  cpm(0,0) = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )

  if ( abs ( x ) == 1.0D+00 .and. y == 0.0D+00 ) then

     do i = 1, n
        cpm(0,i) = x ** i
        cpd(0,i) = 0.5D+00 * i * ( i + 1 ) * x ** ( i + 1 )
     end do

     do j = 1, n
        do i = 1, m
           if ( i == 1 ) then
              cpd(i,j) = cmplx ( 1.0D+30, 0.0D+00, kind = 8 )
           else if ( i == 2 ) then
              cpd(i,j) = -0.25D+00 &
                   * ( j + 2 ) * ( j + 1 ) * j * ( j - 1 ) * x ** ( j + 1 )
           end if
        end do
     end do

     return

  end if

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

  zq = sqrt ( ls * ( 1.0D+00 - z * z ) )
  zs = ls * ( 1.0D+00 - z * z )
  do i = 1, m
     cpm(i,i) = -ls * ( 2.0D+00 * i - 1.0D+00 ) * zq * cpm(i-1,i-1)
  end do
  do i = 0, m
     cpm(i,i+1) = ( 2.0D+00 * i + 1.0D+00 ) * z * cpm(i,i)
  end do

  do i = 0, m
     do j = i + 2, n
        cpm(i,j) = ( ( 2.0D+00 * j - 1.0D+00 ) * z * cpm(i,j-1) &
             - ( i + j - 1.0D+00 ) * cpm(i,j-2) ) / ( j - i )
     end do
  end do

  cpd(0,0) = cmplx ( 0.0D+00, 0.0D+00, kind = 8 )
  do j = 1, n
     cpd(0,j) = ls * j * ( cpm(0,j-1) - z * cpm(0,j) ) / zs
  end do

  do i = 1, m
     do j = i, n
        cpd(i,j) = ls * i * z * cpm(i,j) / zs &
             + ( j + i ) * ( j - i + 1.0D+00 ) / zq * cpm(i-1,j)
     end do
  end do

  return
end subroutine clpmn