function vcoul (q)
c *** see r h landau, program lpott, 1981
c     this sub determines Coulomb potential of sphere or realistic rho(r)
c     rcut is the cut off (matching radius) r
c     all units are mev-1
      implicit real*8 (a-h,k,m,o-z)
      real*8 Rcoul, Rcut
      common /params/   hbarc, pi, Mp, MN, nz, na, nes, nwaves
      common /switch/   nifty
      common /ranges/ Rcoul, Rcut
      dimension ff(4), nifty(20)
c
         vcoul=0.0
c    ********* Below pt coulomb **********************
          eps = 1.d-4
          rrcut = Rcut/hbarc
          qr = q * rrcut
          q2 = q * q
          if (qr .lt. eps) then
            vcoul = 0.5d0 * rrcut**2
          else
            vcoul = (1.d0 - DCOS(qr))/q
          endif
         vcoul =  0.5d0*nz/137.036/pi/pi*vcoul
c        return
c         Above pt coulomb switch **********************
         if(nifty(1).eq.9)return
         zalpha=nz/137.036
         if (q.gt.1.) go to 30
c     special q=0 form --i think this is ok
c        vcoul = rcut*rcut*zalpha/4./pi/pi/hbarc/hbarc
************* Mr. Lu form of above *******************
         rcut2 = Rcut*Rcut/hbarc/hbarc
         rcoul2 = Rcoul*Rcoul/hbarc/hbarc
         vcoul = zalpha/4./pi/pi*(rcut2 - rcoul2 * .2)
******************************************************
         go to 50
 30      continue
c
c must be define rcoul to use uniform sphere
         rhl = q*Rcoul/hbarc
         s = sin(rhl)
         c = cos(rhl)
         if (nifty(10).eq.4) go to 40
         ff(1)=(3./rhl/rhl)*(s/rhl-c)
         go to 45
c     realistic form factor
 40      continue
         nnf = 1
         q2 = q*q
         call ffact (q2,ff,nnf)
***** tvm 10/27/92 hbarc included in line below ******************
 45      vcoul = (0.5/pi/pi)*(zalpha/q/q)*(ff(1)-cos(q*Rcut/hbarc))
 50      continue
      return
      end