c     QGAUSL
c
c     Source: Press, et. al.
c
      subroutine qgausl(n,x1,x2,x,w)
      implicit double precision (a-h,o-z)
      double precision x(n),w(n),x1,x2
      eps=1.0e-8
      m=(n+1)/2
      xm=0.5*(x2+x1)
      xl=0.5*(x2-x1)
      do 12 i=1,m
        z = cos(3.141592654*(i-.25)/(n+.5))
1       continue
          p1 = 1.0
          p2 = 0.0
          do 11 j=1,n
            p3 = p2
            p2 = p1
            p1 = ((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j
 11       continue
          pp = n*(z*p1-p2)/(z*z-1.0)
          z1 = z
          z = z1-p1/pp
          if (dabs(z-z1).gt.eps) go to 1
          x(i) = xm-xl*z
          x(n+1-i) = xm+xl*z
          w(i) = 2.0*xl/((1.0-z*z)*pp*pp)
          w(n+1-i) = w(i)
 12   continue
      return
      end