subroutine VCoul(z1, z2, l, npts, p, w, ldv, vc)
c
c***********************************************************************
c     VCoul calculates coulomb potential matrix from l'th partial wave
c     between two particles z1, z2, in a way suitable for numerical 
c     solution of Lippman-Schwinger Equations.
c
c     On input:
c     z1, z2: The charge number of two interaction particles
c     l:      The partial wave number
c     npts:   The number of grid points for Lippman-Schwinger Eq.
c     p, w:   The grid points and weight for LS Eq.
c     ldv:    The leading dimensions of vc
c     On output:
c     vc(ldv, 1dv) The Coulomb potential.
c***********************************************************************
c
      integer z1, z2, l, npts, ldv
      double precision p(npts), w(npts), vc(ldv, 1)
      double precision alpha, pi, two, zero, charge2, temp
      integer ipt, jpt
c
      parameter (alpha = 7.29735307639485d-3, pi = 3.141592653589793d0)
      parameter (two = 2.0d0, zero = 0.0d0)
c
c***********************************************************************
c
      charge2 = z1*z2*alpha
c
      if (l .eq. 0) then
c     construct upper half of vc
         do jpt = 1, npts
            do ipt = 1, jpt-1
               vc(ipt, jpt) = charge2/(two*p(ipt)*p(jpt))
     &              *log((p(jpt)+p(ipt))/(p(jpt)-p(ipt)))
            enddo
         enddo
c     lower half
         do jpt = 1, npts
            do ipt = jpt+1, npts
               vc(ipt, jpt) = vc(jpt, ipt)
            enddo
         enddo
c     diagnal elements
         do jpt = 1, npts
            vc(jpt, jpt) = zero
            temp = zero
            do ipt = 1, npts
               temp = temp + w(ipt)*vc(jpt, ipt)
            enddo
            vc(jpt, jpt) = (charge2*pi*pi/(two*two*p(jpt))-temp)/w(jpt)
         enddo
      else
         write(8, *) 'only l = 0 has been implimented'
         stop
      endif
c
      return
      end