subroutine LegendreGauss(N,nodes,weights)
implicit none
integer :: N
real(real64) :: nodes(0:N)
real(real64) :: weights(0:N)
! Local
real(real64) :: nodes_local(0:N)
real(real64) :: weights_local(0:N)
real(real64) :: lN1,dlN1
real(real64) :: delta
integer :: j,kIt
if(N == 0) then
nodes_local(0) = 0.0_real64
weights_local(0) = 2.0_real64
elseif(N == 1) then
nodes_local(0) = -sqrt(1.0_real64/3.0_real64)
weights_local(0) = 1.0_real64
nodes_local(1) = -nodes(0)
weights_local(1) = weights(0)
else
do j = 0,((N+1)/2)
nodes_local(j) = -cos((2.0_real64*real(j,real64)+1.0_real64)*pi/(2.0_real64*real(N,real64)+1.0_real64))
do kIt = 1,newtonMax
call LegendrePolynomial(N+1,nodes_local(j),lN1,dlN1)
delta = -lN1/dlN1
nodes_local(j) = nodes_local(j)+delta
if(abs(delta) <= TOL*nodes_local(j)) exit
enddo
call LegendrePolynomial(N+1,nodes_local(j),lN1,dlN1)
weights_local(j) = 2.0_real64/((1.0_real64-nodes_local(j)*nodes_local(j))*dlN1*dlN1)
weights_local(N-j) = weights_local(j)
nodes_local(N-j) = -nodes_local(j)
enddo
endif
if(mod(real(N,real64),2.0_real64) == 0.0_real64) then
call LegendrePolynomial(N+1,0.0_real64,lN1,dlN1)
nodes_local(N/2) = 0.0_real64
weights_local(N/2) = 2.0/(dlN1*dlN1)
endif
do j = 0,N
nodes(j) = real(nodes_local(j),real64)
weights(j) = real(weights_local(j),real64)
enddo
endsubroutine LegendreGauss