code: def _dens(self,R,z,phi=0.,t=0.):
          
          if self.alpha == 0.:
              return 1./4./nu.pi/self.q2*((2.*self.q2+1.)*self.core2+R**2.\
                                         +(2.-1./self.q2)*z**2.)/\
                                         (R**2.+z**2./self.q2+self.core2)**2.
          else:
              m2= self.core2+R**2.+z**2./self.q2
              return 1./self.q2*(self.core2*(1.+2.*self.q2)+R**2.*(1.-self.alpha*self.q2)+z**2.*(2.-(1.+self.alpha)/self.q2))*m2**(-self.alpha/2.-2.)/4./nu.pi
docs: NAME:
             _dens
          PURPOSE:
             evaluate the density force for this potential
          INPUT:
             R - Galactocentric cylindrical radius
             z - vertical height
             phi - azimuth
             t - time
          OUTPUT:
             the density
          HISTORY:
             2013-01-09 - Written - Bovy (IAS)