      Subroutine vcoul (k,kp,a,rc,r,vl,nls,nif)
c *** see RH Landau, Program lpott, 1981
c     this sub determines partial wave decomp of coulmb pot of spherf er
c     rc and cut off(matching radius) r
c     all units are mev-1
c     a= 3*q1*q2/(2 rc)
      Implicit Real*8 (a-h,k,m,o-z)
      Real*8 xcos(48),wt(48)
      Dimension vl(30), pl(30), ff(4)
c
      data nit/0/
      b = -a/(3*rc*rc)
      If (nit.ne.0) GoTo 10
c     1st call to vcoul. setup ptsw wts
      ncoss = 48
      nit = 1
      Call Gauss (ncoss,011,-1.d0,1.d0,xcos,wt)
 10   rhl1 = kp*kp+k*k
      rhl2 = 2*kp*k
      Do 20 i=1,nls
         vl(i) = 0.
 20   Continue
      Do 70 ncos=1,ncoss
         x = xcos(ncos)
         q = Sqrt(abs(rhl1-rhl2*x))
         If (q.gt.1.) GoTo 30
c     special q=0 form --i think this is ok
         vofq = (0.5/3.141593**2)*(a*(rc**3)/3+b*(rc**5)/5+(rc*a/3)*(r*r
     1   -rc*rc))
         GoTo 50
 30      Continue
         rhl = q*rc
         s = sin(rhl)
         c = cos(rhl)
         q2 = q*q
         If (nif.eq.4) GoTo 40
         q3 = q2*q
         q5 = q2*q3
         vofq = (0.5/3.141593**2)*(a*(s-rhl*c)/q3+b*((3*q2*rc*rc-6)*s-
     1   (q3*rc**3-6*rhl)*c)/q5+2*rc*a*(c-cos(q*r))/(3*q2))
         GoTo 50
c     Realistic form factor
 40      Continue
         nnf = 1
         Call ffact (q2,ff,nnf)
         vofq = (0.5/3.141593**2)*(a*2./3./q/q)*rc*(ff(1)-cos(q*r))
 50      Continue
         pl(1) = wt(ncos)
         pl(2) = x*wt(ncos)
         vl(1) = vl(1)+pl(1)*vofq
         vl(2) = vl(2)+pl(2)*vofq
         ii = nls-1
         If (ii.lt.2) GoTo 70
         Do 60 i=2,ii
            g = x*pl(i)
            pl(i+1) = g-pl(i-1)+g-(g-pl(i-1))/i
            vl(i+1) = vl(i+1)+pl(i+1)*vofq
 60      Continue
 70   Continue
      Do 80 i=1,nls
         vl(i) = vl(i)*(3.141593**2)
 80   Continue
      Return
      End
