code: def function(self, x, y, amp, sigma, center_x=0, center_y=0):
          
          x_ = x - center_x
          y_ = y - center_y
          r = np.sqrt(x_**2 + y_**2)
          sigma_x, sigma_y = sigma, sigma
          c = 1. / (2 * sigma_x * sigma_y)
          num_int = self._num_integral(r, c)
          amp_density = self._amp2d_to_3d(amp, sigma_x, sigma_y)
          amp2d = amp_density / (np.sqrt(np.pi) * np.sqrt(sigma_x * sigma_y * 2))
          amp2d *= 2 * 1. / (2 * c)
          return num_int * amp2d
docs: returns Gaussian