Julia Error

punchout.jl

# Punchout measurement example
using Schrodinger

const n = 4                # number of qubit states
const N = 5                # number of cavity fock states

function JC_gen(fc=5, Ej=14, Ec=0.4, ϕ=0)
    ωc  = fc * 2π          # cavity frequency
    Ej′ = Ej*abs(cos(π*ϕ)) # Biased Ej
    ωq  = √(8*Ej′*Ec)*2π   # transmon "base" frequency
                           # real freq is √(8*Ej′*Ec)-Ec for 0-1 transition
    χ   = Ec/12 * 2π       # transmon anharmonicity
    g   = 0.2 * 2π         # coupling strength
    γ   = 0.044            # atom dissipation rate
    κ   = 40E-6            # cavity dissipation rate

    a = qeye(n) ⊗ destroy(N)
    b = destroy(n) ⊗ qeye(N)

    H0 = ωq*b'*b - χ*(b+b')^4 + ωc*a'*a + 0.5*g*(a'*b + a*b')
    c_ops = [sqrt(γ)*a, sqrt(κ)*b]

    return H0, a, c_ops
end

drivem(t,p) = p[1]*cis(+p[2]*t)
drivep(t,p) = p[1]*cis(-p[2]*t)

function run_steady(nϵ, nf, samples=16)
    g = Density(basis(n,0) ⊗ basis(N,0))
    ϵ_drive = logspace(-4,-1,nϵ)
    f_drive = linspace(4.8,5.2,nf)
    # Generate ops
    H0 , a, c_ops = JC_gen()
    N_op = qeye(n)⊗numberop(N)
    I_op = qeye(n)⊗(a+a')/sqrt(2)
    Q_op = qeye(n)⊗(a-a')*1im/sqrt(2)
    res = map(Base.product(ϵ_drive, f_drive)) do x
        ϵ, f = x; T = 1/f
        H = [H0, (a,drivem,[2π*ϵ,2π*f]), (a',drivep,[2π*ϵ,2π*f])]
        U = LindbladProp(H,c_ops,T,samples)
        return psteady(U,g,1000*5,[I_op,Q_op,N_op])
    end
    return ϵ_drive,f_drive,res
end

resB = run_steady(2,3)