code: def _evaluate(self,R,z,phi=0.,t=0.):
          
          #Calculate relevant time
          if t < self._tform:
              smooth= 0.
          elif t < self._tsteady:
              deltat= t-self._tform
              xi= 2.*deltat/(self._tsteady-self._tform)-1.
              smooth= (3./16.*xi**5.-5./8*xi**3.+15./16.*xi+.5)
          else: #bar is fully on
              smooth= 1.
          r2= R**2.+z**2.
          r= numpy.sqrt(r2)
          if r <= self._rb:
              return self._af*smooth*numpy.cos(2.*(phi-self._omegab*t-self._barphi))\
                  *((r/self._rb)**3.-2.)*R**2./r2
          else:
              return -self._af*smooth*numpy.cos(2.*(phi-self._omegab*t-
                                                    self._barphi))\
                                                    *(self._rb/r)**3.\
                                                    *R**2./r2
docs: NAME:
             _evaluate
          PURPOSE:
             evaluate the potential at R,phi,t
          INPUT:
             R - Galactocentric cylindrical radius
             z - vertical height
             phi - azimuth
             t - time
          OUTPUT:
             Phi(R,z,phi,t)
          HISTORY:
             2010-11-24 - Started - Bovy (NYU)