from qutip import *
from pylab import *
%pylab inline

def mno(a,b):
    return b*(a*b)

up=basis(2,0)
dn=basis(2,1)
s0=qeye(2)
sx=sigmax()
sy=sigmay()
sz=sigmaz()

up

sx

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

rpp

rmp

BB0=(tensor(rpp,rmm)+tensor(rpm,rmp)+tensor(rmp,rpm)+tensor(rmm,rpp))/2.0
BB1=(tensor(rpp,mno(rmm,sx))+tensor(rpm,mno(rmp,sx))+tensor(rmp,mno(rpm,sx))+tensor(rmm,mno(rpp,sx)))/2.0
BB2=(tensor(rpp,mno(rmm,sy))+tensor(rpm,mno(rmp,sy))+tensor(rmp,mno(rpm,sy))+tensor(rmm,mno(rpp,sy)))/2.0
BB3=(tensor(rpp,mno(rmm,sz))+tensor(rpm,mno(rmp,sz))+tensor(rmp,mno(rpm,sz))+tensor(rmm,mno(rpp,sz)))/2.0

def telep(Q,R):
    p0=(BB0*R).tr()
    p1=(BB1*R).tr()
    p2=(BB2*R).tr()
    p3=(BB3*R).tr()
    Out=p0*mno(Q,s0)+p1*mno(Q,sx)+p2*mno(Q,sy)+p3*mno(Q,sz)
    return Out

p=0.0
R=(1.0-p)*BB0+p/4.0*tensor(qeye(2),qeye(2))

tpp=telep(rpp,R)
tpm=telep(rpm,R)
tmp=telep(rmp,R)
tmm=telep(rmm,R)

theta=pi/2.0
phi=0.0
In=(cos(theta/2)**2)*rpp+1/2.0*sin(theta)*exp(1j*phi)*rpm+(sin(theta/2)**2)*rmm+1/2.0*sin(theta)*exp(-1j*phi)*rmp

Out=(cos(theta/2)**2)*tpp+1/2.0*sin(theta)*exp(1j*phi)*tpm+(sin(theta/2)**2)*tmm+1/2.0*sin(theta)*exp(-1j*phi)*tmp

Out-mno(In,s0) # R=B0
#Out-mno(In,sx) # R=B1
#Out-mno(In,sy) # R=B2
#Out-mno(In,sz) # R=B3

F=(In*Out).tr()
F

plist=linspace(0.0,1.0,20.0)
flist=zeros(len(plist))
theta=pi/4.0
phi=0.0

for i1 in xrange(len(plist)):
    p=plist[i1]
    R=(1.0-p)*BB0+p/4.0*tensor(qeye(2),qeye(2))
    tpp=telep(rpp,R)
    tpm=telep(rpm,R)
    tmp=telep(rmp,R)
    tmm=telep(rmm,R)
    In=(cos(theta/2)**2)*rpp+1/2.0*sin(theta)*exp(1j*phi)*rpm+(sin(theta/2)**2)*rmm+1/2.0*sin(theta)*exp(-1j*phi)*rmp
    Out=(cos(theta/2)**2)*tpp+1/2.0*sin(theta)*exp(1j*phi)*tpm+(sin(theta/2)**2)*tmm+1/2.0*sin(theta)*exp(-1j*phi)*tmp
    flist[i1]=(In*Out).tr()

plot(plist,flist,'-.r',label='$\theta=\pi/4$', linewidth=4)
xlabel('p',fontsize=20)
ylabel('F',fontsize=20)
show()