docs: def gauss_fit(map_data, chs=None, mode='deg', amplitude=1, x_mean=0, y_mean=0, x_stddev=None, y_stddev=None, theta=None, cov_matrix=None, noise=0, **kwargs):
      """make a 2D Gaussian model and fit the observed data with the model.
  
      Args:
          map_data (xarray.Dataarray): Dataarray of cube or single chs.
          chs (list of int): in prep.
          mode (str): Coordinates for the fitting
              'pix'
              'deg'
          amplitude (float or None): Initial amplitude value of Gaussian fitting.
          x_mean (float): Initial value of mean of the fitting Gaussian in x.
          y_mean (float): Initial value of mean of the fitting Gaussian in y.
          x_stddev (float or None): Standard deviation of the Gaussian in x before rotating by theta.
          y_stddev  (float or None): Standard deviation of the Gaussian in y before rotating by theta.
          theta (float, optional or None): Rotation angle in radians.
          cov_matrix (ndarray, optional): A 2x2 covariance matrix. If specified, overrides the ``x_stddev``, ``y_stddev``, and ``theta`` defaults.
  
      Returns:
          decode cube (xarray cube) with fitting results in array and attrs.
      """
  
      if chs is None:
          chs = np.ogrid[0:63] # the number of channels would be changed
  
      if len(chs) > 1:
          for n,ch in enumerate(chs):
              subdata = np.transpose(np.full_like(map_data[:, :, ch], map_data.values[:, :, ch]))
              subdata[np.isnan(subdata)] = 0
  
              if mode == 'deg':
                  mX, mY = np.meshgrid(map_data.x, map_data.y)
  
              elif mode == 'pix':
                  mX, mY = np.mgrid[0:len(map_data.y), 0:len(map_data.x)]
  
              g_init = models.Gaussian2D(amplitude=np.nanmax(subdata),
                                         x_mean=x_mean, y_mean=y_mean,
                                         x_stddev=x_stddev, y_stddev=y_stddev, theta=theta,
                                         cov_matrix=cov_matrix, **kwargs) + models.Const2D(noise)
              fit_g  = fitting.LevMarLSQFitter()
              g      = fit_g(g_init, mX, mY, subdata)
  
              g_init2 = models.Gaussian2D(amplitude=np.nanmax(subdata-g.amplitude_1),
                                          x_mean=x_mean, y_mean=y_mean,
                                          x_stddev=x_stddev, y_stddev=y_stddev, theta=theta,
                                          cov_matrix=cov_matrix, **kwargs)
              fit_g2  = fitting.LevMarLSQFitter()
              g2      = fit_g2(g_init2, mX, mY, subdata)
  
              if n == 0:
                  results   = np.array([g2(mX,mY)])
                  peaks     = np.array([g2.amplitude.value])
                  x_means   = np.array([g2.x_mean.value])
                  y_means   = np.array([g2.y_mean.value])
                  x_stddevs = np.array([g2.x_stddev.value])
                  y_stddevs = np.array([g2.y_stddev.value])
                  thetas    = np.array([g2.theta.value])
                  if fit_g2.fit_info['param_cov'] is None:
                      unserts = nop.array([0])
                  else:
                      error   = np.diag(fit_g2.fit_info['param_cov'])**0.5
                      uncerts = np.array([error[0]])
  
              else:
                  results   = np.append(results, [g2(mX,mY)], axis=0)
                  peaks     = np.append(peaks, [g2.amplitude.value], axis=0)
                  x_means   = np.append(x_means, [g2.x_mean.value], axis=0)
                  y_means   = np.append(y_means, [g2.y_mean.value], axis=0)
                  x_stddevs = np.append(x_stddevs, [g2.x_stddev.value], axis=0)
                  y_stddevs = np.append(y_stddevs, [g2.y_stddev.value], axis=0)
                  thetas    = np.append(thetas, [g2.theta.value], axis=0)
  
                  if fit_g2.fit_info['param_cov'] is None:
                      uncerts = np.append(uncerts, [0], axis=0)
                  else:
                      error   = np.diag(fit_g2.fit_info['param_cov'])**0.5
                      uncerts = np.append(uncerts, [error[0]], axis=0)
  
          result = map_data.copy()
          result.values = np.transpose(results)
          result.attrs.update({'peak': peaks, 'x_mean': x_means, 'y_mean': y_means,
                               'x_stddev': x_stddevs, 'y_stddev': y_stddevs,
                               'theta': thetas, 'uncert': uncerts})
  
      else:
          subdata = np.transpose(np.full_like(map_data[:, :, 0], map_data.values[:, :, 0]))
          subdata[np.isnan(subdata)] = 0
  
          if mode == 'deg':
              mX, mY = np.meshgrid(map_data.x, map_data.y)
  
          elif mode == 'pix':
              mX, mY = np.mgrid[0:len(map_data.y), 0:len(map_data.x)]
  
          g_init = models.Gaussian2D(amplitude=np.nanmax(subdata),
                                     x_mean=x_mean, y_mean=y_mean,
                                     x_stddev=x_stddev, y_stddev=y_stddev, theta=theta,
                                     cov_matrix=cov_matrix, **kwargs) + models.Const2D(noise)
          fit_g  = fitting.LevMarLSQFitter()
          g      = fit_g(g_init, mX, mY, subdata)
  
          g_init2 = models.Gaussian2D(amplitude=np.nanmax(subdata-g.amplitude_1),
                                      x_mean=x_mean, y_mean=y_mean,
                                      x_stddev=x_stddev, y_stddev=y_stddev, theta=theta,
                                      cov_matrix=cov_matrix, **kwargs)
          fit_g2  = fitting.LevMarLSQFitter()
          g2      = fit_g2(g_init2, mX, mY, subdata)
  
          results   = np.array([g2(mX, mY)])
          peaks     = np.array([g2.amplitude.value])
          x_means   = np.array([g2.x_mean.value])
          y_means   = np.array([g2.y_mean.value])
          x_stddevs = np.array([g2.x_stddev.value])
          y_stddevs = np.array([g2.y_stddev.value])
          thetas    = np.array([g2.theta.value])
          error     = np.diag(fit_g2.fit_info['param_cov'])**0.5
          uncerts   = np.array(error[0])
  
          result = map_data.copy()
          result.values = np.transpose(results)
          result.attrs.update({'peak': peaks, 'x_mean': x_means, 'y_mean': y_means,
                               'x_stddev': x_stddevs, 'y_stddev': y_stddevs,
                               'theta': thetas, 'uncert': uncerts})
  
      return result