c  From: "A SURVEY OF COMPUTATIONAL PHYSICS" 
c  by RH Landau, MJ Paez, and CC BORDEIANU 
c  Copyright Princeton University Press, Princeton, 2008.
c  Electronic Materials copyright: R Landau, Oregon State Univ, 2008;
c  MJ Paez, Univ Antioquia, 2008; and CC BORDEIANU, Univ Bucharest, 2008.
c  Support by National Science Foundation   
c
c  gauss.f: Points and weights for Gaussian quadrature  
c           error message if subroutine called without a main
c           this file must reside in same directory as integ.c 
c     npts = number of points
c     job  = 0  rescalling uniformly between (a, b)
c      1  for integral (0, b) with 50% points inside (0, ab/(a+b))
c      2  for integral (a, inf) with 50% inside (a, b+2a)
c     x, w     output integration points and weights.
c
      subroutine gauss(npts, job, a, b, x, w) 
			
        integer npts, job, m, i, j 
        real*8 x(npts), w(npts), a, b, xi
        real*8 t, t1, pp, p1, p2, p3, aj
        real*8 eps, pi, zero, two, one, half, quarter
        parameter (pi =3.14159265358979323846264338328, eps = 3.E-14)
        parameter (zero = 0.d0, one = 1.d0, two = 2.d0)
        parameter (half = 0.5d0, quarter = 0.25d0)
				
        m = (npts+1)/2
        do 1020 i = 1, m
          t = cos(pi*(i-quarter)/(npts+half))
 1000   continue
        p1 = one
        p2 = zero
        aj = zero
        do 1010 j = 1, npts
          p3 = p2
          p2 = p1
          aj = aj+one
          p1 = ((two*aj-one)*t*p2-(aj-one)*p3)/aj
 1010   continue
        pp = npts*(t*p1-p2)/(t*t-one)
        t1 = t
        t = t1-p1/pp
        if ( abs(t-t1) .gt. eps ) goto 1000
        x(i) = -t
        x(npts+1-i) = t
        w(i) = two/((one-t*t)*pp*pp)
        w(npts+1-i) = w(i)
 1020   continue
c                 rescales the gauss-legendre grid points and weights
        if (job .eq. 0) then
c                                           scale to (a, b) uniformly
          do 1030 i = 1, npts
            x(i) = x(i)*(b-a)/two+(b+a)/two
            w(i) = w(i)*(b-a)/two
 1030     continue
          elseif ( job .eq. 1 ) then
c                scale to (0, b) with 50% points inside (0, ab/(a+b))
            do 1040 i = 1, npts
             xi = x(i)
              x(i) = a*b*(one+xi)/(b+a-(b-a)*xi)
              w(i) = w(i)*two*a*b*b/((b+a-(b-a)*xi)*(b+a-(b-a)*xi))
 1040       continue
          elseif ( job .eq. 2 ) then
c                  scale to (a, inf) with 50% points inside (a, b+2a)
          do 1050 i = 1, npts
            xi = x(i)
            x(i) = (b*xi+b+a+a)/(one-xi)
            w(i) = w(i)*two*(a+b)/((one-xi)*(one-xi))
 1050     continue
          else
            pause 'Wrong value of job'
        endif
        Return
        end