%run init.ipynb

rx,ry,rz = symbols('r_x r_y r_z')
rho1qb = (1/2)*(id(2) + rx*pauli(1) + ry*pauli(2) + rz*pauli(3))
rho1qb

pI,pX,pZ,pY = symbols('p_I p_X p_Z p_Y')
K0 = sqrt(pI)*id(2); K1 = sqrt(pX)*pauli(1); K2 = sqrt(pZ)*pauli(3); K3 = sqrt(pY)*pauli(2)

rho_p = K0*rho1qb*K0 + K1*rho1qb*K1 + K2*rho1qb*K2 + K3*rho1qb*K3
simplify(rho_p)

p = symbols('p')
pI = (1+3*p)/4; pX = (1-p)/4; pY = pX; pZ = pX
rho_p = K0*rho1qb*K0 + K1*rho1qb*K1 + K2*rho1qb*K2 + K3*rho1qb*K3
simplify(rho_p)

import numpy as np
from matplotlib import pyplot as plt

p = np.arange(0,1.1,0.1) # PD channel
C = np.sqrt(1-p)
plt.plot(p,C)
plt.show()

from sympy import *

N,p = symbols('N p')
K0 = sqrt(1-N)*Matrix([[1,0],[0,sqrt(1-p)]])
K1 = sqrt(1-N)*Matrix([[0,sqrt(p)],[0,0]])
K2 = sqrt(N)*Matrix([[sqrt(1-p),0],[0,1]])
K3 = sqrt(N)*Matrix([[0,0],[sqrt(p),0]])
#K0

r00,r01,r10,r11 = symbols('r_{00} r_{01} r_{10} r_{11}')#; r00
rho = Matrix([[r00,r01],[r10,r11]])#; rho

rho_gad = K0*rho*K0 + K1*rho*K1.T + K2*rho*K2 + K3*rho*K3.T
simplify(rho_gad)

