subroutine qgausl (n,x1,x2, x,w)
!
!  Qgausl computes abscissas and weights for Legendre Gaussian quadrature.
!  Inputs are the number of quadrature points n and the interval of
!  integration: [x1,x2].
!
   integer, intent(in)                  :: n
   real, intent(in)                     :: x1,x2
   real, dimension(n), intent(out)      :: x,w
   real                                 :: eps= 1.0e-8, pi=3.141592654, &
                                           xm,xl,p1,p2,p3,pp,z,z1
   integer                              :: i,j,m

   m  = (n+1)/2
   xm = 0.5*(x2+x1)
   xl = 0.5*(x2-x1)
   do i=1,m
     z = cos(pi*(i-.25)/(n+.5))
     newton: do
       p1  = 1.0
       p2  = 0.0
       do j=1,n
         p3 = p2
         p2 = p1
         p1 = ((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j
       end do
       pp  = n*(z*p1-p2)/(z*z-1.0)
       z1  = z
       z   = z1-p1/pp
       if (abs(z-z1) <= eps) exit newton
     end do newton
     x(i)     = xm - xl * z
     x(n-i+1) = xm + xl * z
     w(i)     = 2.0*xl/((1.0-z*z)*pp*pp)
     w(n-i+1) = w(i)

   end do

end subroutine qgausl