code: def bivconvolve (sx_a, sy_a, cxy_a, sx_b, sy_b, cxy_b):
      
      _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)
docs: 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.