docs: def daofind_marginal_fit(self, axis=0):
          """
          Fit 1D Gaussians, defined from the marginal x/y kernel
          distributions, to the marginal x/y distributions of the original
          (unconvolved) image.
  
          These fits are used calculate the star centroid and roundness
          ("GROUND") properties.
  
          Parameters
          ----------
          axis : {0, 1}, optional
              The axis for which the marginal fit is performed:
  
              * 0: for the x axis
              * 1: for the y axis
  
          Returns
          -------
          dx : float
              The fractional shift in x or y (depending on ``axis`` value)
              of the image centroid relative to the maximum pixel.
  
          hx : float
              The height of the best-fitting Gaussian to the marginal x or
              y (depending on ``axis`` value) distribution of the
              unconvolved source data.
          """
  
          # define triangular weighting functions along each axis, peaked
          # in the middle and equal to one at the edge
          x = self.xcenter - np.abs(np.arange(self.nx) - self.xcenter) + 1
          y = self.ycenter - np.abs(np.arange(self.ny) - self.ycenter) + 1
          xwt, ywt = np.meshgrid(x, y)
  
          if axis == 0:    # marginal distributions along x axis
              wt = xwt[0]    # 1D
              wts = ywt    # 2D
              size = self.nx
              center = self.xcenter
              sigma = self.kernel.xsigma
              dxx = center - np.arange(size)
          elif axis == 1:    # marginal distributions along y axis
              wt = np.transpose(ywt)[0]    # 1D
              wts = xwt    # 2D
              size = self.ny
              center = self.ycenter
              sigma = self.kernel.ysigma
              dxx = np.arange(size) - center
  
          # compute marginal sums for given axis
          wt_sum = np.sum(wt)
          dx = center - np.arange(size)
  
          # weighted marginal sums
          kern_sum_1d = np.sum(self.kernel.gaussian_kernel_unmasked * wts,
                               axis=axis)
          kern_sum = np.sum(kern_sum_1d * wt)
          kern2_sum = np.sum(kern_sum_1d**2 * wt)
  
          dkern_dx = kern_sum_1d * dx
          dkern_dx_sum = np.sum(dkern_dx * wt)
          dkern_dx2_sum = np.sum(dkern_dx**2 * wt)
          kern_dkern_dx_sum = np.sum(kern_sum_1d * dkern_dx * wt)
  
          data_sum_1d = np.sum(self.data * wts, axis=axis)
          data_sum = np.sum(data_sum_1d * wt)
          data_kern_sum = np.sum(data_sum_1d * kern_sum_1d * wt)
          data_dkern_dx_sum = np.sum(data_sum_1d * dkern_dx * wt)
          data_dx_sum = np.sum(data_sum_1d * dxx * wt)
  
          # perform linear least-squares fit (where data = sky + hx*kernel)
          # to find the amplitude (hx)
          # reject the star if the fit amplitude is not positive
          hx_numer = data_kern_sum - (data_sum * kern_sum) / wt_sum
          if hx_numer <= 0.:
              return np.nan, np.nan
  
          hx_denom = kern2_sum - (kern_sum**2 / wt_sum)
          if hx_denom <= 0.:
              return np.nan, np.nan
  
          # compute fit amplitude
          hx = hx_numer / hx_denom
          # sky = (data_sum - (hx * kern_sum)) / wt_sum
  
          # compute centroid shift
          dx = ((kern_dkern_dx_sum -
                 (data_dkern_dx_sum - dkern_dx_sum*data_sum)) /
                (hx * dkern_dx2_sum / sigma**2))
  
          hsize = size / 2.
          if abs(dx) > hsize:
              if data_sum == 0.:
                  dx = 0.0
              else:
                  dx = data_dx_sum / data_sum
                  if abs(dx) > hsize:
                      dx = 0.0
  
          return dx, hx