code for matrices for tomography of two entangled photons

tomographyTwoPhotons.py

import numpy
import math

knownProbNames = ["HH","HV","VV","VH","RH","RV","DV","DH","DR","DD","RD","HD","VD","VL","HL","RL"]

StateLetters = ["D","A","R","L","H","V"]
StateVectors = [
                    1/math.sqrt(2)*numpy.array([[1, 1]]), 
                    1/math.sqrt(2)*numpy.array([[1, -1]]),
                    1/math.sqrt(2)*numpy.array([[1, -1j]]), 
                    1/math.sqrt(2)*numpy.array([[1, 1j]]),
                    numpy.array([[1, 0]]), 
                    numpy.array([[0, 1]])
                ]              

pauliMatrices = numpy.array((
                                ((1, 0),
                                 (0, 1)),
                                 
                                ((0, 1),
                                (1, 0)),
                                
                                ((0, -1j),
                                (1j, 0)),
                                
                                ((1, 0),
                                (0, -1)),
                            ))

def getKroneckerOfPauli(i1, i2):
    return numpy.kron(pauliMatrices[i1], pauliMatrices[i2])

def getStateVectorFromNameString(Statename):
    firstStateIdx = StateLetters.index(Statename[0])
    secondStateIdx = StateLetters.index(Statename[1])
    return numpy.outer(StateVectors[firstStateIdx], StateVectors[secondStateIdx]).reshape((4,1))
    
def getB(PsiStateName, sigma16idx):
    statePhiVec4 = getStateVectorFromNameString(PsiStateName)
    statePhiHerConj = getStateVectorFromNameString(PsiStateName).conj().T
    gamma = getKroneckerOfPauli(sigma16idx%4,sigma16idx//4)
    return statePhiHerConj.dot(gamma.dot(statePhiVec4))[0][0]

def getMatBinv():
    bRet = numpy.empty((16, 16), dtype=complex)
    for currentStateIdx, currentStateName in enumerate(knownProbNames):
        for sigmaIdx in range(16):
            bRet[currentStateIdx][sigmaIdx] = getB(currentStateName, sigmaIdx)
    return numpy.linalg.inv(bRet)

def getGammaVec():
    gammaRet = [getKroneckerOfPauli(sigma16idx%4, sigma16idx//4) for sigma16idx in range(16)]
    return numpy.array([gammaRet]).reshape((16,4,4))

Binv = getMatBinv()
def getMmatrixForProb(knownProbIdx):
    matReturn = numpy.zeros((4, 4), dtype=complex)
    for sigma16idx in range(16):       
        matReturn += getKroneckerOfPauli(sigma16idx%4, sigma16idx//4)*Binv[sigma16idx][knownProbIdx]
    return matReturn
    
for stateIdx, stateName in enumerate(knownProbNames):
    print(f"Matrix coefficient for measurement {stateName}")
    print(f"{getMmatrixForProb(stateIdx)} \n")