lpmn Subroutine

subroutine lpmn ( mm, m, n, x, pm, pd )

  !*****************************************************************************80
  !
  !! LPMN computes associated Legendre functions Pmn(X) and derivatives P'mn(x).
  !
  !  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:
  !
  !    19 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 ) MM, the leading dimension of PM and PD.
  !
  !    Input, integer ( kind = 4 ) M, the order of Pmn(x).
  !
  !    Input, integer ( kind = 4 ) N, the degree of Pmn(x).
  !
  !    Input, real ( kind = 8 ) X, the argument of Pmn(x).
  !
  !    Output, real ( kind = 8 ) PM(0:MM,0:N), PD(0:MM,0:N), the
  !    values of Pmn(x) and Pmn'(x).
  !
  implicit none

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

  integer ( kind = 4 ) i
  integer ( kind = 4 ) j
  integer ( kind = 4 ) ls
  integer ( kind = 4 ) m
  real ( kind = 8 ) pd(0:mm,0:n)
  real ( kind = 8 ) pm(0:mm,0:n)
  real ( kind = 8 ) x
  real ( kind = 8 ) xq
  real ( kind = 8 ) xs

  do i = 0, n
     do j = 0, m
        pm(j,i) = 0.0D+00
        pd(j,i) = 0.0D+00
     end do
  end do

  pm(0,0) = 1.0D+00

  if ( abs ( x ) == 1.0D+00 ) then

     do i = 1, n
        pm(0,i) = x ** i
        pd(0,i) = 0.5D+00 * i * ( i + 1.0D+00 ) * x ** ( i + 1 )
     end do

     do j = 1, n
        do i = 1, m
           if ( i == 1 ) then
              pd(i,j) = 1.0D+300
           else if ( i == 2 ) then
              pd(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 ( x ) ) then
     ls = -1
  else
     ls = +1
  end if

  xq = sqrt ( ls * ( 1.0D+00 - x * x ) )
  xs = ls * ( 1.0D+00 - x * x )
  do i = 1, m
     pm(i,i) = - ls * ( 2.0D+00 * i - 1.0D+00 ) * xq * pm(i-1,i-1)
  end do

  do i = 0, m
     pm(i,i+1) = ( 2.0D+00 * i + 1.0D+00 ) * x * pm(i,i)
  end do

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

  pd(0,0) = 0.0D+00
  do j = 1, n
     pd(0,j) = ls * j * ( pm(0,j-1) - x * pm(0,j) ) / xs
  end do

  do i = 1, m
     do j = i, n
        pd(i,j) = ls * i * x * pm(i,j) / xs + ( j + i ) &
             * ( j - i + 1.0D+00 ) / xq * pm(i-1,j)
     end do
  end do

  return
end subroutine lpmn