from qutip import *
from pylab import *
%pylab inline

H=sigmaz()

D=sqrt(0.1)*sigmax()

psi00=(basis(2,0)+basis(2,1)).unit()
psi0=ket2dm(psi00)

tlist=linspace(0,5.0,10.0)
outlist=zeros(len(tlist))

qq=mesolve(H, psi0, tlist, [D], [])

for i1 in range(len(tlist)):
    rho=qq.states[i1]
    outlist[i1]=entropy_vn(rho)

plot(tlist,outlist,'-',label='', linewidth=4)
xlabel('t',fontsize=20)
ylabel('SvN',fontsize=20)
show()

N=10

H=create(N)*destroy(N)

D=sqrt(0.1)*destroy(N)

psi0=fock_dm(N,2)

tlist=linspace(0,5.0,10.0)
outlist=zeros(len(tlist))

qq=mesolve(H, psi0, tlist, [D], [])

for i1 in range(len(tlist)):
    rho=qq.states[i1]
    outlist[i1]=(rho*num(N)).tr()

plot(tlist,outlist,'-',label='', linewidth=4)
xlabel('t',fontsize=20)
ylabel('<N>',fontsize=20)
show()

H=tensor(sigmaz(),sigmaz())

d1=sqrt(0.1)*sigmax()
D1=tensor(qeye(2),d1)
D2=tensor(d1,qeye(2))

psi00=(tensor(basis(2,0),basis(2,0))+tensor(basis(2,1),basis(2,1))).unit()
psi0=ket2dm(psi00)

tlist=linspace(0,5.0,10.0)
outlist=zeros(len(tlist))

qq=mesolve(H, psi0, tlist, [D1,D2], [])

for i1 in range(len(tlist)):
    rho=qq.states[i1]
    outlist[i1]=concurrence(rho)

plot(tlist,outlist,'-',label='', linewidth=4)
xlabel('t',fontsize=20)
ylabel('Concurrence',fontsize=20)
show()