code: def bivnorm (sx, sy, cxy):
      
      _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
docs: 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