Given the parameters of a Gaussian bivariate distribution, compute the
    correct normalization for the equivalent 2D Gaussian. It's 1 / (2 pi sqrt
    (sx**2 sy**2 - cxy**2). This function adds a lot of sanity checking.

    Inputs:

    * sx: standard deviation (not variance) of x var
    * sy: standard deviation (not variance) of y var
    * cxy: covariance (not correlation coefficient) of x and y

    Returns: the scalar normalization

def bivnorm (sx, sy, cxy):
    """Given the parameters of a Gaussian bivariate distribution, compute the
    correct normalization for the equivalent 2D Gaussian. It's 1 / (2 pi sqrt
    (sx**2 sy**2 - cxy**2). This function adds a lot of sanity checking.

    Inputs:

    * sx: standard deviation (not variance) of x var
    * sy: standard deviation (not variance) of y var
    * cxy: covariance (not correlation coefficient) of x and y

    Returns: the scalar normalization

    """
    _bivcheck (sx, sy, cxy)
    from numpy import pi, sqrt

    t = (sx * sy)**2 - cxy**2
    if t <= 0:
        raise ValueError ('covariance just barely out of bounds '
                          '(sx=%.10e, sy=%.10e, cxy=%.10e, cxy/sxsy=%.16f)' %
                          (sx, sy, cxy, cxy / (sx * sy)))
    return (2 * pi * sqrt (t))**-1