code: def function(self, x, y, amp, R_sersic, Re, n_sersic, gamma, e1, e2, center_x=0, center_y=0, alpha=3.):
          
          phi_G, q = param_util.ellipticity2phi_q(e1, e2)
          Rb = R_sersic
          x_shift = x - center_x
          y_shift = y - center_y
  
          cos_phi = np.cos(phi_G)
          sin_phi = np.sin(phi_G)
  
          xt1 = cos_phi*x_shift+sin_phi*y_shift
          xt2 = -sin_phi*x_shift+cos_phi*y_shift
          xt2difq2 = xt2/(q*q)
          R_ = np.sqrt(xt1*xt1+xt2*xt2difq2)
          #R_ = R_.astype(np.float32)
          if isinstance(R_, int) or isinstance(R_, float):
              R_ = max(self._smoothing, R_)
          else:
              R_[R_ < self._smoothing] = self._smoothing
          if isinstance(R_, int) or isinstance(R_, float):
              R = max(self._smoothing, R_)
          else:
              R=np.empty_like(R_)
              _R = R_[R_ > self._smoothing]  #in the SIS regime
              R[R_ <= self._smoothing] = self._smoothing
              R[R_ > self._smoothing] = _R
  
          k, bn = self.k_bn(n_sersic, Re)
          result = amp * (1 + (Rb / R) ** alpha) ** (gamma / alpha) * np.exp(-bn * (((R ** alpha + Rb ** alpha) / Re ** alpha) ** (1. / (alpha * n_sersic)) - 1.))
          return np.nan_to_num(result)
docs: returns Core-Sersic function