code: def _dens(self,R,z,phi=0.,t=0.):
          
          if self.isNonAxi:
              R2= R**2.
              Rt2= R2*(1.-self._1m1overb2*nu.sin(phi)**2.)
              denom= 1./(Rt2+(z/self._q)**2.+self._core2)
              denom2= denom**2.
              return 1./4./nu.pi\
                  *(2.*Rt2/R2*(denom-Rt2*denom2)\
                        +denom/self._q**2.-2.*z**2.*denom**2./self._q**4.\
                        -self._1m1overb2\
                        *(2.*R2*nu.sin(2.*phi)**2./4.*self._1m1overb2\
                              *denom**2.+denom*nu.cos(2.*phi)))
          else:
              return 1./4./nu.pi/self._q**2.*((2.*self._q**2.+1.)*self._core2+R**2.\
                                                  +(2.-self._q**-2.)*z**2.)/\
                                                  (R**2.+(z/self._q)**2.+self._core2)**2.
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
          HISTORY: