Make the plot generation more robust

.gitignore

@@ -3,4 +3,5 @@ src/documentation
 src/_benchmarks-sourcecode
 src/_codesnippets/build
 src/images/codesnippets
-src/benchmark-data
+src/benchmark-data
+Manifest.toml

Project.toml

@@ -0,0 +1,3 @@
+[deps]
+PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
+QuantumOptics = "6e0679c1-51ea-5a7c-ac74-d61b76210b0c"

make.jl

@@ -1,3 +1,9 @@
+
+if normpath(Base.active_project()) != normpath(joinpath(@__DIR__, "Project.toml"))
+    #If the current directory isn't the active Julia project, make it so
+    import Pkg
+    Pkg.activate(@__DIR__)
+end
+
 # Build Code Snippets
-cd("src/_codesnippets")
-run(`julia make.jl`)
+include("src/_codesnippets/make.jl")

src/_codesnippets/make.jl

@@ -2,13 +2,15 @@ sourcedir = "src"
 builddir = "build"
 imagedir = "../../images/codesnippets"
 
+cd(@__DIR__)
 cp(sourcedir, builddir; force=true)
 cd(builddir)
 
-names = filter(name->endswith(name, ".jl"), readdir("."))
+names = filter(endswith(".jl"), readdir("."))
 for name in names
     println("Executing $name")
-    run(`julia $name`)
+    include(joinpath(builddir, "$name"))
+    PyPlot.clf()
 end
 
 imagenames = filter(name->endswith(name, ".svg")||endswith(name, ".png"), readdir("."))

src/_codesnippets/src/composite.jl

@@ -1,18 +1,19 @@
 using QuantumOptics
-b_spin = SpinBasis(1//2)
-b_fock = FockBasis(200)
-sp = sigmap(b_spin)
-sm = sigmam(b_spin)
-a = destroy(b_fock)
-at = create(b_fock)
-H = sp⊗a + sm⊗at
-T = [0:0.01:50;]
-ψ0 = spindown(b_spin) ⊗ coherentstate(b_fock, 6)
-tout, ψt = timeevolution.schroedinger(T, ψ0, H)
 
-using PyPlot
-plot(tout, expect(1, sp*sm, ψt))
-xlabel("Time")
-ylabel("Spin excitation")
-tight_layout()
-savefig("composite.svg")
+let b_spin = SpinBasis(1 // 2), b_fock = FockBasis(200)
+    sp = sigmap(b_spin)
+    sm = sigmam(b_spin)
+    a = destroy(b_fock)
+    at = create(b_fock)
+    H = sp ⊗ a + sm ⊗ at
+    T = [0:0.01:50;]
+    ψ0 = spindown(b_spin) ⊗ coherentstate(b_fock, 6)
+    tout, ψt = timeevolution.schroedinger(T, ψ0, H)
+
+    using PyPlot
+    plot(tout, expect(1, sp * sm, ψt))
+    xlabel("Time")
+    ylabel("Spin excitation")
+    tight_layout()
+    savefig("composite.svg")
+end

src/_codesnippets/src/cooling.jl

@@ -1,25 +1,26 @@
 using QuantumOptics
-bc = FockBasis(16)
-ba = SpinBasis(1//2)
-a = destroy(bc) ⊗ one(ba)
-sm = one(bc) ⊗ sigmam(ba)
-H0 = 0.5*sm'*sm + a + a'
-Hx = 0.5*(a'*sm + a*sm')
-J = [a, sqrt(2)*sm]
-Jdagger = dagger.(J)
 
-fquantum(t, ψ, u) = H0 + Hx*cos(u[1]), J, Jdagger
-function fclassical(du, u, ψ, t)
-  du[1] = 0.3*real(u[2])
-  du[2] = sin(real(u[1]))*real(expect(dagger(a)*sm, ψ))
-end
-u0 = ComplexF64[sqrt(2), 6.0]
-ψ0 = semiclassical.State(fockstate(bc, 0)⊗spindown(ba), u0)
-tout, p = semiclassical.master_dynamic([0:1.0:200;], ψ0, fquantum,
-      fclassical; fout=(t,rho)->rho.classical[2])
+let bc = FockBasis(16), ba = SpinBasis(1 // 2)
+  a = destroy(bc) ⊗ one(ba)
+  sm = one(bc) ⊗ sigmam(ba)
+  H0 = 0.5 * sm' * sm + a + a'
+  Hx = 0.5 * (a' * sm + a * sm')
+  J = [a, sqrt(2) * sm]
+  Jdagger = dagger.(J)
 
-using PyPlot
-plot(tout, p.^2)
-ylabel("Kinetic energy")
-xlabel("Time")
-savefig("cooling.svg")
+  fquantum(t, ψ, u) = H0 + Hx * cos(u[1]), J, Jdagger
+  function fclassical(du, u, ψ, t)
+    du[1] = 0.3 * real(u[2])
+    du[2] = sin(real(u[1])) * real(expect(dagger(a) * sm, ψ))
+  end
+  u0 = ComplexF64[sqrt(2), 6.0]
+  ψ0 = semiclassical.State(fockstate(bc, 0) ⊗ spindown(ba), u0)
+  tout, p = semiclassical.master_dynamic([0:1.0:200;], ψ0, fquantum,
+    fclassical; fout=(t, rho) -> rho.classical[2])
+
+  using PyPlot
+  plot(tout, p .^ 2)
+  ylabel("Kinetic energy")
+  xlabel("Time")
+  savefig("cooling.svg")
+end

src/_codesnippets/src/fock.jl

@@ -1,17 +1,19 @@
 using QuantumOptics
-b = FockBasis(50)
-a = destroy(b)
-at = create(b)
-H = 0.5*(a^2 + at^2)
-psi0 = fockstate(b, 3)
-tout, psit = timeevolution.schroedinger([0:0.25:1;], psi0, H)
 
-using PyPlot
-x = [-5:0.1:5;]
-for i in 1:4
-    subplot(2, 2, i)
-    Q = qfunc(psit[i], x, x)
-    pcolor(x, x, Q)
-end
-tight_layout()
-savefig("fock.png")
+let b = FockBasis(50)
+    a = destroy(b)
+    at = create(b)
+    H = 0.5 * (a^2 + at^2)
+    psi0 = fockstate(b, 3)
+    tout, psit = timeevolution.schroedinger([0:0.25:1;], psi0, H)
+
+    using PyPlot
+    x = [-5:0.1:5;]
+    for i in 1:4
+        subplot(2, 2, i)
+        Q = qfunc(psit[i], x, x)
+        pcolor(x, x, Q)
+    end
+    tight_layout()
+    savefig("fock.png")
+end

src/_codesnippets/src/gross-pitaevskii.jl

@@ -1,23 +1,24 @@
 using QuantumOptics
-b_x = PositionBasis(-10, 10, 300)
-b_p = MomentumBasis(b_x)
-Tpx = transform(b_p, b_x)
-Txp = transform(b_x, b_p)
-x = position(b_x)
-p = momentum(b_p)
 
-Hkin = LazyProduct(Txp, p^2/2, Tpx)
-Hpsi = diagonaloperator(b_x, Ket(b_x).data)
-H = LazySum(Hkin, Hpsi)
+let b_x = PositionBasis(-10, 10, 300), b_p = MomentumBasis(b_x)
+    Tpx = transform(b_p, b_x)
+    Txp = transform(b_x, b_p)
+    x = position(b_x)
+    p = momentum(b_p)
 
-fquantum(t, ψ) = (Hpsi.data.nzval .= -50 .* abs2.(ψ.data); H)
-ψ₀ = gaussianstate(b_x,-2π,2,1.5) + gaussianstate(b_x,2π,-2,1.5)
-normalize!(ψ₀)
-tout, ψₜ = timeevolution.schroedinger_dynamic([0:0.01:6;], ψ₀, fquantum)
-density = [abs2(ψ.data[j]) for ψ=ψₜ, j=1:length(b_x)]
+    Hkin = LazyProduct(Txp, p^2 / 2, Tpx)
+    Hpsi = diagonaloperator(b_x, Ket(b_x).data)
+    H = LazySum(Hkin, Hpsi)
 
-using PyPlot
-contourf(samplepoints(b_x),tout,density,50)
-xlabel("x")
-ylabel("Time")
-savefig("gross_pitaevskii.png")
+    fquantum(t, ψ) = (Hpsi.data.nzval .= -50 .* abs2.(ψ.data); H)
+    ψ₀ = gaussianstate(b_x, -2π, 2, 1.5) + gaussianstate(b_x, 2π, -2, 1.5)
+    normalize!(ψ₀)
+    tout, ψₜ = timeevolution.schroedinger_dynamic([0:0.01:6;], ψ₀, fquantum)
+    density = [abs2(ψ.data[j]) for ψ = ψₜ, j = 1:length(b_x)]
+
+    using PyPlot
+    contourf(samplepoints(b_x), tout, density, 50)
+    xlabel("x")
+    ylabel("Time")
+    savefig("gross_pitaevskii.png")
+end

src/_codesnippets/src/nlevel.jl

@@ -1,20 +1,22 @@
 using QuantumOptics
-b = NLevelBasis(3)
-t12 = transition(b, 1, 2)
-t23 = transition(b, 2, 3)
-t31 = transition(b, 1, 3)
-H = 10*(t31 + dagger(t31))
-J = [1.2*t23, 0.6*t12]
-psi0 = basisstate(b, 1)
-T = [0:0.01:10;]
-tout, psit = timeevolution.mcwf(T, psi0, H, J; seed=2)
 
-using PyPlot
-plot(tout, expect(dm(basisstate(b, 3)), psit), label=L"$|3\rangle$")
-plot(tout, expect(dm(basisstate(b, 2)), psit), label=L"$|2\rangle$")
-plot(tout, expect(dm(basisstate(b, 1)), psit), label=L"$|1\rangle$")
-xlabel("Time")
-ylabel("Probability")
-legend()
-tight_layout()
-savefig("nlevel.svg")
+let b = NLevelBasis(3)
+    t12 = transition(b, 1, 2)
+    t23 = transition(b, 2, 3)
+    t31 = transition(b, 1, 3)
+    H = 10 * (t31 + dagger(t31))
+    J = [1.2 * t23, 0.6 * t12]
+    psi0 = basisstate(b, 1)
+    T = [0:0.01:10;]
+    tout, psit = timeevolution.mcwf(T, psi0, H, J; seed=2)
+
+    using PyPlot
+    plot(tout, expect(dm(basisstate(b, 3)), psit), label=L"$|3\rangle$")
+    plot(tout, expect(dm(basisstate(b, 2)), psit), label=L"$|2\rangle$")
+    plot(tout, expect(dm(basisstate(b, 1)), psit), label=L"$|1\rangle$")
+    xlabel("Time")
+    ylabel("Probability")
+    legend()
+    tight_layout()
+    savefig("nlevel.svg")
+end

src/_codesnippets/src/particle.jl

@@ -1,18 +1,20 @@
 using QuantumOptics
-basis = PositionBasis(-2, 2, 200)
-x = position(basis)
-p = momentum(basis)
-H = p^2/4 + 2*DenseOperator(x^2)
-energies, states = eigenstates((H+dagger(H))/2, 5)
 
-using PyPlot
-xpoints = samplepoints(basis)
-plot(xpoints, 2*xpoints.^2)
-fill_between(xpoints, 0., 2*xpoints.^2, alpha=0.5)
-for i=1:length(states)
-    plot(xpoints, abs2.(states[i].data).*40 .+ energies[i])
-end
-xlabel("Position")
-ylabel("Energy")
-tight_layout()
-savefig("particle.svg")
+let basis = PositionBasis(-2, 2, 200)
+    x = position(basis)
+    p = momentum(basis)
+    H = p^2 / 4 + 2 * DenseOperator(x^2)
+    energies, states = eigenstates((H + dagger(H)) / 2, 5)
+
+    using PyPlot
+    xpoints = samplepoints(basis)
+    plot(xpoints, 2 * xpoints .^ 2)
+    fill_between(xpoints, 0.0, 2 * xpoints .^ 2, alpha=0.5)
+    for i in eachindex(states)
+        plot(xpoints, abs2.(states[i].data) .* 40 .+ energies[i])
+    end
+    xlabel("Position")
+    ylabel("Energy")
+    tight_layout()
+    savefig("particle.svg")
+end

src/_codesnippets/src/spin.jl

@@ -1,22 +1,24 @@
 using QuantumOptics
-b = SpinBasis(3//2)
-sm = sigmam(b)
-H = 2*sigmaz(b)
-J = [sm]
-τ = [0:0.025:5;]
-ω = [-5:0.05:25;]
-ρ0 = dm(spinup(b))
-corr = timecorrelations.correlation(τ, ρ0, H, J, sigmap(b), sm)
-ω, S = timecorrelations.spectrum(ω, H, J, sm)
 
-using PyPlot
-subplot(2, 1, 1)
-plot(τ, corr)
-xlabel(L"\tau")
-ylabel(L"\langle \sigma_+(\tau) \sigma_-(0)\rangle")
-subplot(2, 1, 2)
-plot(ω, S)
-xlabel(L"\omega")
-ylabel(L"S(\omega)")
-tight_layout()
-savefig("spin.svg")
+let b = SpinBasis(3 // 2)
+    sm = sigmam(b)
+    H = 2 * sigmaz(b)
+    J = [sm]
+    τ = [0:0.025:5;]
+    ω = [-5:0.05:25;]
+    ρ0 = dm(spinup(b))
+    corr = timecorrelations.correlation(τ, ρ0, H, J, sigmap(b), sm)
+    ω, S = timecorrelations.spectrum(ω, H, J, sm)
+
+    using PyPlot
+    subplot(2, 1, 1)
+    plot(τ, corr)
+    xlabel(L"\tau")
+    ylabel(L"\langle \sigma_+(\tau) \sigma_-(0)\rangle")
+    subplot(2, 1, 2)
+    plot(ω, S)
+    xlabel(L"\omega")
+    ylabel(L"S(\omega)")
+    tight_layout()
+    savefig("spin.svg")
+end