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