def METHOD_NAME(self, x, y, sigma0, Rs, e1, e2, center_x=0, center_y=0):
    """Returns df/dx and df/dy of the function (integral of NFW)"""
    phi_q, q = param_util.ellipticity2phi_q(e1, e2)
    # shift
    x_ = x - center_x
    y_ = y - center_y
    # rotate
    x__, y__ = util.rotate(x_, y_, phi_q)
    f__x, f__y = self.cse_major_axis_set.METHOD_NAME(
        x__ / Rs, y__ / Rs, self._a_list, self._s_list, q
    )
    # rotate deflections back
    f_x, f_y = util.rotate(f__x, f__y, -phi_q)
    const = self._normalization(sigma0, Rs, q) / Rs
    return const * f_x, const * f_y