%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from IPython.display import Image

from qutip import *

Image(filename='images/qobj.png')

q = Qobj([[1], [0]])

q

# the dimension, or composite Hilbert state space structure
q.dims

# the shape of the matrix data representation
q.shape

# the matrix data itself. in sparse matrix format. 
q.data

# get the dense matrix representation
q.full()

# some additional properties
q.isherm, q.type 

sy = Qobj([[0,-1j], [1j,0]])  # the sigma-y Pauli operator

sy

sz = Qobj([[1,0], [0,-1]]) # the sigma-z Pauli operator

sz

# some arithmetic with quantum objects

H = 1.0 * sz + 0.1 * sy

print("Qubit Hamiltonian = \n")
H

# The hermitian conjugate
sy.dag()

# The trace
H.tr()

# Eigen energies
H.eigenenergies()

# Fundamental basis states (Fock states of oscillator modes)

N = 2 # number of states in the Hilbert space
n = 1 # the state that will be occupied

basis(N, n)    # equivalent to fock(N, n)

fock(4, 2) # another example

# a coherent state
coherent(N=10, alpha=1.0)

# a fock state as density matrix
fock_dm(5, 2) # 5 = hilbert space size, 2 = state that is occupied

# coherent state as density matrix
coherent_dm(N=8, alpha=1.0)

# thermal state
n = 1 # average number of thermal photons
thermal_dm(8, n)

# Pauli sigma x
sigmax()

# Pauli sigma y
sigmay()

# Pauli sigma z
sigmaz()

#  annihilation operator

destroy(N=8) # N = number of fock states included in the Hilbert space

# creation operator

create(N=8) # equivalent to destroy(5).dag()

# the position operator is easily constructed from the annihilation operator
a = destroy(8)

x = a + a.dag()

x

def commutator(op1, op2):
    return op1 * op2 - op2 * op1

a = destroy(5)

commutator(a, a.dag())

x =       (a + a.dag())/sqrt(2)
p = -1j * (a - a.dag())/sqrt(2)

commutator(x, p)

commutator(sigmax(), sigmay()) - 2j * sigmaz()

-1j * sigmax() * sigmay() * sigmaz()

sigmax()**2 == sigmay()**2 == sigmaz()**2 == qeye(2)

sz1 = tensor(sigmaz(), qeye(2))

sz1

psi1 = tensor(basis(N,1), basis(N,0)) # excited first qubit
psi2 = tensor(basis(N,0), basis(N,1)) # excited second qubit

sz1 * psi1 == psi1 # this should not be true, because sz1 should flip the sign of the excited state of psi1

sz1 * psi2 == psi2 # this should be true, because sz1 should leave psi2 unaffected

sz2 = tensor(qeye(2), sigmaz())

sz2

tensor(sigmax(), sigmax())

epsilon = [1.0, 1.0]
g = 0.1

sz1 = tensor(sigmaz(), qeye(2))
sz2 = tensor(qeye(2), sigmaz())

H = epsilon[0] * sz1 + epsilon[1] * sz2 + g * tensor(sigmax(), sigmax())

H

wc = 1.0 # cavity frequency
wa = 1.0 # qubit/atom frenqency
g = 0.1  # coupling strength

# cavity mode operator
a = tensor(destroy(5), qeye(2))

# qubit/atom operators
sz = tensor(qeye(5), sigmaz())   # sigma-z operator
sm = tensor(qeye(5), destroy(2)) # sigma-minus operator

# the Jaynes-Cumming Hamiltonian
H = wc * a.dag() * a - 0.5 * wa * sz + g * (a * sm.dag() + a.dag() * sm)

H

a = tensor(destroy(3), qeye(2))
sp = tensor(qeye(3), create(2))

a * sp

tensor(destroy(3), create(2))

# Hamiltonian
H = sigmax()

# initial state
psi0 = basis(2, 0)

# list of times for which the solver should store the state vector
tlist = np.linspace(0, 10, 100)

result = mesolve(H, psi0, tlist, [], [])

result

len(result.states)

result.states[-1] # the finial state

expect(sigmaz(), result.states[-1])

expect(sigmaz(), result.states)

fig, axes = plt.subplots(1,1)

axes.plot(tlist, expect(sigmaz(), result.states))

axes.set_xlabel(r'$t$', fontsize=20)
axes.set_ylabel(r'$\left<\sigma_z\right>$', fontsize=20);

result = mesolve(H, psi0, tlist, [], [sigmax(), sigmay(), sigmaz()])

fig, axes = plt.subplots(1,1)

axes.plot(tlist, result.expect[2], label=r'$\left<\sigma_z\right>$')
axes.plot(tlist, result.expect[1], label=r'$\left<\sigma_y\right>$')
axes.plot(tlist, result.expect[0], label=r'$\left<\sigma_x\right>$')

axes.set_xlabel(r'$t$', fontsize=20)
axes.legend(loc=2);

w = 1.0               # oscillator frequency
kappa = 0.1           # relaxation rate
a = destroy(10)       # oscillator annihilation operator
rho0 = fock_dm(10, 5) # initial state, fock state with 5 photons
H = w * a.dag() * a   # Hamiltonian

# A list of collapse operators
c_ops = [sqrt(kappa) * a]

tlist = np.linspace(0, 50, 100)

# request that the solver return the expectation value of the photon number state operator a.dag() * a
result = mesolve(H, rho0, tlist, c_ops, [a.dag() * a]) 

fig, axes = plt.subplots(1,1)
axes.plot(tlist, result.expect[0])
axes.set_xlabel(r'$t$', fontsize=20)
axes.set_ylabel(r"Photon number", fontsize=16);

from qutip.ipynbtools import version_table

version_table()