Given two independent bivariate distributions, compute a bivariate
    distribution corresponding to their convolution.

    I'm sure this is worked out in a ton of places, but I got the equations
    from Pineau+ (2011A&A...527A.126P).

    Returns: (sx_c, sy_c, cxy_c), the parameters of the convolved
    distribution.

def bivconvolve (sx_a, sy_a, cxy_a, sx_b, sy_b, cxy_b):
    """Given two independent bivariate distributions, compute a bivariate
    distribution corresponding to their convolution.

    I'm sure this is worked out in a ton of places, but I got the equations
    from Pineau+ (2011A&A...527A.126P).

    Returns: (sx_c, sy_c, cxy_c), the parameters of the convolved
    distribution.

    """
    _bivcheck (sx_a, sy_a, cxy_a)
    _bivcheck (sx_b, sy_b, cxy_b)

    sx_c = np.sqrt (sx_a**2 + sx_b**2)
    sy_c = np.sqrt (sy_a**2 + sy_b**2)
    cxy_c = cxy_a + cxy_b

    return _bivcheck (sx_c, sy_c, cxy_c)