subroutine LegendreQandL(N,x,q,qprime,lN)
implicit none
integer :: N
real(real64) :: x
real(real64) :: lN,q,qprime
! Local
real(real64) :: lNm1,lNm2,dlNm1,dlNm2,dlN,lN1,dlN1
integer :: i
lNm2 = 1.0_real64
lNm1 = x
dlNm2 = 0.0_real64
dlNm1 = 1.0_real64
do i = 2,N
lN = (2.0_real64*i-1.0_real64)/(real(i,real64))*x*lNm1-(real(i,real64)-1.0_real64)/(real(i,real64))*lNm2
dlN = dlNm2+(2.0_real64*real(i,real64)-1.0_real64)*lNm1
lNm2 = lNm1
lNm1 = lN
dlNm2 = dlNm1
dlNm1 = dlN
enddo
i = N+1
lN1 = (2.0_real64*i-1.0_real64)/(real(i,real64))*x*lN-(real(i,real64)-1.0_real64)/(real(i,real64))*lNm2
dlN1 = dlNm2+(2.0_real64*real(i,real64)-1.0_real64)*lNm1
q = lN1-lNm2
qprime = dlN1-dlNm2
endsubroutine LegendreQandL