docs: def combine_psf(kernel_list_new, kernel_old, sigma_bkg, factor=1, stacking_option='median', symmetry=1):
          """
          updates psf estimate based on old kernel and several new estimates
          :param kernel_list_new: list of new PSF kernels estimated from the point sources in the image
          :param kernel_old: old PSF kernel
          :param sigma_bkg: estimated background noise in the image
          :param factor: weight of updated estimate based on new and old estimate, factor=1 means new estimate,
          factor=0 means old estimate
          :param stacking_option: option of stacking, mean or median
          :param symmetry: imposed symmetry of PSF estimate
          :return: updated PSF estimate and error_map associated with it
          """
  
          n = int(len(kernel_list_new) * symmetry)
          angle = 360. / symmetry
          kernelsize = len(kernel_old)
          kernel_list = np.zeros((n, kernelsize, kernelsize))
          i = 0
          for kernel_new in kernel_list_new:
              for k in range(symmetry):
                  kernel_rotated = image_util.rotateImage(kernel_new, angle * k)
                  kernel_norm = kernel_util.kernel_norm(kernel_rotated)
                  kernel_list[i, :, :] = kernel_norm
                  i += 1
  
          kernel_old_rotated = np.zeros((symmetry, kernelsize, kernelsize))
          for i in range(symmetry):
              kernel_old_rotated[i, :, :] = kernel_old
  
          kernel_list_new = np.append(kernel_list, kernel_old_rotated, axis=0)
          if stacking_option == 'median':
              kernel_new = np.median(kernel_list_new, axis=0)
          elif stacking_option == 'mean':
              kernel_new = np.mean(kernel_list_new, axis=0)
          else:
              raise ValueError(" stack_option must be 'median' or 'mean', %s is not supported." % stacking_option)
          kernel_new[kernel_new < 0] = 0
          kernel_new = kernel_util.kernel_norm(kernel_new)
          kernel_return = factor * kernel_new + (1.-factor)* kernel_old
  
          kernel_bkg = copy.deepcopy(kernel_return)
          kernel_bkg[kernel_bkg < sigma_bkg] = sigma_bkg
          error_map = np.var(kernel_list_new, axis=0) / kernel_bkg**2 / 2.
          return kernel_return, error_map