from qutip import *
from pylab import *
%pylab inline

up=basis(2,0)
dn=basis(2,1)

rpp=up*up.dag()
rpm=up*dn.dag()
rmp=rpm.trans()
rmm=dn*dn.dag()

def depol(rho,p):
    K0=sqrt(1.0-p)*qeye(2)
    K1=sqrt(p/3.0)*sigmax()
    K2=sqrt(p/3.0)*sigmay()
    K3=sqrt(p/3.0)*sigmaz()
    x=K0*rho*K0.dag()+K1*rho*K1.dag()+K2*rho*K2.dag()+K3*rho*K3.dag()
    return x

p=0.8
K0=sqrt(1.0-p)*qeye(2)
K1=sqrt(p/3.0)*sigmax()
K2=sqrt(p/3.0)*sigmay()
K3=sqrt(p/3.0)*sigmaz()
K0*K0.dag()+K1*K1.dag()+K2*K2.dag()+K3*K3.dag()

def damp(rho,p):
    K0=rpp+sqrt(1.0-p)*rmm
    K1=sqrt(p)*rpm
    x=K0*rho*K0.dag()+K1*rho*K1.dag()
    return x

def flip(rho,p):
    K0=sqrt(p)*sigmaz()
    K1=sqrt(1.0-p)*qeye(2)
    x=K0*rho*K0.dag()+K1*rho*K1.dag()
    return x

psi=(up+dn).unit()

rho=ket2dm(psi)

plist=linspace(0,1.0,10.0)
slist=zeros(len(plist))

for i1 in range(len(plist)):
    p=plist[i1]
    rho1=depol(rho,p)
    slist[i1]=entropy_vn(rho1)

plot(plist,slist,'-',label='', linewidth=4)
xlabel('p',fontsize=20)
ylabel('SvN',fontsize=20)
show()

BB0=(tensor(rpp,rmm)+tensor(rpm,rmp)+tensor(rmp,rpm)+tensor(rmm,rpp))/2.0

p=0.4
rpp1=depol(rpp,p)
rpm1=depol(rpm,p)
rmp1=depol(rmp,p)
rmm1=depol(rmm,p)

BB0_E=(tensor(rpp1,rmm)+tensor(rpm1,rmp)+tensor(rmp1,rpm)+tensor(rmm1,rpp))/2.0 

concurrence(BB0_E)