docs: def gaussian_convolve (maj1, min1, pa1, maj2, min2, pa2):
      """Convolve two Gaussians analytically.
  
      Given the shapes of two 2-dimensional Gaussians, this function returns
      the shape of their convolution.
  
      Arguments:
  
      maj1
        Major axis of input Gaussian 1.
      min1
        Minor axis of input Gaussian 1.
      pa1
        Orientation angle of input Gaussian 1, in radians.
      maj2
        Major axis of input Gaussian 2.
      min2
        Minor axis of input Gaussian 2.
      pa2
        Orientation angle of input Gaussian 2, in radians.
  
      The return value is ``(maj3, min3, pa3)``, with the same format as the
      input arguments. The axes can be measured in any units, so long as they're
      consistent.
  
      Implementation copied from MIRIAD’s ``gaufac``.
  
      """
      c1 = np.cos (pa1)
      s1 = np.sin (pa1)
      c2 = np.cos (pa2)
      s2 = np.sin (pa2)
  
      a = (maj1*c1)**2 + (min1*s1)**2 + (maj2*c2)**2 + (min2*s2)**2
      b = (maj1*s1)**2 + (min1*c1)**2 + (maj2*s2)**2 + (min2*c2)**2
      g = 2 * ((min1**2 - maj1**2) * s1 * c1 + (min2**2 - maj2**2) * s2 * c2)
  
      s = a + b
      t = np.sqrt ((a - b)**2 + g**2)
      maj3 = np.sqrt (0.5 * (s + t))
      min3 = np.sqrt (0.5 * (s - t))
  
      if abs (g) + abs (a - b) == 0:
          pa3 = 0.
      else:
          pa3 = 0.5 * np.arctan2 (-g, a - b)
  
      # "Amplitude of the resulting Gaussian":
      # f = pi / (4 * np.log (2)) * maj1 * min1 * maj2 * min2 \
      #    / np.sqrt (a * b - 0.25 * g**2)
  
      return maj3, min3, pa3