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 dav_pp(A,G,p,W,t):
    xx=(1-(1-p)*(1-exp(-A*t)))*rpp+(1-p)*(1-exp(-A*t))*rmm
    return xx
def dav_pm(A,G,p,W,t):
    xx=exp(-1j*W*t)*exp(-G*t)*rpm
    return xx
def dav_mp(A,G,p,W,t):
    xx=exp(1j*W*t)*exp(-G*t)*rmp
    return xx
def dav_mm(A,G,p,W,t):
    xx=(1-((1-p)*(1-exp(-A*t)))*p/(1-p))*rmm+(1-p)*(1-exp(-A*t))*p/(1-p)*rpp
    return xx

def davis(rho,A,G,p,W,t):
    #
    dpp=dav_pp(A,G,p,W,t)
    dpm=dav_pm(A,G,p,W,t)
    dmp=dav_mp(A,G,p,W,t)
    dmm=dav_mm(A,G,p,W,t)
    #
    xpp=(rpp*rho.dag()).tr()
    xmp=(rmp*rho.dag()).tr()
    xpm=(rpm*rho.dag()).tr()
    xmm=(rmm*rho.dag()).tr()
    out=xpp*dpp+xmm*dmm+xmp*dmp+xpm*dpm
    return out
    

A=1.0
G=1.0
p=1.0/2.0
q=1.0/2.0
W=0.0

vv=(up+dn).unit()
rho=ket2dm(vv)
rho

tlist=linspace(0.0,3.0,20.0)
outlist=zeros(len(tlist))

for i1 in range(len(tlist)):
    t=tlist[i1]
    out=davis(rho,A,G,p,W,t)
    outlist[i1]=entropy_vn(out)

plot(tlist,outlist,'-.r',label='xx', linewidth=4)
xlabel('t',fontsize=20)
ylabel('S_vN',fontsize=20)
show()

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

A=1.0
G=1.0
p=1.0/2.0
q=1.0/2.0
W=0.0

t=0.1
dpp=dav_pp(A,G,p,W,t)
dpm=dav_pm(A,G,p,W,t)
dmp=dav_mp(A,G,p,W,t)
dmm=dav_mm(A,G,p,W,t)

Env=(tensor(dpp,rpp)+tensor(dpm,rpm)+tensor(dmp,rmp)+tensor(dmm,rmm))/2.0

concurrence(Env)

tlist=linspace(0.0,3.0,20.0)
outlist=zeros(len(tlist))
for i1 in range(len(tlist)):
    t=tlist[i1]
    #
    dpp=dav_pp(A,G,p,W,t)
    dpm=dav_pm(A,G,p,W,t)
    dmp=dav_mp(A,G,p,W,t)
    dmm=dav_mm(A,G,p,W,t)
    #
    out=(tensor(dpp,dpp)+tensor(dpm,dpm)+tensor(dmp,dmp)+tensor(dmm,dmm))/2.0
    outlist[i1]=concurrence(out)

plot(tlist,outlist,'-.r',label='xx', linewidth=4)
xlabel('t',fontsize=20)
ylabel('C',fontsize=20)
show()