using JaynesCummings, PlotlyJS N = 120 ω_q = 1.0 ω_r = 1.0 w_samples = 100 gnumber = 251 couplings = linspace(0,2.5,gnumber); resonator_excitation = Array{Float64}(gnumber) qubit_excitation = Array{Float64}(gnumber) groundstate_energy = Array{Float64}(gnumber) adaga = qeye(2)⊗(a_dagger(N)*a(N)) σpσm = (σ_plus*σ_minus)⊗qeye(N) for (i,g) in enumerate(couplings) H = gen_hamiltonian(N,ω_q,ω_r,g,rwa=false) evals, ekets = eig(H) groundstate_energy[i] = evals[1] ψ = ekets[:,1] resonator_excitation[i] = real(ψ'*adaga*ψ)[1] qubit_excitation[i] = real(ψ'*σpσm*ψ)[1] end plot([scatter(;x=couplings,y=resonator_excitation,name="Resonator Occupation"), scatter(;x=couplings,y=qubit_excitation,name="Qubit Occupation")], Layout(;xaxis=attr(title="Coupling"))) plot(scatter(;x=couplings,y=groundstate_energy/real(ħ)),Layout(;xaxis=attr(title="Coupling"),yaxis=attr(title="Energy (ħ)"))) H = gen_hamiltonian(N,ω_q,ω_r,couplings[end],rwa=false) _, ekets = eig(H) ρ = ekets[:,1]*ekets[:,1]' ρ_resonator = ptrace(ρ,[2,N],1) plot_densitymatrix(ρ_resonator) W = calc_wignerfunction_resonator(ρ_resonator,w_samples,maxdisp=7.5); plot_wignerfunction3d(W,maxdisp=7.5)