LegendrePolynomial Subroutine

  subroutine LegendrePolynomial(N,x,lAtX,dLdxAtX)
    implicit none
    integer     :: N
    real(real64)    :: x
    real(real64)    :: lAtX,dLdxAtX
    ! Local
    real(real64) :: lNm1,lNm2,dlNm1,dlNm2
    integer  :: i

    if(N == 0) then

      lAtX = 1.0_real64
      dLdxAtX = 0.0_real64

    elseif(N == 1) then

      lAtX = x
      dLdxAtX = 1.0_real64

    else

      lnM2 = 1.0_real64
      lnM1 = x
      dlnM2 = 0.0_real64
      dlnM1 = 1.0_real64

      do i = 2,N

        lAtX = ((2.0_real64*real(i,real64)-1.0_real64)*x*lnM1- &
                (real(i,real64)-1.0_real64)*lnM2)/(real(i,real64))

        dldxAtX = dlnM2+(2.0_real64*real(i,real64)-1.0_real64)*lnM1
        lnM2 = lnM1
        lnM1 = lAtX
        dlnM2 = dlnM1
        dlnM1 = dldxAtX

      enddo

    endif

  endsubroutine LegendrePolynomial