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``.

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