docs: 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)