Given a 2D Gaussian expressed as an ellipse (major, minor, pa), compute the
    equivalent parameters for a Gaussian bivariate distribution. We assume
    that the ellipse is normalized so that the functions evaluate identicall
    for major = minor.

    Inputs:

    * mjr: major axis (sigma not FWHM) of the Gaussian
    * mnr: minor axis (sigma not FWHM) of the Gaussian
    * pa: position angle (from +x to +y) of the Gaussian, radians

    Returns:

    * 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

def ellbiv (mjr, mnr, pa):
    """Given a 2D Gaussian expressed as an ellipse (major, minor, pa), compute the
    equivalent parameters for a Gaussian bivariate distribution. We assume
    that the ellipse is normalized so that the functions evaluate identicall
    for major = minor.

    Inputs:

    * mjr: major axis (sigma not FWHM) of the Gaussian
    * mnr: minor axis (sigma not FWHM) of the Gaussian
    * pa: position angle (from +x to +y) of the Gaussian, radians

    Returns:

    * 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

    """
    _ellcheck (mjr, mnr, pa)

    cpa, spa = np.cos (pa), np.sin (pa)
    q = np.asarray ([[cpa, -spa],
                     [spa, cpa]])
    cov = np.dot (q, np.dot (np.diag ([mjr**2, mnr**2]), q.T))
    sx = np.sqrt (cov[0,0])
    sy = np.sqrt (cov[1,1])
    cxy = cov[0,1]

    return _bivcheck (sx, sy, cxy)