code: def _dens(self,R,z,phi=0.,t=0.):
          
          x,y,z= self._compute_xyz(R,phi,z,t)
          m2 = x**2/self._a2+y**2/self._b2+z**2/self._c2
          if m2 < 1:
              return self._rhoc_M*(1.-m2/self.a**2)**self.n
          else:
              return 0.
docs: NAME:
             _dens
          PURPOSE:
             evaluate the density for this potential
          INPUT:
             R - Galactocentric cylindrical radius
             z - vertical height
             phi - azimuth
             t - time
          OUTPUT:
             the density