c***********************************************************************
c
c       Oregon State University Nuclear Theory Group
c       Multi-Channel Anti Kaon -- Nucleon Code
c       Guangliang He & Rubin H. Landau
c
c***********************************************************************
c***********************************************************************
c     gauss.F Version 1.2
c     Latest Revision Date 22:38:03 9/12/90
c***********************************************************************
c
      subroutine gauss2(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)
c     x, w     output grid points and weights.
c***********************************************************************
c
      integer*4 npts,job,m,i,j 
      real*8 x(*),w(*),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.0d-14)
      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
         pause 'Wrong value of job'
      endif
c
      return
      end