subroutine Gauss(npts,job,a,b,x,w)
c
c***********************************************************************
c     rescale rescals the Gauss-Legendre grid points and weights
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), with
c              rational scaling
c     job = 3  for integral (a, inf) with 50% inside (a, b+2a), with
c              atan() scaling. (a+b) has to great than zero.
c
c     x, w     output grid points and weights.
c***********************************************************************
c
      integer npts,job,m,i,j 
      double precision x(npts),w(npts),a,b,xi
      double precision t,t1,pp,p1,p2,p3,aj
      double precision eps,pi,zero,two,one,half,quarter,pio4
      parameter (pi = 3.14159265358979323846264338328d0, eps = 3.0d-15)
      parameter (zero=0.0d0,one=1.0d0,two=2.0d0)
      parameter (half=0.5d0,quarter=0.25d0)
c
c *** FIRST EXECTUABLE *************************************************
c
      m=(npts+1)/2
      do 10020 i=1,m
         t=cos(pi*(i-quarter)/(npts+half))
10000 continue
         p1=one
         p2=zero
         aj=zero
         do 10010 j=1,npts
            p3=p2
            p2=p1
            aj=aj+one
            p1=((two*aj-one)*t*p2-(aj-one)*p3)/aj
10010    continue
         pp=npts*(t*p1-p2)/(t*t-one)
         t1=t
         t=t1-p1/pp
c
         if(abs(t-t1).gt.eps) goto 10000
c
         x(i)=-t
         x(npts+1-i)=t
         w(i)=two/((one-t*t)*pp*pp)
         w(npts+1-i)=w(i)
10020 continue
c
c***********************************************************************
c     rescale the grid points 
c***********************************************************************
c
      if (job.eq.0) then
c     scale to (a,b) uniformly
         do 10030 i=1,npts
            x(i)=x(i)*(b-a)/two+(b+a)/two
            w(i)=w(i)*(b-a)/two
10030    continue
c
      elseif (job.eq.1) then
c     scale to (0,b) with 50% points inside (0,ab/(a+b))
         do 10040 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))
10040    continue
c
      elseif (job.eq.2) then
c     scale to (a,inf) with 50% points inside (a,b+2a)
         do 10050 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))
10050    continue
      else if (job .eq. 3) then
c     scale to (a, inf) with 50% points inside (a, b+2a)
         pio4 = atan(one)
         if (a+b .le. zero) stop 'Gauss: a+b less than zero'
         do 10060 i = 1, npts
            xi = x(i)
            x(i) = a+(a+b)*tan(pio4*(one+xi))
            w(i) = w(i)*(a+b)*pio4
     &           /(cos(pio4*(one+xi))*cos(pio4*(one+xi)))
10060    continue
      else
         stop 'Gauss: Wrong value of job'
      endif
c
      return
      end