clqmn Subroutine

subroutine clqmn ( mm, m, n, x, y, cqm, cqd )

  !*****************************************************************************80
  !
  !! CLQMN: associated Legendre functions and derivatives for complex argument.
  !
  !  Discussion:
  !
  !    This procedure computes the associated Legendre functions of the second 
  !    kind, Qmn(z) and Qmn'(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:
  !
  !    02 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 CQM and CQD.
  !
  !    Input, integer ( kind = 4 ) M, N, the order and degree of Qmn(z).
  !
  !    Input, real ( kind = 8 ) X, Y, the real and imaginary parts of the 
  !    argument Z.
  !
  !    Output, complex ( kind = 8 ) CQM(0:MM,0:N), CQD(0:MM,0:N), the values of
  !    Qmn(z) and Qmn'(z).
  !
  implicit none

  integer ( kind = 4 ) mm
  integer ( kind = 4 ) n 

  complex ( kind = 8 ) cq0
  complex ( kind = 8 ) cq1
  complex ( kind = 8 ) cq10
  complex ( kind = 8 ) cqf
  complex ( kind = 8 ) cqf0
  complex ( kind = 8 ) cqf1
  complex ( kind = 8 ) cqf2
  complex ( kind = 8 ) cqm(0:mm,0:n)
  complex ( kind = 8 ) cqd(0:mm,0:n)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) j
  integer ( kind = 4 ) k
  integer ( kind = 4 ) km
  integer ( kind = 4 ) ls
  integer ( kind = 4 ) m
  real ( kind = 8 ) x
  real ( kind = 8 ) xc
  real ( kind = 8 ) y
  complex ( kind = 8 ) z
  complex ( kind = 8 ) zq
  complex ( kind = 8 ) zs

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

  if ( abs ( x ) == 1.0D+00 .and. y == 0.0D+00 ) then
     do i = 0, m
        do j = 0, n
           cqm(i,j) = cmplx ( 1.0D+30, 0.0D+00, kind = 8 )
           cqd(i,j) = cmplx ( 1.0D+30, 0.0D+00, kind = 8 )
        end do
     end do
     return
  end if

  xc = abs ( z )

  if ( imag ( z ) == 0.0D+00 .or. xc < 1.0D+00 ) then
     ls = 1
  end if

  if ( 1.0D+00 < xc ) then
     ls = -1
  end if

  zq = sqrt ( ls * ( 1.0D+00 - z * z ) )
  zs = ls * ( 1.0D+00 - z * z )
  cq0 = 0.5D+00 * log ( ls * ( 1.0D+00 + z ) / ( 1.0D+00 - z ) )

  if ( xc < 1.0001D+00 ) then

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

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

  else

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

     cqf2 = cmplx ( 0.0D+00, 0.0D+00, kind = 8 )
     cqf1 = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
     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
           cqm(0,k) = cqf0
        end if
        cqf2 = cqf1
        cqf1 = cqf0
     end do

     do k = 0, n
        cqm(0,k) = cq0 * cqm(0,k) / cqf0
     end do

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

     cq10 = -1.0D+00 / zq
     do k = 0, n 
        cqm(1,k) = cq10 * cqm(1,k) / cqf0
     end do

     do j = 0, n
        cq0 = cqm(0,j)
        cq1 = cqm(1,j)
        do i = 0, m - 2
           cqf = -2.0D+00 * ( i + 1 ) * z / zq * cq1 &
                + ( j - i ) * ( j + i + 1.0D+00 ) * cq0
           cqm(i+2,j) = cqf
           cq0 = cq1
           cq1 = cqf
        end do
     end do

  end if

  cqd(0,0) = ls / zs
  do j = 1, n
     cqd(0,j) = ls * j * ( cqm(0,j-1) - z * cqm(0,j) ) / zs
  end do

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

  return
end subroutine clqmn