Merge 68b915b into 014c86b

src/QuantumOptics.jl

@@ -20,6 +20,7 @@ export bases, Basis, GenericBasis, CompositeBasis,
         subspace, SubspaceBasis, projector,
         particle, PositionBasis, MomentumBasis, gaussianstate,
                 positionoperator, momentumoperator, laplace_x, laplace_p,
+        nlevel, NLevelBasis, transition, nlevelstate,
         nparticlebasis, BosonicNParticleBasis, FermionicNParticleBasis, nparticleoperator_1, nparticleoperator_2,
         metrics, tracedistance,
         spectralanalysis, operatorspectrum, operatorspectrum_hermitian,
@@ -42,6 +43,7 @@ include("spin.jl")
 include("fock.jl")
 include("subspace.jl")
 include("particle.jl")
+include("nlevel.jl")
 include("nparticles.jl")
 include("metrics.jl")
 include("ode_dopri.jl")
@@ -71,6 +73,7 @@ using .spin
 using .fock
 using .subspace
 using .particle
+using .nlevel
 using .nparticles
 using .metrics
 using .spectralanalysis

src/nlevel.jl

@@ -0,0 +1,34 @@
+module nlevel
+
+import Base.==
+
+using ..bases, ..states, ..operators, ..operators_sparse
+
+export NLevelBasis, transition, nlevelstate
+
+
+type NLevelBasis <: Basis
+    shape::Vector{Int}
+    N::Int
+    function NLevelBasis(N::Int)
+        if N < 1
+            throw(DimensionMismatch())
+        end
+        new([N], N)
+    end
+end
+
+==(b1::NLevelBasis, b2::NLevelBasis) = b1.N == b2.N
+
+function transition(b::NLevelBasis, to::Int, from::Int)
+    op = SparseOperator(b)
+    op.data[to, from] = 1.
+    op
+end
+
+function nlevelstate(b::NLevelBasis, n::Int)
+    @assert n <= b.N
+    basis_ket(b, n)
+end
+
+end # module

src/steadystate.jl

@@ -98,7 +98,15 @@ function eigenvector(L::DenseSuperOperator)
 end
 
 function eigenvector(L::SparseSuperOperator)
-    d, v, nconv, niter, nmult, resid = Base.eigs(L.data; nev=1, sigma=1e-30)
+    d, v, nconv, niter, nmult, resid = try
+      Base.eigs(L.data; nev=1, sigma=1e-30)
+    catch err
+      if isa(err, LinAlg.SingularException)
+        error("Base.LinAlg.eigs() algorithm failed; try using DenseOperators")
+      else
+        rethrow(err)
+      end
+    end
     data = reshape(v[:,1], length(L.basis_r[1]), length(L.basis_r[2]))
     op = DenseOperator(L.basis_r[1], L.basis_r[2], data)
     return op/trace(op)

test/test_nlevel.jl

@@ -0,0 +1,11 @@
+using Base.Test
+using QuantumOptics
+
+energies = [0., 1.1, 2.6]
+N = 3
+b = NLevelBasis(N)
+
+@test b == NLevelBasis(N)
+@test b != NLevelBasis(N+1)
+@test full(transition(b, 2, 1)) == basis_ket(b, 2) ⊗ basis_bra(b, 1)
+@test nlevelstate(b, 2) == Ket(b, [0, 1, 0])

test/test_steadystate.jl

@@ -66,3 +66,7 @@ nss = expect(create(fockbasis)*destroy(fockbasis), ρss)
 ρss = steadystate.eigenvector(full(Hp), map(full, Jp))
 nss = expect(create(fockbasis)*destroy(fockbasis), ρss)
 @test n_an - real(nss) < 1e-3
+
+
+# Test error messages
+@test_throws ErrorException steadystate.eigenvector(sx, [sm])