      Subroutine intgl (nch)
c *** see RH Landau, Program lpott, 1981
      Implicit  Real*8   (a-h,o-z)
      Common /enpaul/ pf,ak,akd,cskkd,xi,e,xon
      Common /enint/ a,b,lPauli
      Common /deint/ ans
      Common /enqav/ xip
      Real*8 x(32),w(32),sum,dum,ampi,ampisq,amn,ek,k,ksq
      Real*8 amu,const,won,xx(16),ww(16)
c     Real dlog
c     Real Sqrt
      Complex*16 ans,ci
c     Complex cmplx
      logical lin,lPauli,l1,l2
      data lin,ci/.false.,(0.0,1.0)/
      If (lin) GoTo 10
      lin = .true.
      nhf = 16
c     rhl change of points so can run cex at osu
      nhf = 4
      mesh = nhf
      Call Gauss (mesh,012,-1.d0,1.d0,xx,ww)
c     here we have set up mesh on (-1,+1)
      ampi = 137.3/197.33e0
      ampisq = ampi**2
      amn = 938.91/197.33e0
      pi = 3.14159265
 10   Continue
c
c     scale quad pts..
      along = (b-a)/2.0
      amid = (a+b)/2.0
      Do 20 i=1,mesh
         x(i) = xx(i)*along+amid
         w(i) = ww(i)*along
 20   Continue
c
c     now test If we have a singularity..
      If (e.lt.0.0) GoTo 110
      l1 = xon.gt.a
      l2 = xon.lt.b
      If (l1.and.l2) GoTo 30
      GoTo 110
c
 30   Continue
c     here we have a pole in the integral //////  1st Do principal value
      won = xon**2
      won = Sqrt(ampisq+won)
      amu = amn*won/(amn+won)
      const = xon*amu*g(xon,nch)**2
      If (lPauli) const = const*qav(xon)
c
c     changed version with redistributed quadrature points for this cas
      dxon = xon
      If (xon.lt.amid) GoTo 60
      bpr = b
      apr = 2.*dxon-bpr
      alpr = (bpr-apr)/2.
      arst = a
      brst = apr
      alrst = (brst-arst)/2.
      amrst = (arst+brst)/2.
      Do 40 i=1,mesh
         x(i) = xx(i)*alrst+amrst
         w(i) = ww(i)*alrst
 40   Continue
      Do 50 i=1,mesh
         ii = i+mesh
         x(ii) = xx(i)*alpr+dxon
         w(ii) = ww(i)*alpr
 50   Continue
      GoTo 90
 60   Continue
      apr = a
      bpr = 2.*dxon-apr
      alpr = (bpr-apr)/2.
      arst = bpr
      brst = b
      alrst = (brst-arst)/2.
      amrst = (arst+brst)/2.
      Do 70 i=1,mesh
         x(i) = xx(i)*alpr+dxon
         w(i) = ww(i)*alpr
 70   Continue
      Do 80 i=1,mesh
         ii = i+mesh
         x(ii) = xx(i)*alrst+amrst
         w(ii) = ww(i)*alrst
 80   Continue
 90   nmesh = 2*mesh
      sum = 0.0e0
      nmesh = 2*mesh
      Do 100 i=1,nmesh
         k = x(i)
         ksq = k*k
         y = k
         ek = Sqrt(ampisq+ksq)-ampi+ksq/(2.0*amn)
         dum = ksq*g(y,nch)**2
         If (lPauli) dum = dum*qav(y)
         dum = dum/(e-ek)
         dum = dum-const/(xon-k)
         sum = sum+w(i)*dum
 100  Continue
c
c     add the correction integral..
      dum = xon-a
      dum = dum/(b-xon)
      y = sum+log(dum)*const
c
c     and finally the pole term..
      y1 = -pi*const
      ans = cmplx(y,y1)
      Return
c
c     here there is no singularity..
 110  Continue
      sum = 0.0
      Do 120 i=1,mesh
         k = x(i)
         ksq = k**2
         y = k
         ek = Sqrt(ampisq+ksq)-ampi+ksq/(2.0*amn)
         dum = ksq*g(y,nch)**2
         If (lPauli) dum = dum*qav(y)
         sum = sum+w(i)*dum/(e-ek)
 120  Continue
      y = sum
      xzero=0.d0
      ans = cmplx(y,xzero)
      Return
      End
