code: def _evaluate(self, R, z, phi=0, t=0):
          
          Ks = self._K(R)
          Bs = self._B(R)
          Ds = self._D(R)
  
          return -self._H * np.exp(-(R-self._r_ref) / self._Rs) \
                 * np.sum(self._Cs / Ks / Ds * np.cos(self._ns * self._gamma(R, phi - self._omega * t)) / np.cosh(Ks * z / Bs) ** Bs)
docs: NAME:
              _evaluate
          PURPOSE:
              Evaluate the potential at the given coordinates. (without the amp factor; handled by super class)
          INPUT:
              :param R: galactocentric cylindrical radius (must be scalar, not array)
              :param z: vertical height (must be scalar, not array)
              :param phi: azimuth (must be scalar, not array)
              :param t: time (must be scalar, not array)
          OUTPUT:
              :return: Phi(R, z, phi, t)
          HISTORY:
              2017-05-12  Jack Hong (UBC)