lpn Subroutine

subroutine lpn ( n, x, pn, pd )

  !*****************************************************************************80
  !
  !! LPN computes Legendre polynomials Pn(x) and derivatives Pn'(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:
  !
  !    07 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 maximum degree.
  !
  !    Input, real ( kind = 8 ) X, the argument.
  !
  !    Output, real ( kind = 8 ) PN(0:N), PD(0:N), the values and derivatives
  !    of the polyomials of degrees 0 to N at X.
  !
  implicit none

  integer ( kind = 4 ) n

  integer ( kind = 4 ) k
  real ( kind = 8 ) p0
  real ( kind = 8 ) p1
  real ( kind = 8 ) pd(0:n)
  real ( kind = 8 ) pf
  real ( kind = 8 ) pn(0:n)
  real ( kind = 8 ) x

  pn(0) = 1.0D+00
  pn(1) = x
  pd(0) = 0.0D+00
  pd(1) = 1.0D+00
  p0 = 1.0D+00
  p1 = x

  do k = 2, n

     pf = ( 2.0D+00 * k - 1.0D+00 ) / k * x * p1 &
          - ( k - 1.0D+00 ) / k * p0
     pn(k) = pf

     if ( abs ( x ) == 1.0D+00 ) then
        pd(k) = 0.5D+00 * x ** ( k + 1 ) * k * ( k + 1.0D+00 )
     else
        pd(k) = k * ( p1 - x * pf ) / ( 1.0D+00 - x * x )
     end if

     p0 = p1
     p1 = pf

  end do

  return
end subroutine lpn