lpni Subroutine

subroutine lpni ( n, x, pn, pd, pl )

  !*****************************************************************************80
  !
  !! LPNI computes Legendre polynomials Pn(x), derivatives, and integrals.
  !
  !  Discussion:
  !
  !    This routine computes Legendre polynomials Pn(x), Pn'(x)
  !    and the integral of Pn(t) from 0 to 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:
  !
  !    13 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), PL(0:N), the values, 
  !    derivatives and integrals of the polyomials of degrees 0 to N at X.
  !
  implicit none

  integer ( kind = 4 ) n

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

  pn(0) = 1.0D+00
  pn(1) = x
  pd(0) = 0.0D+00
  pd(1) = 1.0D+00
  pl(0) = x
  pl(1) = 0.5D+00 * x * x
  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

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

     if ( k /= 2 * int ( k / 2 ) ) then

        r = 1.0D+00 / ( k + 1.0D+00 )
        n1 = ( k - 1 ) / 2
        do j = 1, n1
           r = ( 0.5D+00 / j - 1.0D+00 ) * r
        end do
        pl(k) = pl(k) + r

     end if

  end do

  return
end subroutine lpni