Given the parameters of a Gaussian bivariate distribution, compute the
    parameters for the equivalent 2D Gaussian in ellipse form (major, minor,
    pa).

    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

    Outputs:

    * mjr: major axis of equivalent 2D Gaussian (sigma, not FWHM)
    * mnr: minor axis
    * pa: position angle, rotating from +x to +y

    Lots of sanity-checking because it's obnoxiously easy to have numerics
    that just barely blow up on you.

def bivell (sx, sy, cxy):
    """Given the parameters of a Gaussian bivariate distribution, compute the
    parameters for the equivalent 2D Gaussian in ellipse form (major, minor,
    pa).

    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

    Outputs:

    * mjr: major axis of equivalent 2D Gaussian (sigma, not FWHM)
    * mnr: minor axis
    * pa: position angle, rotating from +x to +y

    Lots of sanity-checking because it's obnoxiously easy to have numerics
    that just barely blow up on you.

    """
    # See CfA notebook #1 pp. 129-133.
    _bivcheck (sx, sy, cxy)
    from numpy import arctan2, sqrt

    sx2, sy2, cxy2 = sx**2, sy**2, cxy**2

    pa = 0.5 * arctan2 (2 * cxy, sx2 - sy2)
    h = sqrt ((sx2 - sy2)**2 + 4*cxy2)

    t = 2 * (sx2 * sy2 - cxy2) / (sx2 + sy2 - h)
    if t < 0:
        raise ValueError ('covariance just barely out of bounds [1] '
                          '(sx=%.10e, sy=%.10e, cxy=%.10e, cxy/sxsy=%.16f)' %
                          (sx, sy, cxy, cxy / (sx * sy)))
    mjr = sqrt (t)

    t = 2 * (sx2 * sy2 - cxy2) / (sx2 + sy2 + h)
    if t < 0: # if we got this far, shouldn't happen, but ...
        raise ValueError ('covariance just barely out of bounds [2] '
                          '(sx=%.10e, sy=%.10e, cxy=%.10e, cxy/sxsy=%.16f)' %
                          (sx, sy, cxy, cxy / (sx * sy)))
    mnr = sqrt (t)

    return ellnorm (mjr, mnr, pa)