docs: def gaussian_deconvolve (smaj, smin, spa, bmaj, bmin, bpa):
      """Deconvolve two Gaussians analytically.
  
      Given the shapes of 2-dimensional “source” and “beam” Gaussians, this
      returns a deconvolved “result” Gaussian such that the convolution of
      “beam” and “result” is “source”.
  
      Arguments:
  
      smaj
        Major axis of source Gaussian.
      smin
        Minor axis of source Gaussian.
      spa
        Orientation angle of source Gaussian, in radians.
      bmaj
        Major axis of beam Gaussian.
      bmin
        Minor axis of beam Gaussian.
      bpa
        Orientation angle of beam Gaussian, in radians.
  
      The return value is ``(rmaj, rmin, rpa, status)``. The first three values
      have the same format as the input arguments. The *status* result is one of
      "ok", "pointlike", or "fail". A "pointlike" status indicates that the
      source and beam shapes are difficult to distinguish; a "fail" status
      indicates that the two shapes seem to be mutually incompatible (e.g.,
      source and beam are very narrow and orthogonal).
  
      The axes can be measured in any units, so long as they're consistent.
  
      Ideally if::
  
        rmaj, rmin, rpa, status = gaussian_deconvolve (smaj, smin, spa, bmaj, bmin, bpa)
  
      then::
  
        smaj, smin, spa = gaussian_convolve (rmaj, rmin, rpa, bmaj, bmin, bpa)
  
      Implementation derived from MIRIAD’s ``gaudfac``. This function currently
      doesn't do a great job of dealing with pointlike sources, i.e. ones where
      “source” and “beam” are nearly indistinguishable.
  
      """
      # I've added extra code to ensure ``smaj >= bmaj``, ``smin >= bmin``, and
      # increased the coefficient in front of "limit" from 0.1 to 0.5. Feel a
      # little wary about that first change.
  
      from numpy import cos, sin, sqrt, min, abs, arctan2
  
      if smaj < bmaj:
          smaj = bmaj
      if smin < bmin:
          smin = bmin
  
      alpha = ((smaj * cos (spa))**2 + (smin * sin (spa))**2 -
               (bmaj * cos (bpa))**2 - (bmin * sin (bpa))**2)
      beta = ((smaj * sin (spa))**2 + (smin * cos (spa))**2 -
              (bmaj * sin (bpa))**2 - (bmin * cos (bpa))**2)
      gamma = 2 * ((smin**2 - smaj**2) * sin (spa) * cos (spa) -
                   (bmin**2 - bmaj**2) * sin (bpa) * cos (bpa))
  
      s = alpha + beta
      t = sqrt ((alpha - beta)**2 + gamma**2)
      limit = 0.5 * min ([smaj, smin, bmaj, bmin])**2
      status = 'ok'
  
      if alpha < 0 or beta < 0 or s < t:
          dmaj = dmin = dpa = 0
  
          if 0.5 * (s - t) < limit and alpha > -limit and beta > -limit:
              status = 'pointlike'
          else:
              status = 'fail'
      else:
          dmaj = sqrt (0.5 * (s + t))
          dmin = sqrt (0.5 * (s - t))
  
          if abs (gamma) + abs (alpha - beta) == 0:
              dpa = 0
          else:
              dpa = 0.5 * arctan2 (-gamma, alpha - beta)
  
      return dmaj, dmin, dpa, status