code: def function(self, x, y, theta_E, gamma, e1, e2, center_x=0, center_y=0):
          
          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)
docs: :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