def random_cov(d, diff=None):
    """Generate random covariance matrix.
    
    Generates a random covariance matrix, or two dependent covariance matrices
    if the argument `diff` is given.
    
    """
    S = 0.8*np.random.randn(d,d)
    copy_triu_to_tril(S)
    np.fill_diagonal(S,0)
    mineig = linalg.eigvalsh(S, eigvals=(0,0))[0]
    drand = 0.8*np.random.randn(d)
    if mineig < 0:
        S += np.diag(np.exp(drand)-mineig)
    else:
        S += np.diag(np.exp(drand))
    if not diff:
        return S.T
    S2 = S * np.random.randint(2, size=(d,d))*np.exp(diff*np.random.randn(d,d))
    copy_triu_to_tril(S2)
    np.fill_diagonal(S2,0)
    mineig = linalg.eigvalsh(S2, eigvals=(0,0))[0]
    drand += diff*np.random.randn(d)
    if mineig < 0:
        S2 += np.diag(np.exp(drand)-mineig)
    else:
        S2 += np.diag(np.exp(drand))
    return S.T, S2.T

def _validate_covars(covars, covariance_type, n_components):
    """Do basic checks on matrix covariance sizes and values
    """
    from scipy import linalg

    if covariance_type == "spherical":
        if len(covars) != n_components:
            raise ValueError("'spherical' covars have length n_components")
        elif np.any(covars <= 0):
            raise ValueError("'spherical' covars must be non-negative")
    elif covariance_type == "tied":
        if covars.shape[0] != covars.shape[1]:
            raise ValueError("'tied' covars must have shape (n_dim, n_dim)")
        elif not np.allclose(covars, covars.T) or np.any(linalg.eigvalsh(covars) <= 0):
            raise ValueError("'tied' covars must be symmetric, " "positive-definite")
    elif covariance_type == "diag":
        if len(covars.shape) != 2:
            raise ValueError("'diag' covars must have shape" "(n_components, n_dim)")
        elif np.any(covars <= 0):
            raise ValueError("'diag' covars must be non-negative")
    elif covariance_type == "full":
        if len(covars.shape) != 3:
            raise ValueError("'full' covars must have shape " "(n_components, n_dim, n_dim)")
        elif covars.shape[1] != covars.shape[2]:
            raise ValueError("'full' covars must have shape " "(n_components, n_dim, n_dim)")
        for n, cv in enumerate(covars):
            if not np.allclose(cv, cv.T) or np.any(linalg.eigvalsh(cv) <= 0):
                raise ValueError("component %d of 'full' covars must be " "symmetric, positive-definite" % n)
    else:
        raise ValueError("covariance_type must be one of " + "'spherical', 'tied', 'diag', 'full'")

 def _compute_eigen(self, eigenvalues_only=False, b=None):
     if eigenvalues_only:
         if b is None:
             if self.__e_val is not None: return self.__e_val
             else:
                 self.__e_val = eigvalsh(self._array)
                 return self.__e_val
         else:
             if self.__e_val_b.has_key(hash(str(b))):
                 return self.__e_val_b[hash(str(b))]
             else:
                 self.__e_val_b[hash(str(b))]=eigvalsh(self._array, b)
                 return self.__e_val_b[hash(str(b))]
     else:
         if b is None:
             if self.__e_val is not None and self.__e_vec is not None:
                 return self.__e_val, self.__e_vec
             else:
                 self.__e_val, self.__e_vec = eigh(self._array)
                 return self.__e_val, self.__e_vec
         else:
             if (self.__e_val_b.has_key(hash(str(b))) and
                 self.__e_vec_b.has_key(hash(str(b)))):
                 return self.__e_val_b[hash(str(b))], self.__e_vec_b[hash(str(b))]
             else:
                 self.__e_val_b[hash(str(b))], self.__e_vec_b[hash(str(b))] \
                    = eigvalsh(self._array, b)
                 self.__e_vec_t_b = self.__e_vec_b[hash(str(b))].T
                 return self.__e_val_b[hash(str(b))], self.__e_vec_b[hash(str(b))]

    def test1(self,n,p,prec,matrix_type,method):
        '''
        check the quality of tridiagonalization.
        '''
        if p==2:
            assert(method=='qr' or method=='sqrtm')
        else:
            assert(method=='qr')
        if prec is not None:
            gmpy2.get_context().precision=prec
        assert(matrix_type in ['array','sparse','mpc'])

        v0=self.gen_randv0(n,p,matrix_type='array')
        H0=self.gen_randomH(n,p,matrix_type=matrix_type)
        if p is not None:
            if method=='qr':
                data,offset=tridiagonalize_qr(H0,q=v0,prec=prec)
            else:
                data,offset=tridiagonalize2(H0,q=v0,prec=prec)
        else:
            data,offset=tridiagonalize(H0,q=v0,prec=prec)
        B=construct_tridmat(data,offset).toarray()
        if sps.issparse(H0):
            H0=H0.toarray()
        H0=complex128(H0)
        e1=eigvalsh(H0)
        e2=eigvalsh(B)
        assert_allclose(e1,e2)

def _validate_covars(covars, cvtype, nmix, n_dim):
    from scipy import linalg
    if cvtype == 'spherical':
        if len(covars) != nmix:
            raise ValueError("'spherical' covars must have length nmix")
        elif np.any(covars <= 0):
            raise ValueError("'spherical' covars must be non-negative")
    elif cvtype == 'tied':
        if covars.shape != (n_dim, n_dim):
            raise ValueError("'tied' covars must have shape (n_dim, n_dim)")
        elif (not np.allclose(covars, covars.T)
              or np.any(linalg.eigvalsh(covars) <= 0)):
            raise ValueError("'tied' covars must be symmetric, "
                             "positive-definite")
    elif cvtype == 'diag':
        if covars.shape != (nmix, n_dim):
            raise ValueError("'diag' covars must have shape (nmix, n_dim)")
        elif np.any(covars <= 0):
            raise ValueError("'diag' covars must be non-negative")
    elif cvtype == 'full':
        if covars.shape != (nmix, n_dim, n_dim):
            raise ValueError("'full' covars must have shape "
                             "(nmix, n_dim, n_dim)")
        for n, cv in enumerate(covars):
            if (not np.allclose(cv, cv.T)
                or np.any(linalg.eigvalsh(cv) <= 0)):
                raise ValueError("component %d of 'full' covars must be "
                                 "symmetric, positive-definite" % n)

def _dense_eigs(data, isherm, vecs, N, eigvals, num_large, num_small):
    """
    Internal functions for computing eigenvalues and eigenstates for a dense
    matrix.
    """
    if debug:
        logger.debug(inspect.stack()[0][3] + ": vectors = " + str(vecs))

    evecs = None

    if vecs:
        if isherm:
            if eigvals == 0:
                evals, evecs = la.eigh(data)
            else:
                if num_small > 0:
                    evals, evecs = la.eigh(
                        data, eigvals=[0, num_small - 1])
                if num_large > 0:
                    evals, evecs = la.eigh(
                        data, eigvals=[N - num_large, N - 1])
        else:
            evals, evecs = la.eig(data)
    else:
        if isherm:
            if eigvals == 0:
                evals = la.eigvalsh(data)
            else:
                if num_small > 0:
                    evals = la.eigvalsh(data, eigvals=[0, num_small - 1])
                if num_large > 0:
                    evals = la.eigvalsh(data, eigvals=[N - num_large, N - 1])
        else:
            evals = la.eigvals(data)

    _zipped = list(zip(evals, range(len(evals))))
    _zipped.sort()
    evals, perm = list(zip(*_zipped))

    if vecs:
        evecs = np.array([evecs[:, k] for k in perm])

    if not isherm and eigvals > 0:
        if vecs:
            if num_small > 0:
                evals, evecs = evals[:num_small], evecs[:num_small]
            elif num_large > 0:
                evals, evecs = evals[(N - num_large):], evecs[(N - num_large):]
        else:
            if num_small > 0:
                evals = evals[:num_small]
            elif num_large > 0:
                evals = evals[(N - num_large):]

    return np.array(evals), np.array(evecs)

def electronicEnergy(bonds, pa):
    """Compute the electronic contribution to the ground state energy of a
    polyene.

    Input arguments:
        bonds: a tuple of length equal to the number of carbon atoms in the
            polyene, containing the lengths of the bonds measured relative
            to their equilibrium lengths (i.e., bond[i] = u_{i+2} - u_{i+1}).
        pa: a dictionary containing the values of the system parameters.  These
        parameters must include,
            t0: the "bare" hopping (at equilibrium bond length), in eV
            alpha: the electron-phonon coupling, in units of eV per Angstrom
            chainlength: the number of carbon atoms in the polyene
            boundary: the type of boundary conditions ('periodic' or 'open')

    Returns: the ground state electronic energy of the specified polyene.
            
    Note that the length of the bonds tuple should be equal to the number of
    atoms even in the 'open' boundary condition case, in which the number of
    bonds in the molecule is one less than the number of atoms.  In this last
    case, the last entry of the bonds tuple should be set to 0."""
    
    if not set(['t0','alpha','chainlength','boundary']).issubset(set(pa.keys())):
        raise MissingParameters(set(pa.keys()))
    if len(bonds) != pa['chainlength']:
        raise BondNumberNotEqualChainLength(len(bonds), pa['chainlength'])
    if pa['boundary'] is 'open' and bonds[-1] != 0:
        raise OpenBCImplyLastBondAbsent(bonds[-1])
    energies = eigvalsh(hamiltonianMatrix(bonds, pa))
    return groundStateEnergy(energies)

def dos2d_ewindow(h,energies=np.linspace(-1.,1.,30),delta=None,info=False,
                    use_green=True,nk=300,mode="adaptive"):
  """Calculate the density of states in certain eenrgy window"""
  ys = [] # density of states
  if delta is None: # pick a good delta value
    delta = 0.1*(max(energies) - min(energies))/len(energies)
  if use_green:
    from green import bloch_selfenergy
    for energy in energies:
      (g,selfe) = bloch_selfenergy(h,nk=nk,energy=energy, delta=delta,
                   mode=mode)
      ys.append(-g.trace()[0,0].imag)
      if info: print("Done",energy)
    write_dos(energies,ys) # write in file
    return
  else: # do not use green function    
    from dosf90 import calculate_dos # import fortran library
    import scipy.linalg as lg
    kxs = np.linspace(0.,1.,nk)
    kys = np.linspace(0.,1.,nk)
    hkgen= h.get_hk_gen() # get hamiltonian generator
    ys = energies*0.
    weight = 1./(nk*nk)
    for ix in kxs:
      for iy in kys:
        k = np.array([ix,iy,0.]) # create kpoint
        hk = hkgen(k) # get hk hamiltonian
        evals = lg.eigvalsh(hk) # get eigenvalues
        ys += weight*calculate_dos(evals,energies,delta) # add this contribution
      if info: print("Done",ix)
    write_dos(energies,ys) # write in file
    return

def laplacian_spectrum(G, weight='weight'):
    """Return eigenvalues of the Laplacian of G

    Parameters
    ----------
    G : graph
       A NetworkX graph

    weight : string or None, optional (default='weight')
       The edge data key used to compute each value in the matrix.
       If None, then each edge has weight 1.

    Returns
    -------
    evals : NumPy array
      Eigenvalues

    Notes
    -----
    For MultiGraph/MultiDiGraph, the edges weights are summed.
    See to_numpy_matrix for other options.

    See Also
    --------
    laplacian_matrix
    """
    from scipy.linalg import eigvalsh
    return eigvalsh(nx.laplacian_matrix(G, weight=weight).todense())

def boolean_fermi_surface(h,write=True,output_file="BOOL_FERMI_MAP.OUT",
                    e=0.0,nk=50,nsuper=1,reciprocal=False,
                    delta=None):
  """Calculates the Fermi surface of a 2d system"""
  if h.dimensionality!=2: raise  # continue if two dimensional
  hk_gen = h.get_hk_gen() # gets the function to generate h(k)
  kxs = np.linspace(-nsuper,nsuper,nk)  # generate kx
  kys = np.linspace(-nsuper,nsuper,nk)  # generate ky
  kdos = [] # empty list
  kxout = []
  kyout = []
  if reciprocal: R = h.geometry.get_k2K() # get matrix
  # setup a reasonable value for delta
  if delta is None:
    delta = 8./np.max(np.abs(h.intra))/nk
  for x in kxs:
    for y in kxs:
      r = np.matrix([x,y,0.]).T # real space vectors
      k = np.array((R*r).T)[0] # change of basis
      hk = hk_gen(k) # get hamiltonian
      evals = lg.eigvalsh(hk) # diagonalize
      de = np.abs(evals - e) # difference with respect to fermi
      de = de[de<delta] # energies close to fermi
      if len(de)>0: kdos.append(1.0) # add to the list
      else: kdos.append(0.0) # add to the list
      kxout.append(x)
      kyout.append(y)
  if write:  # optionally, write in file
    f = open(output_file,"w") 
    for (x,y,d) in zip(kxout,kyout,kdos):
      f.write(str(x)+ "   "+str(y)+"   "+str(d)+"\n")
    f.close() # close the file
  return (kxout,kyout,d) # return result

def quad_form(x, P):
    x,P = map(Expression.cast_to_const, (x,P))
    # Check dimensions.
    n = P.size[0]
    if P.size[1] != n or x.size != (n,1):
        raise Exception("Invalid dimensions for arguments.")
    if x.curvature.is_constant():
        return x.T*P*x
    elif P.curvature.is_constant():
        np_intf = intf.get_matrix_interface(np.ndarray)
        P = np_intf.const_to_matrix(P.value)
        # Replace P with symmetric version.
        P = (P + P.T)/2
        # Check if P is PSD.
        eigvals = LA.eigvalsh(P)
        if min(eigvals) > 0:
            P_sqrt = Constant(LA.sqrtm(P).real)
            return square(norm2(P_sqrt*x))
        elif max(eigvals) < 0:
            P_sqrt = Constant(LA.sqrtm(-P).real)
            return -square(norm2(P_sqrt*x))
        else:
            raise Exception("P has both positive and negative eigenvalues.")
    else:
        raise Exception("At least one argument to quad_form must be constant.")

    def unweighted_likelihood_(self, x):
        """
        return the likelihood of each data for each component
        the values are not weighted by the component weights

        Parameters
        ----------
        x: array of shape (n_samples,self.dim)
           the data used in the estimation process

        Returns
        -------
        like, array of shape(n_samples,self.k)
          unweighted component-wise likelihood
        """
        n = x.shape[0]
        like = np.zeros((n, self.k))

        for k in range(self.k):
            # compute the data-independent factor first
            w = - np.log(2 * np.pi) * self.dim
            m = np.reshape(self.means[k], (1, self.dim))
            b = self.precisions[k]
            if self.prec_type == 'full':
                w += np.log(eigvalsh(b)).sum()
                dx = m - x
                q = np.sum(np.dot(dx, b) * dx, 1)
            else:
                w += np.sum(np.log(b))
                q = np.dot((m - x) ** 2, b)
            w -= q
            w /= 2
            like[:, k] = np.exp(w)
        return like

def _h_gen_roots_and_weights(n, mu, factor, func):
    """Compute the roots and weights for Gaussian-Hermite quadrature.
    Internal function.
    """
    if n < 1:
        raise ValueError("n must be positive.")
    
    bn = np.sqrt(np.arange(1,n, dtype=np.float64)/factor)
    c = np.diag(bn, -1)
    x = linalg.eigvalsh(c, overwrite_a=True)

    # improve roots by one application of Newton's method
    dy = func(n, x)
    df = factor*n*func(n-1, x)
    x -= dy/df

    df /= df.max()
    w = 1 / (df * df)

    # symmetrize
    w = (w + w[::-1])/2
    x = (x - x[::-1])/2

    # scale w correctly
    w *= np.sqrt(2.0*np.pi/factor) / w.sum()

    if mu:
        return [x, w, mu]
    else:
        return x, w

    def evaluate_eigenvalues_at(self, nodes, component=None, as_matrix=False):
        r"""
        Evaluate the eigenvalues :math:`\lambda_i\left(x\right)` at some grid nodes :math:`\gamma`.

        :param nodes: The grid nodes :math:`\gamma` we want to evaluate the eigenvalues at.
        :param component: The index :math:`i` of the eigenvalue :math:`\lambda_i` that gets evaluated.
        :param as_matrix: Returns the whole matrix :math:`\Lambda` instead of only a list with the eigenvalues :math:`\lambda_i`.
        :return: A sorted list with :math:`N` entries for all the eigenvalues evaluated at the nodes. Or a single value if a component was specified.
        """
        result = []

        # Hack to see if we evaluate at a single value
        if type(nodes) == numpy.ndarray:
            # Max to get rid of singular dimensions
            n = max(nodes.shape)
        else:
            try:
                n = len(nodes)
            except TypeError:
                n = len([nodes])

        # Memory for storing temporary values
        tmppot = numpy.ndarray((n,self.number_components,self.number_components), dtype=numpy.floating)
        tmpew = numpy.ndarray((n,self.number_components), dtype=numpy.floating)

        # evaluate potential
        values = self.evaluate_at(nodes)

        # fill in values
        for row in xrange(0, self.number_components):
            for col in xrange(0, self.number_components):
                tmppot[:, row, col] = values[row*self.number_components + col]

        # calculate eigenvalues assuming hermitian matrix (eigvalsh for stability!)
        for i in xrange(0, n):
            ew = linalg.eigvalsh(tmppot[i,:,:])
            # Sorting the eigenvalues biggest first.
            ew.sort()
            tmpew[i,:] = ew[::-1]

        tmp = [ tmpew[:,index] for index in xrange(0, self.number_components) ]

        if component is not None:
            (row, col) = component
            if row == col:
                result = tmp[row]
            else:
                result = numpy.zeros(tmp[row].shape, dtype=numpy.floating)
        elif as_matrix is True:
            result = []
            for row in xrange(self.number_components):
                for col in xrange(self.number_components):
                    if row == col:
                        result.append(tmp[row])
                    else:
                        result.append( numpy.zeros(tmp[row].shape, dtype=numpy.floating) )
        else:
            result = tmp

        return result

def disp_rel(kx_vec, ky_vec, m1=1):
    m, n = kx_vec.shape
    omega = np.zeros((4, m, n))
    for i in range(m):
        for j in range(n):
            kx = kx_vec[i,j]
            ky = ky_vec[i,j]

            km = 0.5*(sqrt(3)*ky - kx)
            kp = 0.5*(sqrt(3)*ky + kx)
            skx = sin(kx)
            skp = sin(kp)
            skm = sin(km)
            ckx = cos(kx)
            ckp = cos(kp)
            ckm = cos(km)
            sq3 = sqrt(3)

            K = np.array([
                [4, 0, 2*1j*skx - 1j*skp + 1j*skm - 2*ckx - ckp - ckm, 0],
                [0, 2*sq3, 0, -sq3*1j*skp + sq3*1j*skm - sq3*ckp - sq3*ckm],
                [-2*1j*skx + 1j*skp-1j*skm - 2*ckx - ckp - ckm, 0, 4, 0],
                [0, sq3*1j*skp - sq3*1j*skm - sq3*ckp - sq3*ckm, 0, 2*sq3]])



            M = np.array([[1., 0., 0., 0],
                          [0., 1., 0., 0],
                          [0., 0., m1, 0],
                          [0., 0., 0., m1]])
            vals = LA.eigvalsh(K, M)
            omega[:,i,j] = vals

    return omega

def eigbh(cm,bm,return_vecs=True):
    '''
    Get the eigenvalues and eigenvectors for matrice with block structure specified by block marker.

    Parameters:
        :cm: csr_matrix/bsr_matrix, the input matrix.
        :return_vecs: bool, return the eigenvectors or not.
        :bm: <BlockMarkerBase>, the block marker.

    Return:
        (eigenvalues,eigenvectors) if return vecs==True.
        (eigenvalues) if return vecs==False.
    '''
    EL,UL=[],[]
    is_sparse=sps.issparse(cm)
    for i in xrange(bm.nblock):
        mi=bm.extract_block(cm,(i,i))
        if return_vecs:
            if is_sparse: mi=mi.toarray()
            ei,ui=eigh(mi)
            EL.append(ei)
            UL.append(ui)
        else:
            ei=eigvalsh(mi)
            EL.append(ei)
    if return_vecs:
        return concatenate(EL),sps.block_diag(UL)
    else:
        return concatenate(EL)

    def test_nonint(self):
        #get the exact solution.
        h_exact=self.model_exact.hgen.H()
        E_excit=eigvalsh(h_exact)
        Emin_exact=sum(E_excit[E_excit<0])

        #the solution in occupation representation.
        h_occ=self.model_occ.hgen.H()
        Emin=eigsh(h_occ,which='SA',k=1)[0]
        print 'The Ground State Energy for hexagon(t = %s, t2 = %s) is %s, tolerence %s.'%(self.t,self.t2,Emin,Emin-Emin_exact)
        assert_almost_equal(Emin_exact,Emin)

        #the solution through updates
        H_serial=op2collection(op=self.model_occ.hgen.get_opH())
        H=get_H(H=H_serial,hgen=self.expander)
        H2,bm2=get_H_bm(H=H_serial,hgen=self.expander2,bstr='QM')
        Emin=eigsh(H,k=1,which='SA')[0]
        Emin2=eigsh(H2,k=1,which='SA')[0]
        print 'The Ground State Energy is %s, tolerence %s.'%(Emin,Emin-Emin2)
        assert_almost_equal(Emin_exact,Emin)
        assert_almost_equal(Emin_exact,Emin2)

        #the solution through dmrg.
        bmgen=get_bmgen(self.expander3.spaceconfig,'QM')
        dmrgegn=DMRGEngine(hchain=H_serial,hgen=self.expander3,tol=0,bmg=bmgen,symmetric=True)
        EG2=dmrgegn.run_finite(endpoint=(5,'<-',0),maxN=[10,20,30,40,40],tol=0)[-1]
        assert_almost_equal(Emin_exact,EG2*H_serial.nsite,decimal=4)

def _gauss_lobatto_mesh(n):
    from numpy import arange, asarray, diag, sqrt, zeros
    from scipy.linalg import eigvalsh

    if n < 2:
        raise ValueError("n must be > 1")
    if n == 2:
        return asarray((-1., 1.)), asarray((1., 1.))

    xi = zeros(n)
    wi = zeros(n)
    Pn = zeros(n)
    i = arange(1, n - 2)
    b = sqrt((i * (2. + i)) / (3. + 4.*i * (2. + i)))  # coeff for Jacobi Poly with a=b=1

    M = diag(b, -1) + diag(b, 1)
    xi[1:n - 1] = eigvalsh(M)
    xi[0] = -1.; xi[-1] = 1.

    Pim2 = 1.  # P_{i-2}
    Pim1 = xi  # P_{i-1}
    for i in range(2, n):  # want P_{n-1}
        wi = (1. / i) * ((2 * i - 1) * xi * Pim1 - (i - 1) * Pim2)
        Pim2 = Pim1
        Pim1 = wi
    wi = 2. / (n * (n - 1) * wi ** 2)
    wi[0] = wi[-1] = 2. / (n * (n - 1))
    return xi, wi

def check_dfunc(dfunc,nband,D,Gap=0.,method='eval'):
    '''
    Checking for dfunc.

    dfunc:
        the hybridization function.
    nband:
        the number of bands.
    method:
        the method for checking.
        * `pauli` -> pauli decomposition for 2D.
        * `eval` -> check for eigenvalue.
    '''
    ion()
    if nband!=2 and method=='pauli':
        warnings.warn('Checking pauli components is not allowed for non-2band bath!')
        method='eval'
    if ndim(D)==0:
        D=[-D,D]
    wlist=linspace(D[0],D[1],1000)
    if nband==1:
        dl=array([dfunc(w) for w in wlist])
    elif method=='eval':
        dl=array([eigvalsh(dfunc(w)) for w in wlist])
    elif method=='pauli':
        dl=array([s2vec(dfunc(w)) for w in wlist])
    plot(wlist,dl)
    if method=='eval':
        legend([r'$\rho_%s$'%i for i in xrange(nband)])
    elif method=='pauli':
        legend([r'$\rho_0$',r'$\rho_x$',r'$\rho_y$',r'$\rho_z$'])
    pdb.set_trace()

def plot_mohr3d(S):
    r"""Plot 3D Mohr circles."""
    
    S3, S2, S1 = eigvalsh(S)
 
    R_maj = 0.5*(S1 - S3)
    cent_maj = 0.5*(S1+S3)
    
    R_min = 0.5*(S2 - S3)
    cent_min = 0.5*(S2 + S3)
    
    R_mid = 0.5*(S1 - S2)
    cent_mid = 0.5*(S1 + S2)
    
    circ1 = plt.Circle((cent_maj,0), R_maj, facecolor='#cce885', lw=3,
                       edgecolor='#5c8037')
    circ2 = plt.Circle((cent_min,0), R_min, facecolor='w', lw=3,
                       edgecolor='#15a1bd')
    circ3 = plt.Circle((cent_mid,0), R_mid, facecolor='w', lw=3,
                       edgecolor='#e4612d')
    plt.axis('image')
    ax = plt.gca()
    ax.add_artist(circ1)
    ax.add_artist(circ2)
    ax.add_artist(circ3)
    ax.set_xlim(S3 - .1*R_maj, S1 + .1*R_maj)
    ax.set_ylim(-1.1*R_maj, 1.1*R_maj)
    plt.xlabel(r"$\sigma$", size=18)
    plt.ylabel(r"$\tau$", size=18)
    #plt.savefig('Mohr_circle_3D.svg')
    plt.show()