<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 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