def _spectrum_es(H, wlist, c_ops, a_op, b_op):
    """
    Internal function for calculating the spectrum of the correlation function
    :math:`\left<A(\\tau)B(0)\\right>`.
    """
    if debug:
        print(inspect.stack()[0][3])

    # construct the Liouvillian
    L = liouvillian(H, c_ops)

    # find the steady state density matrix and a_op and b_op expecation values
    rho0 = steadystate(L)

    a_op_ss = expect(a_op, rho0)
    b_op_ss = expect(b_op, rho0)

    # eseries solution for (b * rho0)(t)
    es = ode2es(L, b_op * rho0)

    # correlation
    corr_es = expect(a_op, es)

    # covariance
    cov_es = corr_es - np.real(np.conjugate(a_op_ss) * b_op_ss)

    # spectrum
    spectrum = esspec(cov_es, wlist)

    return spectrum

def test_ho_lgmres():
    "Steady state: Thermal HO - iterative-lgmres solver"
    # thermal steadystate of an oscillator: compare numerics with analytical
    # formula
    a = destroy(40)
    H = 0.5 * 2 * np.pi * a.dag() * a
    gamma1 = 0.05

    wth_vec = np.linspace(0.1, 3, 20)
    p_ss = np.zeros(np.shape(wth_vec))

    for idx, wth in enumerate(wth_vec):

        n_th = 1.0 / (np.exp(1.0 / wth) - 1)  # bath temperature
        c_op_list = []
        rate = gamma1 * (1 + n_th)
        c_op_list.append(np.sqrt(rate) * a)
        rate = gamma1 * n_th
        c_op_list.append(np.sqrt(rate) * a.dag())
        rho_ss = steadystate(H, c_op_list, method='iterative-lgmres')
        p_ss[idx] = np.real(expect(a.dag() * a, rho_ss))

    p_ss_analytic = 1.0 / (np.exp(1.0 / wth_vec) - 1)
    delta = sum(abs(p_ss_analytic - p_ss))
    assert_equal(delta < 1e-3, True)

def test_qubit_power():
    "Steady state: Thermal qubit - power solver"
    # thermal steadystate of a qubit: compare numerics with analytical formula
    sz = sigmaz()
    sm = destroy(2)

    H = 0.5 * 2 * np.pi * sz
    gamma1 = 0.05

    wth_vec = np.linspace(0.1, 3, 20)
    p_ss = np.zeros(np.shape(wth_vec))

    for idx, wth in enumerate(wth_vec):

        n_th = 1.0 / (np.exp(1.0 / wth) - 1)  # bath temperature
        c_op_list = []
        rate = gamma1 * (1 + n_th)
        c_op_list.append(np.sqrt(rate) * sm)
        rate = gamma1 * n_th
        c_op_list.append(np.sqrt(rate) * sm.dag())
        rho_ss = steadystate(H, c_op_list, method='power')
        p_ss[idx] = expect(sm.dag() * sm, rho_ss)

    p_ss_analytic = np.exp(-1.0 / wth_vec) / (1 + np.exp(-1.0 / wth_vec))
    delta = sum(abs(p_ss_analytic - p_ss))
    assert_equal(delta < 1e-5, True)

    def steady(N = 20):  # number of basis states to consider
        n=num(N)
        a = destroy(N)
        H = a.dag() * a
        print H.eigenstates()
        #psi0 = basis(N, 10)  # initial state
        kappa = 0.1  # coupling to oscillator
        c_op_list = []
        n_th_a = 2  # temperature with average of 2 excitations
        rate = kappa * (1 + n_th_a)
        c_op_list.append(sqrt(rate) * a)  # decay operators
        rate = kappa * n_th_a
        c_op_list.append(sqrt(rate) * a.dag())  # excitation operators
        final_state = steadystate(H, c_op_list)
        fexpt = expect(a.dag() * a, final_state)

        #tlist = linspace(0, 100, 100)

        #mcdata = mcsolve(H, psi0, tlist, c_op_list, [a.dag() * a], ntraj=100)

        #medata = mesolve(H, psi0, tlist, c_op_list, [a.dag() * a])
        #plot(tlist, mcdata.expect[0],
        #plt.plot(tlist, medata.expect[0], lw=2)
        plt.axhline(y=fexpt, color='r', lw=1.5) # ss expt. value as horiz line (= 2)
        plt.ylim([0, 10])
        plt.show()

    def find_expect(self, vg, pwr, fd):
        if self.power_plot:
            phi, pwr=vg
        else:
            phi, fd=vg
        pwr_fridge=pwr-self.atten
        lin_pwr=0.001*10**(pwr_fridge/10.0)
        Omega=sqrt(lin_pwr/h*2.0)

        gamma, Delta=self._get_GammaDelta(fd=fd, f0=self.f0, Np=self.Np, gamma=self.gamma)
        g_el=self._get_Gamma_C(fd=fd)
        wTvec=self._get_fTvec(phi=phi, gamma=gamma, Delta=Delta, fd=fd, Psaw=lin_pwr)

        if self.acoustic_plot:
            Om=Omega*sqrt(gamma/fd)
        else:
            Om=Omega*sqrt(g_el/fd)
        wT = wTvec-fd*self.nvec #rotating frame of gate drive \omega_m-m*\omega_\gate
        transmon_levels = Qobj(diag(wT[range(self.N_dim)]))
        rate1 = (gamma+g_el)*(1.0 + self.N_gamma)
        rate2 = (gamma+g_el)*self.N_gamma
        c_ops=[sqrt(rate1)*self.a_op, sqrt(rate2)*self.a_dag]#, sqrt(rate3)*self.a_op, sqrt(rate4)*self.a_dag]
        Omega_vec=-0.5j*(Om*self.a_dag - conj(Om)*self.a_op)
        H=transmon_levels +Omega_vec
        final_state = steadystate(H, c_ops) #solve master equation
        fexpt=expect(self.a_op, final_state) #expectation value of relaxation operator
        #return fexpt
        if self.acoustic_plot:
            return 1.0*gamma/Om*fexpt
        else:
            return 1.0*sqrt(g_el*gamma)/Om*fexpt

    def rhos(self, nslice=None):
        """rho
        return steadystate density matrices"""
        self.precalc = True

        if self.noisy:
            print("N = {}, len() = {}".format(self.N_field_levels, len(self.long_range)))

        def progress(*args):
            if self.noisy:
                print("|", sep="", end="", flush=True)
            return args

        if nslice is not None:
            return np.asarray(
                [progress(qt.steadystate(ham, self._c_ops())) for ham in list(self.hamiltonian())[nslice]]
            ).T
        else:
            return np.asarray([progress(qt.steadystate(ham, self._c_ops())) for ham in self.hamiltonian()]).T

def test_dqd_current():
    "Counting statistics: current and current noise in a DQD model"

    G = 0
    L = 1
    R = 2

    sz = projection(3, L, L) - projection(3, R, R)
    sx = projection(3, L, R) + projection(3, R, L)
    sR = projection(3, G, R)
    sL = projection(3, G, L)

    w0 = 1
    tc = 0.6 * w0
    GammaR = 0.0075 * w0
    GammaL = 0.0075 * w0
    nth = 0.00
    eps_vec = np.linspace(-1.5*w0, 1.5*w0, 20)

    J_ops = [GammaR * sprepost(sR, sR.dag())]

    c_ops = [np.sqrt(GammaR * (1 + nth)) * sR,
             np.sqrt(GammaR * (nth)) * sR.dag(),
             np.sqrt(GammaL * (nth)) * sL,
             np.sqrt(GammaL * (1 + nth)) * sL.dag(),
             ]

    I = np.zeros(len(eps_vec))
    S = np.zeros(len(eps_vec))

    for n, eps in enumerate(eps_vec):
        H = (eps/2 * sz + tc * sx)
        L = liouvillian(H, c_ops)
        rhoss = steadystate(L)
        I[n], S[n] = countstat_current_noise(L, [], rhoss=rhoss, J_ops=J_ops)

        I2 = countstat_current(L, rhoss=rhoss, J_ops=J_ops)
        assert_(abs(I[n] - I2) < 1e-8)

        I2 = countstat_current(L, c_ops, J_ops=J_ops)
        assert_(abs(I[n] - I2) < 1e-8)

    Iref = tc**2 * GammaR / (tc**2 * (2 + GammaR/GammaL) +
                             GammaR**2/4 + eps_vec**2)
    Sref = 1 * Iref * (
        1 - 8 * GammaL * tc**2 *
        (4 * eps_vec**2 * (GammaR - GammaL) +
         GammaR * (3 * GammaL * GammaR + GammaR**2 + 8*tc**2)) /
        (4 * tc**2 * (2 * GammaL + GammaR) + GammaL * GammaR**2 +
         4 * eps_vec**2 * GammaL)**2
    )

    assert_allclose(I, Iref, 1e-4)
    assert_allclose(S, Sref, 1e-4)

def find_expect(vg): #phi=0.1, Omega_vec=3.0):
    phi, Omega=vg#.shape
    Omega_vec=- 0.5j*(Omega*adag - conj(Omega)*a)

    Ej = Ejmax*absolute(cos(pi*phi)) #Josephson energy as function of Phi.

    wTvec = -Ej + sqrt(8.0*Ej*Ec)*(nvec+0.5)+Ecvec #\omega_m

    wT = wTvec-wdvec #rotating frame of gate drive \omega_m-m*\omega_\gate
    transmon_levels = Qobj(diag(wT[range(N)]))
    H=transmon_levels +Omega_vec #- 0.5j*(Omega_true*adag - conj(Omega_true)*a)
    final_state = steadystate(H, c_op_list) #solve master equation

    return expect( a, final_state) #expectation value of relaxation operator

    def find_expect(vg, self=a): #phi=0.1, Omega_vec=3.0):
            phi, Omega=vg#.shape
            Omega_vec=-0.5j*(Omega*self.a_dag - conj(Omega)*self.a_op)

            Ej = self.Ejmax*absolute(cos(pi*phi)) #Josephson energy as function of Phi.

            wTvec = (-Ej + sqrt(8.0*Ej*self.Ec)*(self.nvec+0.5)+self.Ecvec)/h #\omega_m

            wT = wTvec-a.fdvec #rotating frame of gate drive \omega_m-m*\omega_\gate
            transmon_levels = Qobj(diag(wT[range(self.N_dim)]))
            H=transmon_levels +Omega_vec #- 0.5j*(Omega_true*adag - conj(Omega_true)*a)
            final_state = steadystate(H, self.c_ops) #solve master equation

            return expect( self.a_op, final_state) #expectation value of relaxation operator

def test_driven_cavity_power_bicgstab():
    "Steady state: Driven cavity - power-bicgstab solver"

    N = 30
    Omega = 0.01 * 2 * np.pi
    Gamma = 0.05

    a = destroy(N)
    H = Omega * (a.dag() + a)
    c_ops = [np.sqrt(Gamma) * a]
    M = build_preconditioner(H, c_ops, method='power')
    rho_ss = steadystate(H, c_ops, method='power-bicgstab', M=M, use_precond=1)
    rho_ss_analytic = coherent_dm(N, -1.0j * (Omega)/(Gamma/2))
    assert_((rho_ss - rho_ss_analytic).norm() < 1e-4)

def test_driven_cavity_lgmres():
    "Steady state: Driven cavity - iterative-lgmres solver"

    N = 30
    Omega = 0.01 * 2 * np.pi
    Gamma = 0.05

    a = destroy(N)
    H = Omega * (a.dag() + a)
    c_ops = [np.sqrt(Gamma) * a]

    rho_ss = steadystate(H, c_ops, method='iterative-lgmres')
    rho_ss_analytic = coherent_dm(N, -1.0j * (Omega)/(Gamma/2))

    assert_((rho_ss - rho_ss_analytic).norm() < 1e-4)

    def jc_steadystate(self, N, wc, wa, g, kappa, gamma,
                       pump, psi0, use_rwa, tlist):

        # Hamiltonian
        a = tensor(destroy(N), identity(2))
        sm = tensor(identity(N), destroy(2))

        if use_rwa:
            # use the rotating wave approxiation
            H = wc * a.dag(
            ) * a + wa * sm.dag() * sm + g * (a.dag() * sm + a * sm.dag())
        else:
            H = wc * a.dag() * a + wa * sm.dag() * sm + g * (
                a.dag() + a) * (sm + sm.dag())

        # collapse operators
        c_op_list = []

        n_th_a = 0.0  # zero temperature

        rate = kappa * (1 + n_th_a)
        c_op_list.append(np.sqrt(rate) * a)

        rate = kappa * n_th_a
        if rate > 0.0:
            c_op_list.append(np.sqrt(rate) * a.dag())

        rate = gamma
        if rate > 0.0:
            c_op_list.append(np.sqrt(rate) * sm)

        rate = pump
        if rate > 0.0:
            c_op_list.append(np.sqrt(rate) * sm.dag())

        # find the steady state
        rho_ss = steadystate(H, c_op_list)

        return expect(a.dag() * a, rho_ss), expect(sm.dag() * sm, rho_ss)

def _spectrum_pi(H, wlist, c_ops, a_op, b_op, use_pinv=False):
    """
    Internal function for calculating the spectrum of the correlation function
    :math:`\left<A(\\tau)B(0)\\right>`.
    """

    #print issuper(H)

    L = H if issuper(H) else liouvillian(H, c_ops)

    tr_mat = tensor([qeye(n) for n in L.dims[0][0]])
    N = prod(L.dims[0][0])

    A = L.full()
    b = spre(b_op).full()
    a = spre(a_op).full()

    tr_vec = transpose(mat2vec(tr_mat.full()))

    rho_ss = steadystate(L)
    rho = transpose(mat2vec(rho_ss.full()))

    I = identity(N * N)
    P = kron(transpose(rho), tr_vec)
    Q = I - P

    spectrum = zeros(len(wlist))

    for idx, w in enumerate(wlist):
        if use_pinv:
            MMR = pinv(-1.0j * w * I + A)
        else:
            MMR = dot(Q, solve(-1.0j * w * I + A, Q))

        s = dot(tr_vec,
                   dot(a, dot(MMR, dot(b, transpose(rho)))))
        spectrum[idx] = -2 * real(s[0, 0])

    return spectrum

def find_expect(phi=0.1, Omega_vec=3.0, wd=wd):

    wdvec=nvec*wd

    Ej = Ejmax*absolute(cos(pi*phi)) #; % Josephson energy as function of Phi.
    wq0=sqrt(8*Ej*Ec)
    X=Np*pi*(wq0-f0)/f0
    Ba=-1*gamma*(sin(2.0*X)-2.0*X)/(2.0*X**2.0)
    Ecp=Ec*(1-2*Ba/wq0)
    Ecvec=-Ecp*(6.0*nvec**2+6.0*nvec+3.0)/12.0

    wTvec = -Ej + sqrt(8.0*Ej*Ecp)*(nvec+0.5)+Ecvec


    #print diff(wTvec)
    if phi==0.4:
        print Ba/wq0
    if 0:
        zr=zeros(N)
        X=Np*pi*(diff(wTvec)-f0)/f0
        Ba=-0*gamma*(sin(2.0*X)-2.0*X)/(2.0*X**2.0)
        ls=-Ba#/(2*Ct)/(2*pi)#/1e6
        zr[1:]=ls
        #print Ga0/(2*Ct)/(2*pi)
        #print zr
        wTvec=wTvec+zr
    wTvec = wTvec-wdvec

    #print wTvec+diff(wTvec)


    transmon_levels = Qobj(diag(wTvec[range(N)]))

    H=transmon_levels +Omega_vec #- 0.5j*(Omega_true*adag - conj(Omega_true)*a)
    final_state = steadystate(H, c_op_list)

    return expect( tm_l, final_state)

wVals = np.linspace(-2 * np.pi, 6 * np.pi, 100)

# Solving the dynamics of a spin 1/2 in B0 and Brf in the lab frame

def H(w):
    """ Hamiltonian in the lab frame. """
    return w0 / 2 * qt.sigmaz() + w1 * np.cos(w) * qt.sigmax()

c1 = np.sqrt(1.0/T1) * qt.destroy(2)
c2 = np.sqrt(1.0/T2) / 2 * qt.destroy(2)
c_op_list = [c1,c2]


wHVals = np.linspace(-2 * np.pi, 6 * np.pi, 51)
MxH = [qt.expect(qt.sigmax(), qt.steadystate(H(w), c_op_list, maxiter=100)) for w in wHVals]
MyH = [qt.expect(qt.sigmay(), qt.steadystate(H(w), c_op_list, maxiter=100)) for w in wHVals]
MzH = [qt.expect(qt.sigmaz(), qt.steadystate(H(w), c_op_list, maxiter=100)) for w in wHVals]


# Plot results of the algebraic formula and the solution given by qutip

plt.ylim([-1.0, 1.0])

plt.plot(wVals, Mx(wVals), label='Mx')
plt.plot(wVals, My(wVals), label='My')
plt.plot(wVals, Mz(wVals), label='Mz')

plt.plot(wHVals, MxH, 'bo', label='MxH')
plt.plot(wHVals, MyH, 'go', label='MyH')
plt.plot(wHVals, MzH, 'ro', label='MzH')

def find_expect(phi=0.1, Omega_vec=3.0, wd=wd, use_pinv=True):

    Ej = (phi**2)/(8*Ec) #Ejmax*absolute(cos(pi*phi)) #; % Josephson energy as function of Phi.

    wTvec = -Ej + sqrt(8.0*Ej*Ec)*(nvec+0.5)+Ecvec

    #wdvec=nvec*sqrt(8.0*Ej*Ec)
    #wdvec=nvec*wd

    wT = wTvec-wdvec
    transmon_levels = Qobj(diag(wT[range(N)]))
    H=transmon_levels +Omega_vec #- 0.5j*(Omega_true*adag - conj(Omega_true)*a)
    L=liouvillian(H, c_op_list) #ends at same thing as L = H_comm + Lindblad_tm ;
    rho_ss = steadystate(L) #same as rho_ss_c but that's in column vector form

    tr_mat = tensor([qeye(n) for n in L.dims[0][0]])
    #N = prod(L.dims[0][0])

    A = L.full()
    #D, U= L.eigenstates()
    #print D.shape, D
    #print diag(D)
    b = spre(p).full()-spost(p).full()
    a2 = spre(a).full()

    tr_vec = transpose(mat2vec(tr_mat.full()))

    rho = transpose(mat2vec(rho_ss.full()))

    I = identity(N * N)
    P = kron(transpose(rho), tr_vec)
    Q = I - P

    #wp=4.9
    #w=-(wp-wd)
    spectrum = zeros(len(wlist2), dtype=complex)
    for idx, w in enumerate(wlist2):

        if use_pinv:
            MMR = pinv(-1.0j * w * I + A) #eig(MMR)[0] is equiv to Dint
        else:
            MMR = dot(Q, solve(-1.0j * w * I + A, Q))
    #print diag(1.0/(1j*(wp-wd)*ones(N**2)+D)) #Dint = diag(1./(1i*(wp-wd)*diag(eye(dim^2)) + diag(D)))

    #print 1.0/(1j*(wp-wd)*ones(N**2)+D) #Dint = diag(1./(1i*(wp-wd)*diag(eye(dim^2)) + diag(D)))
    #U2=squeeze([u.full() for u in U]).transpose()
    #Dint=eig(MMR)[0]
    #print "MMR", eig(MMR)[1]
    #print "Umult", U2*Dint*inv(U2)

        s = dot(tr_vec,
            dot(a2, dot(MMR, dot(b, transpose(rho)))))
        #spectrum[idx] = -2 * real(s[0, 0])
        spectrum[idx]=1j*s[0][0] #matches Chi_temp result #(1/theta_steps)
    return spectrum
    #final_state = steadystate(H, c_op_list)
    #print H.shape
    #print dir(H)
    #U,D = eig(H.full())
    #print D
    #Uinv = Qobj(inv(D))
    #U=Qobj(D)
    # Doing the chi integral (gives the susceptibility)

    #Dint = 1.0/(1.0j*(wp-wd)) #Qobj(1.0/(1.0j*(wp-wd)*diag(qeye(N**2))))# + diag(D))))

    #Hint = H.expm() #U*H*Uinv

    #Chi_temp(i,j) += (1.0/theta_steps)*1j*trace(reshape(tm_l*Lint*(p_l - p_r)*rho_ss_c,dim,dim))
    #exp1=correlation_2op_1t(H, None, wlist2, c_op_list, a, p, solver="es")
    #exp2=correlation_2op_1t(H, None, wlist2, c_op_list, p, a, solver="es", reverse=True)
    #exp1=correlation_2op_1t(H, None, wlist2, c_op_list, a, p_l-p_r, solver="es")

    #exp1=correlation_ss(H, wlist2, c_op_list, a, p)
    #exp2=correlation_ss(H, wlist2, c_op_list, p, a, reverse=True)
    exp1=spectrum(H, wlist2, c_op_list, a, p, solver="pi", use_pinv=False)
    exp2=spectrum(H, wlist2, c_op_list, p, a, solver="pi", use_pinv=False)
    return exp1-exp2
    return expect( a, final_state) #tm_l

print Omega_vec
H=transmon_levels +Omega_vec #- 0.5j*(Omega_true*adag - conj(Omega_true)*a)
print H
#Ht = H.data.T
#from qutip.cy.spmath import zcsr_kron

#from qutip.fastsparse import fast_csr_matrix, fast_identity
#spI=fast_identity(N)
#data = -1j * zcsr_kron(spI, H.data)
#data += 1j * zcsr_kron(Ht, spI)
#print data.toarray()
L=liouvillian(H, c_op_list) #ends at same thing as L = H_comm + Lindblad_tm ;
#print L
#print dir(L)

rho_ss = steadystate(L) #same as rho_ss_c but that's in column vector form
print rho_ss

wlist = linspace(-500.0, 500.0, 201)

use_pinv=True

tr_mat = tensor([qeye(n) for n in L.dims[0][0]])
N = prod(L.dims[0][0])

A = L.full()
D, U= L.eigenstates()
print D.shape, D
print diag(D)
b = spre(p).full()-spost(p).full()
a = spre(a).full()

# Calculating steady-state solutions using the qutip library

def H(w):
    """ Hamiltonian of a spin 1/2 in B0 + Brf using the RWA. """
    return (w0 - w)/2 * qt.sigmaz() + w1/2 * qt.sigmax()

# Collsapse operators which should describe the relaxation of the system
c1 = np.sqrt(1.0/T1) * qt.destroy(2)
c2 = np.sqrt(1.0/T2)/2 * qt.destroy(2)
c3 = np.sqrt(1.0/(T1*T2)) * qt.destroy(2)
c4 = np.sqrt(1.0/T2) /2 * qt.create(2)
c5 = np.sqrt(1.0/(T1*T2)) * qt.create(2)


wHVals = np.linspace(-2 * np.pi, 6 * np.pi, 31)
MxH = [qt.expect(qt.sigmax(), qt.steadystate(H(w), [c1,c2,c3], maxiter=10)) for w in wHVals]
MyH = [qt.expect(qt.sigmay(), qt.steadystate(H(w), [c2,c5], maxiter=10)) for w in wHVals]
MzH = [qt.expect(qt.sigmaz(), qt.steadystate(H(w), [c1], maxiter=10)) for w in wHVals]

wVals = np.linspace(-2 * np.pi, 6 * np.pi, 100)

# Plot results of the algebraic formula and the solution given by qutip

plt.ylim([-1.0, 1.0])

plt.plot(wVals, Mx(wVals), label='Mx')
plt.plot(wVals, My(wVals), label='My')
plt.plot(wVals, Mz(wVals), label='Mz')

plt.plot(wHVals, MxH, 'bo', label='MxH')
plt.plot(wHVals, MyH, 'go', label='MyH')

def two_tone(vg):
    phi, wd=vg
    Ej = Ejmax*absolute(cos(pi*phi))#(phi**2)/(8*Ec) # #; % Josephson energy as function of Phi.

    wTvec = -Ej + sqrt(8.0*Ej*Ec)*(nvec+0.5)+Ecvec
    wdvec=nvec*wd
    wT = wTvec-wdvec

    transmon_levels = Qobj(diag(wT[range(N)]))
    

    for Om in Omega_op:
        H = transmon_levels + Om
        exp1=spectrum(H, wlist2, c_op_list, a, p, solver="pi", use_pinv=False)
        exp2=spectrum(H, wlist2, c_op_list, p, a, solver="pi", use_pinv=False)
        return exp1-exp2


        #H_comm = -1j*(kron(It,H) - kron(transpose(H),It));

        #L = H_comm + Lindblad_tm + Lindblad_tp #+ Lindblad_deph;

        #print L

        #L2 = [reshape(eye(N),1,N**2); L]; %Add the condition trace(rho) = 1

        #rho_ss = L2\[1;zeros(dim^2,1)]; % Steady state density matrix

        #rho_ss_c = rho_ss; %Column vector form
        L = liouvillian(H, c_op_list)

        #D, V=L.eigenstates()
        #print "D", D.shape
        #print "V", V.shape
        #print L.diag()
        #print L.diag().shape
        #raise Exception
        A = L.full()

        rho_ss = steadystate(L)
        rho = transpose(mat2vec(rho_ss.full()))
        P = kron(transpose(rho), tr_vec)
        Q = I - P

        if 1:
            w=wp-wd
            #MMR = pinv(-1.0j * w * I + A)
            MMR = dot(Q, solve(-1.0j * w * I + A, Q))

            #Lexp=L.expm()
            #MMR=Lexp*MMR*Lexp.dag()
            #    MMR = np.dot(Q, np.linalg.solve(-1.0j * w * I + A, Q))

            #print p_l.shape
            #print p_l
            #print transpose(rho).shape

            #print dot(p_l.full(), transpose(rho))
            s = dot(tr_vec,
                       dot(tm_l, dot(MMR, dot(p_l_m_p_r, transpose(rho)))))
            return -2 * real(s[0, 0])
        D,U = eig(L.full())

        Uinv = inv(U)

        Dint = diag(1.0/(1.0j*(wp-wd)*Idiag + diag(D)))

        Lint = U*Dint*Uinv


        #H_comm = -1.0j*(kron(It,H) - kron(transpose(H),It))


        #print H
        #print H.shape, H.full().shape
        #print H_comm, H_comm.shape
        #L = H_comm + Lindblad_tm + Lindblad_tp + Lindblad_deph
        #p_l=spre(p)
        #p_r=spost(p)
        #tm_l=spre(a)
        #print L.shape, A.shape,
        #print dir(L)
        #print L.eigenenergies().shape
        #print L.eigenstates()#.shape
        #L2 = [reshape(eye(dim),1,dim^2);L]; %Add the condition trace(rho) = 1
        #rho_ss = L2\[1;zeros(dim^2,1)]; % Steady state density matrix
        #rho_ss_c = rho_ss; %Column vector form

        #print "tm_l", tm_l.full().shape
        #print "Lint", Lint.shape
        #print "pl", p_l.full().shape
        #print "pr", p_r.full().shape
        #print rho.shape, to_super(rho_ss).full().shape
        #print "rho", to_super(rho_ss) #rho_ss.full().shape #rho.shape
        #print (tm_l.full()*Lint*(p_l.full() - p_r.full())*rho).shape
        #raise Exception
        return 1j*trace(tm_l*Lint*(p_l_m_p_r)*rho_ss)