code: def ellipse_to_cov(sigma_maj, sigma_min, theta):
      
      cth = np.cos(theta)
      sth = np.sin(theta)
      covxx = cth**2 * sigma_maj**2 + sth**2 * sigma_min**2
      covyy = sth**2 * sigma_maj**2 + cth**2 * sigma_min**2
      covxy = cth * sth * sigma_maj**2 - cth * sth * sigma_min**2
      return np.array([[covxx, covxy], [covxy, covyy]])
docs: Compute the covariance matrix in two variables x and y given
      the std. deviation along the semi-major and semi-minor axes and
      the rotation angle of the error ellipse.
  
      Parameters
      ----------
      sigma_maj : float
          Std. deviation along major axis of error ellipse.
  
      sigma_min : float
          Std. deviation along minor axis of error ellipse.
  
      theta : float
          Rotation angle in radians from x-axis to ellipse major axis.