def lindblad_terms(gamma, rho, Ne, verbose=1):
    u"""Return the Lindblad terms for decays gamma in matrix form.

    >>> from sympy import pprint
    >>> aux = define_frequencies(4, explicitly_antisymmetric=True)
    >>> omega_level, omega, gamma = aux
    >>> gamma = gamma.subs({gamma[2, 0]:0, gamma[3, 0]:0, gamma[3, 1]:0})
    >>> pprint(gamma, use_unicode=True)
    ⎡ 0   -γ₂₁   0     0  ⎤
    ⎢                     ⎥
    ⎢γ₂₁   0    -γ₃₂   0  ⎥
    ⎢                     ⎥
    ⎢ 0   γ₃₂    0    -γ₄₃⎥
    ⎢                     ⎥
    ⎣ 0    0    γ₄₃    0  ⎦
    >>> rho = define_density_matrix(4)
    >>> pprint(lindblad_terms(gamma, rho, 4), use_unicode=True)
    ⎡                -γ₂₁⋅ρ₁₂             -γ₃₂⋅ρ₁₃             -γ₄₃⋅ρ₁₄      ⎤
    ⎢ γ₂₁⋅ρ₂₂        ─────────            ─────────            ─────────     ⎥
    ⎢                    2                    2                    2         ⎥
    ⎢                                                                        ⎥
    ⎢-γ₂₁⋅ρ₂₁                          γ₂₁⋅ρ₂₃   γ₃₂⋅ρ₂₃    γ₂₁⋅ρ₂₄   γ₄₃⋅ρ₂₄⎥
    ⎢─────────  -γ₂₁⋅ρ₂₂ + γ₃₂⋅ρ₃₃   - ─────── - ───────  - ─────── - ───────⎥
    ⎢    2                                2         2          2         2   ⎥
    ⎢                                                                        ⎥
    ⎢-γ₃₂⋅ρ₃₁     γ₂₁⋅ρ₃₂   γ₃₂⋅ρ₃₂                         γ₃₂⋅ρ₃₄   γ₄₃⋅ρ₃₄⎥
    ⎢─────────  - ─────── - ───────  -γ₃₂⋅ρ₃₃ + γ₄₃⋅ρ₄₄   - ─────── - ───────⎥
    ⎢    2           2         2                               2         2   ⎥
    ⎢                                                                        ⎥
    ⎢-γ₄₃⋅ρ₄₁     γ₂₁⋅ρ₄₂   γ₄₃⋅ρ₄₂    γ₃₂⋅ρ₄₃   γ₄₃⋅ρ₄₃                     ⎥
    ⎢─────────  - ─────── - ───────  - ─────── - ───────       -γ₄₃⋅ρ₄₄      ⎥
    ⎣    2           2         2          2         2                        ⎦

    Notice that there are more terms than simply adding a decay
    gamma_ij*rho_ij/2 for each coherence.

    """
    # We count the necessary Lindblad operators.
    Nterms = 0
    for i in range(Ne):
        for j in range(i):
            if gamma[i, j] != 0:
                Nterms += 1

    L = zeros(Ne)
    counter = 0
    t0 = time()
    for i in range(Ne):
        for j in range(i):
            if gamma[i, j] != 0:
                counter += 1
                sig = ket(j+1, Ne)*bra(i+1, Ne)
                L += gamma[i, j]*lindblad_operator(sig, rho)
                tn = time()
                if tn-t0 > 1:
                    aux = "Calculated up to i={}, j={}, or {}/{} = {:2.2f} %."
                    if verbose > 0:
                        print(aux.format(i, j, counter, Nterms,
                                         float(counter+1)/Nterms*100))
                    t0 = tn
    return L
def lindblad_terms ( gamma , rho , Ne , verbose = 1 ) : Nterms = 0 for i in range ( Ne ) : for j in range ( i ) : if gamma [ i , j ] != 0 : Nterms += 1 L = zeros ( Ne ) counter = 0 t0 = time ( ) for i in range ( Ne ) : for j in range ( i ) : if gamma [ i , j ] != 0 : counter += 1 sig = ket ( j + 1 , Ne ) * bra ( i + 1 , Ne ) L += gamma [ i , j ] * lindblad_operator ( sig , rho ) tn = time ( ) if tn - t0 > 1 : aux = "Calculated up to i={}, j={}, or {}/{} = {:2.2f} %." if verbose > 0 : print ( aux . format ( i , j , counter , Nterms , float ( counter + 1 ) / Nterms * 100 ) ) t0 = tn return L
u"""Return the Lindblad terms for decays gamma in matrix form.

    >>> from sympy import pprint
    >>> aux = define_frequencies(4, explicitly_antisymmetric=True)
    >>> omega_level, omega, gamma = aux
    >>> gamma = gamma.subs({gamma[2, 0]:0, gamma[3, 0]:0, gamma[3, 1]:0})
    >>> pprint(gamma, use_unicode=True)
    ⎡ 0   -γ₂₁   0     0  ⎤
    ⎢                     ⎥
    ⎢γ₂₁   0    -γ₃₂   0  ⎥
    ⎢                     ⎥
    ⎢ 0   γ₃₂    0    -γ₄₃⎥
    ⎢                     ⎥
    ⎣ 0    0    γ₄₃    0  ⎦
    >>> rho = define_density_matrix(4)
    >>> pprint(lindblad_terms(gamma, rho, 4), use_unicode=True)
    ⎡                -γ₂₁⋅ρ₁₂             -γ₃₂⋅ρ₁₃             -γ₄₃⋅ρ₁₄      ⎤
    ⎢ γ₂₁⋅ρ₂₂        ─────────            ─────────            ─────────     ⎥
    ⎢                    2                    2                    2         ⎥
    ⎢                                                                        ⎥
    ⎢-γ₂₁⋅ρ₂₁                          γ₂₁⋅ρ₂₃   γ₃₂⋅ρ₂₃    γ₂₁⋅ρ₂₄   γ₄₃⋅ρ₂₄⎥
    ⎢─────────  -γ₂₁⋅ρ₂₂ + γ₃₂⋅ρ₃₃   - ─────── - ───────  - ─────── - ───────⎥
    ⎢    2                                2         2          2         2   ⎥
    ⎢                                                                        ⎥
    ⎢-γ₃₂⋅ρ₃₁     γ₂₁⋅ρ₃₂   γ₃₂⋅ρ₃₂                         γ₃₂⋅ρ₃₄   γ₄₃⋅ρ₃₄⎥
    ⎢─────────  - ─────── - ───────  -γ₃₂⋅ρ₃₃ + γ₄₃⋅ρ₄₄   - ─────── - ───────⎥
    ⎢    2           2         2                               2         2   ⎥
    ⎢                                                                        ⎥
    ⎢-γ₄₃⋅ρ₄₁     γ₂₁⋅ρ₄₂   γ₄₃⋅ρ₄₂    γ₃₂⋅ρ₄₃   γ₄₃⋅ρ₄₃                     ⎥
    ⎢─────────  - ─────── - ───────  - ─────── - ───────       -γ₄₃⋅ρ₄₄      ⎥
    ⎣    2           2         2          2         2                        ⎦

    Notice that there are more terms than simply adding a decay
    gamma_ij*rho_ij/2 for each coherence.