lpmns Subroutine

subroutine lpmns ( m, n, x, pm, pd )

  !*****************************************************************************80
  !
  !! LPMNS 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:
  !
  !    18 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 ) M, the order of Pmn(x).
  !
  !    Input, integer ( kind = 4 ) N, the degree of Pmn(x).
  !
  !    Input, real ( kind = 8 ) X, the argument.
  !
  !    Output, real ( kind = 8 ) PM(0:N), PD(0:N), the values and derivatives
  !    of the function from degree 0 to N.
  !
  implicit none

  integer ( kind = 4 ) n

  integer ( kind = 4 ) k
  integer ( kind = 4 ) m
  real ( kind = 8 ) pm(0:n)
  real ( kind = 8 ) pm0
  real ( kind = 8 ) pm1
  real ( kind = 8 ) pm2
  real ( kind = 8 ) pmk
  real ( kind = 8 ) pd(0:n)
  real ( kind = 8 ) x
  real ( kind = 8 ) x0

  do k = 0, n
     pm(k) = 0.0D+00
     pd(k) = 0.0D+00
  end do

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

     do k = 0, n
        if ( m == 0 ) then
           pm(k) = 1.0D+00
           pd(k) = 0.5D+00 * k * ( k + 1.0D+00 )
           if ( x < 0.0D+00 ) then
              pm(k) = ( -1.0D+00 ) ** k * pm(k)
              pd(k) = ( -1.0D+00 ) ** ( k + 1 ) * pd(k)
           end if
        else if ( m == 1 ) then
           pd(k) = 1.0D+300
        else if ( m == 2 ) then
           pd(k) = -0.25D+00 * ( k + 2.0D+00 ) * ( k + 1.0D+00 ) &
                * k * ( k - 1.0D+00 )
           if ( x < 0.0D+00 ) then
              pd(k) = ( -1.0D+00 ) ** ( k + 1 ) * pd(k)
           end if
        end if
     end do
     return
  end if

  x0 = abs ( 1.0D+00 - x * x )
  pm0 = 1.0D+00
  pmk = pm0
  do k = 1, m
     pmk = ( 2.0D+00 * k - 1.0D+00 ) * sqrt ( x0 ) * pm0
     pm0 = pmk
  end do
  pm1 = ( 2.0D+00 * m + 1.0D+00 ) * x * pm0
  pm(m) = pmk
  pm(m+1) = pm1
  do k = m + 2, n
     pm2 = ( ( 2.0D+00 * k - 1.0D+00 ) * x * pm1 &
          - ( k + m - 1.0D+00 ) * pmk ) / ( k - m )
     pm(k) = pm2
     pmk = pm1
     pm1 = pm2
  end do

  pd(0) = ( ( 1.0D+00 - m ) * pm(1) - x * pm(0) ) &
       / ( x * x - 1.0D+00 )  
  do k = 1, n
     pd(k) = ( k * x * pm(k) - ( k + m ) * pm(k-1) ) &
          / ( x * x - 1.0D+00 )
  end do

  return
end subroutine lpmns