def function(self, x, y, theta_E, gamma, e1, e2, center_x=0, center_y=0):
        """
        :param x: set of x-coordinates
        :type x: array of size (n)
        :param theta_E: Einstein radius of lense
        :type theta_E: float.
        :param gamma: power law slope of mass profifle
        :type gamma: <2 float
        :param q: Axis ratio
        :type q: 0<q<1
        :param phi_G: position angel of SES
        :type q: 0<phi_G<pi/2
        :returns:  function
        :raises: AttributeError, KeyError
        """
        phi_G, q = param_util.ellipticity2phi_q(e1, e2)
        gamma, q = self._param_bounds(gamma, q)
        theta_E *= q
        x_shift = x - center_x
        y_shift = y - center_y
        E = theta_E / (((3 - gamma) / 2.) ** (1. / (1 - gamma)) * np.sqrt(q))
        #E = phi_E
        eta = -gamma+3
        xt1 = np.cos(phi_G)*x_shift+np.sin(phi_G)*y_shift
        xt2 = -np.sin(phi_G)*x_shift+np.cos(phi_G)*y_shift
        p2 = xt1**2+xt2**2/q**2
        s2 = 0. # softening
        return 2 * E**2/eta**2 * ((p2 + s2)/E**2)**(eta/2)