subroutine gaumeh(l,i,rt,eta,nx,index)
c-------------------------------------------------------
c  calculates phase shifts using an nx point gauss mehler method
c  with sqrt(1-x**2) as weight function using form of the wkb phase shift
c  formula given in the notes of rt pack of 5 july 73
c  on first pass with index=0 it calculates the points and weights
c-------------------------------------------------------
      dimension w(24),x(24),sq(24),xs(24)
      common/Blank/ e,b,pk,pl,qsq(19,10),qss(19,10)
      if(index.ne.0) go to 30
      pi = 3.14159265
      p = pi/(2.*float(nx)+1.)
      do 20 j=1,nx
      x(j)=cos(float(j)*p)
      xs(j)=x(j)**2
      d =1.-xs(j)
      sq(j)=sqrt(d)
      w(j)=p*d
      write(6,10) j,x(j),w(j)
   10 format(1h ,'j=',i5,'  x(j)=',f10.7,'  w(j)=',f10.7)
   20 continue
      index=1
   30 eta = 0.
      do 40 j=1,nx
      r=rt/x(j)
      call potent(r,l,i)
      f = (qss(l,i)/sq(j) - pk)/xs(j)
   40 eta = eta + w(j)*f
      eta = rt*eta + (pl+0.5-pk*rt)*0.5*pi
      return
      end