diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..3f31d183b87b61f6e60874a76330981b60b864b3
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,2 @@
+*__pycache__
+.vscode/
\ No newline at end of file
diff --git a/__main__.py b/__main__.py
new file mode 100644
index 0000000000000000000000000000000000000000..b3c2324bd3cd9b1f04a3c777e427197c6b2899b1
--- /dev/null
+++ b/__main__.py
@@ -0,0 +1,135 @@
+import numpy as np
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+import elements
+import graphics
+import matrices
+from math import sqrt
+
+print("Hi World!")
+
+## Tests
+'''print("###------ Portes anticommutantes ------")
+source = elements.Source(initial_vector = np.array([(3/4)*(1. + 0j), (4/5)*(1. + 0j)]))
+
+elements.Screen(previous = (source, "left")).show()
+
+BS1 = elements.BeamSplitter(previous = {"left": None, "right": (source, "left")}, matrix=matrices.HALF_MIRROR_LEFT_ACTIVE)
+
+elements.Screen(previous = (BS1, "left")).show("left")
+elements.Screen(previous = (BS1, "right")).show("right")
+
+X_left_1 = elements.Plate(previous = (BS1, "left"), jones_matrix=matrices.X)
+Z_left_2 = elements.Plate(previous = (X_left_1, "left"), jones_matrix = matrices.Z)
+
+
+Z_right_1 = elements.Plate(previous = (BS1, "right"), jones_matrix = matrices.Z)
+X_right_2 = elements.Plate(previous = (Z_right_1, "left"), jones_matrix = matrices.X)
+
+BS2 = elements.BeamSplitter(previous = {"left": (X_right_2, "left"), "right": (Z_left_2, "left")}, matrix = matrices.HALF_MIRROR_RIGHT_ACTIVE)
+
+screen_right = elements.Screen(previous = (BS2, "right"))
+screen_left = elements.Screen(previous = (BS2, "left"))
+
+screen_left.show("left")
+screen_right.show("right")'''
+
+'''print("### ------ Portes commutantes ------")
+
+source2 = elements.Source(initial_vector = np.array([1.+0j, 0.]))
+
+BS12 = elements.BeamSplitter(previous = {"left": (source2, "left"), "right": None}, matrix=matrices.SYM_BS_INACTIVE)
+
+elements.Screen(previous = (BS12, "left")).show()
+elements.Screen(previous = (BS12, "right")).show()
+
+Phase_left_12 = elements.Plate(previous = (BS12, "left"), jones_matrix=matrices.PHASE(np.pi/5))
+Z_left_22 = elements.Plate(previous = (Phase_left_12, "left"), jones_matrix = matrices.Z)
+
+
+Z_right_12 = elements.Plate(previous = (BS12, "right"), jones_matrix = matrices.Z)
+Phase_right_22 = elements.Plate(previous = (Z_right_12, "left"), jones_matrix = matrices.PHASE(np.pi/5))
+
+BS22 = elements.BeamSplitter(previous = {"left": (Z_left_22, "left"), "right": (Phase_right_22, "left")}, matrix = matrices.SYM_BS_INACTIVE)
+
+screen_right2 = elements.Screen(previous = (BS22, "right"))
+screen_left2 = elements.Screen(previous = (BS22, "left"))
+
+screen_left2.show("left")
+screen_right2.show("right")'''
+
+
+'''print("###------ Tests BeamSplitters ------")
+source = elements.Source(initial_vector = np.array([(3/5)*(1. + 0j), (4/5)*(1. + 0j)]))
+
+elements.Screen(previous = (source, "left")).show()
+
+BS1 = elements.BeamSplitter(previous = {"left": None, "right": (source, "left")}, matrix=matrices.HALF_MIRROR_LEFT_ACTIVE)
+
+elements.Screen(previous = (BS1, "left")).show("left")
+elements.Screen(previous = (BS1, "right")).show("right")'''
+
+'''print("###------ Portes classiques ------")
+source = elements.Source(initial_vector = np.array([1. + 0j, 0.]))
+
+elements.Screen(previous = (source, "left")).show()
+
+Nom =   {str(matrices.X) : "X",
+        str(matrices.Y) : "Y",
+        str(matrices.Z) : "Z",
+        str(matrices.H) : "H",
+        str(matrices.S) : "S",
+        str(matrices.T) : "T",}
+
+Portes =    [matrices.X,
+            matrices.Y,
+            matrices.Z,
+            matrices.H,
+            matrices.S,
+            matrices.T]
+
+for porte1 in Portes:
+    P1 = elements.Plate(previous = (source, "left"), jones_matrix=porte1)
+    elements.Screen(previous = (P1, "left")).show( Nom[str(porte1)] + "|ψ⟩")
+    for porte2 in Portes:
+        P2 = elements.Plate(previous = (P1, "left"), jones_matrix=porte2)
+        elements.Screen(previous = (P2, "left")).show(Nom[str(porte2)] + Nom[str(porte1)] + "|ψ⟩")'''
+
+print("###------ Tentative marche aléatoire ------")
+
+
+def marche(n):
+    Photons = [elements.Source(initial_vector = np.array([0., 1.+0j]))]
+    for i in range(n):
+        New_Photons = []
+        new_photons = dict()
+        for source in Photons:
+            PH = elements.Plate(previous = (source, "left"), jones_matrix=matrices.H)
+            BS = elements.BeamSplitter(previous = {"left": (PH, "left"), "right": None}, matrix=matrices.POLARISER_BS)
+            LEFT = elements.Path(previous = (BS, "left"), time = -1)
+            RIGHT = elements.Path(previous = (BS, "right"), time = 1)
+            Screen_left = elements.Screen(previous = (LEFT, "left"))
+            Screen_right = elements.Screen(previous = (RIGHT, "right"))
+            if Screen_left.temporality() in new_photons.keys():
+                new_photons[Screen_left.temporality()][0] += Screen_left.out()[0]
+                new_photons[Screen_left.temporality()][1] += Screen_left.out()[1]
+            else:
+                new_photons[Screen_left.temporality()] = np.array([Screen_left.out()[0], Screen_left.out()[1]])
+            if Screen_right.temporality() in new_photons.keys():
+                new_photons[Screen_right.temporality()][0] += Screen_right.out()[0] 
+                new_photons[Screen_right.temporality()][1] += Screen_right.out()[1]
+            else:
+                new_photons[Screen_right.temporality()] = np.array([Screen_right.out()[0], Screen_right.out()[1]])
+        Photons = [elements.Source(initial_vector = list(new_photons.values())[k] ,time = list(new_photons.keys())[k] ) for k in range(len(new_photons))]
+    repartition = dict()
+    for photon in Photons:
+        if photon.temporality('left') in repartition.keys():
+            repartition[photon.temporality('left')] += photon.out('left')[0]**2 + photon.out('left')[1]**2
+        else:
+            repartition[photon.temporality('left')] = photon.out('left')[0]**2 + photon.out('left')[1]**2
+    plt.plot(list(repartition.keys()),list(repartition.values()))
+    plt.title("Distribution de densité de probabilité")
+    plt.xlabel("Position")
+    plt.ylabel("Probabilité de présence")
+    plt.show()
+marche(100)
\ No newline at end of file
diff --git a/definition.py b/definition.py
deleted file mode 100644
index 4958740d3917f29e60675e4542d7cdf9d215cb17..0000000000000000000000000000000000000000
--- a/definition.py
+++ /dev/null
@@ -1,136 +0,0 @@
-from __future__ import (absolute_import, division, print_function, unicode_literals)
-from cmath import cos, sin, exp, polar, acos
-from math import pi
-import matplotlib as mpl
-from mpl_toolkits.mplot3d import Axes3D
-import numpy as np
-import matplotlib.pyplot as plt
-from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg
-
-mpl.rcParams['legend.fontsize'] = 10
-
-class Plate:
-    """type = ''LR, 'LP', 'CP'
-    Fast axis unchanged convention (pour un LR)
-    orientation = Right, Left (pour un CP)"""
-    def __init__(self, element, theta = 0, delta = 0, orientation = ""):
-        self.element = element
-        self.delta = delta
-        self.theta = theta
-        self.orientation = orientation
-
-class Photon:
-    def __init__(self, state0, state1):
-        self.state0 = [state0]
-        self.state1 = [state1]
-
-    def state(self):
-        print(self.state0[-1], "|0> + ", self.state1[-1], "|1>")
-
-    def round_state(self,decimal,i):
-        s_0r = round((self.state0[i]).real,decimal)
-        s_0i = round((self.state0[i]).imag,decimal)
-        s_1r = round((self.state1[i]).real,decimal)
-        s_1i = round((self.state1[i]).imag,decimal)
-        return(s_0r, s_0i, s_1r, s_1i)
-
-    def is_2D(self,i):
-        s_0r = round((self.state0[i]).real,3)
-        s_0i = round((self.state0[i]).imag,3)
-        s_1r = round((self.state1[i]).real,3)
-        s_1i = round((self.state1[i]).imag,3)
-        return (s_0r*s_0i == 0 and s_1r*s_1i == 0 and s_0r*s_1i == 0 and s_1r*s_0i == 0)
-
-    def pur(self,i):
-        rho0, phi0 = polar(self.state0[i])
-        rho1, phi1 = polar(self.state1[i])
-        return (2*acos(rho0), phi1-phi0)
-
-    def gate(self,plate):
-        alpha = self.state0[-1]
-        beta = self.state1[-1]
-        if plate.element == 'LP':
-            self.state0.append(alpha*cos(plate.theta)*cos(plate.theta) + beta*cos(plate.theta)*sin(plate.theta))
-            self.state1.append(alpha*cos(plate.theta)*sin(plate.theta) + beta*sin(plate.theta)*sin(plate.theta))
-        if plate.element == 'LR':
-            self.state0.append(alpha*(cos(plate.theta)*cos(plate.theta) + exp(1j*plate.delta)*sin(plate.theta)*sin(plate.theta)) + beta*(1-exp(1j*plate.delta))*cos(plate.theta)*sin(plate.theta))
-            self.state1.append(alpha*(1-exp(1j*plate.delta))*cos(plate.theta)*sin(plate.theta) + beta*(sin(plate.theta)*sin(plate.theta) + exp(1j*plate.delta)*cos(plate.theta)*cos(plate.theta)))
-        else:
-            if plate.orientation == "Right":
-                self.state0.append(1/2*(alpha + 1j*beta))
-                self.state1.append(1/2*(beta - 1j*alpha))
-            else:
-                self.state0.append(1/2*(alpha - 1j*beta))
-                self.state1.append(1/2*(beta + 1j*alpha))
-
-    def representation(self,i):
-
-        #Bloch Sphere
-
-        fig = plt.figure()
-        ax = fig.gca(projection='3d')
-        """thismanager = plt.get_current_fig_manager()
-        thismanager.window.SetPosition((500, 0))
-        thismanager.window.wm_geometry("+500+0")"""
-        ax.set_title('Step '+str(i))
-        theta = np.linspace(0, 2 * np.pi, 100)
-        z = np.zeros(100)
-        x = np.sin(theta)
-        y = np.cos(theta)
-        ax.plot(x, y, z, color = 'black', linestyle='dashed', linewidth=0.5, label='sphere')
-        ax.plot(y, z, x, color = 'black', linestyle='dashed', linewidth=0.5)
-        ax.plot(z, x, y, color = 'black', linestyle='dashed', linewidth=0.5)
-        ax.quiver(0, 0, 0, 0, 0, 1, color = 'black', arrow_length_ratio = 0.1)
-        ax.text(0, 0, 1.1, '|0>', color = 'black')
-        ax.quiver(0, 0, 0, 0, 0, -1, color = 'black', arrow_length_ratio = 0.1)
-        ax.text(0, 0, -1.1, '|1>', color = 'black')
-
-        if i>0:
-            theta, phi = self.pur(i-1)
-            ax.quiver(0, 0, 0, sin(theta)*cos(phi), sin(theta)*sin(phi), cos(theta), color = 'red', arrow_length_ratio = 0.1, label ='before')
-            ax.text(sin(theta)*cos(phi)+0.1, sin(theta)*sin(phi)+0.1, cos(theta)+0.1, 'before', color = 'red')
-
-        theta, phi = self.pur(i)
-        ax.quiver(0, 0, 0, sin(theta)*cos(phi), sin(theta)*sin(phi), cos(theta), color = 'green', arrow_length_ratio = 0.1, label ='after')
-        ax.text(sin(theta)*cos(phi)+0.1, sin(theta)*sin(phi)+0.1, cos(theta)+0.1, 'after', color = 'green')
-
-        ax.grid(False)
-        ax.axis(False)
-        ax.legend()
-
-        #2D representation
-
-        if self.is_2D(i):
-            fig2 = plt.figure()
-            ax_2D = fig2.add_subplot(111)
-            theta = np.linspace(0, 2 * np.pi, 100)
-            x = np.sin(theta)
-            y = np.cos(theta)
-            ax_2D.plot(x,y, color = 'black', linestyle='dashed', linewidth=0.5, label='circle')
-            ax_2D.quiver(*[0, 0], [0, 1], [1,0], scale = 1, scale_units = 'xy', color = 'black')
-            ax_2D.text(0, 1.1, '|0>', color = 'black')
-            ax_2D.text(1.1, 0, '|1>', color = 'black')
-            ax_2D.grid(False)
-            ax_2D.axis(False)
-            ax_2D.legend()
-            (s_0r, s_0i, s_1r, s_1i) = self.round_state(3,i)
-            if (s_0r != 0 or s_1r != 0):
-                ax_2D.quiver(*[0, 0], [s_1r], [s_0r], scale = 1, scale_units = 'xy', color = 'red', label = 'Phase nul')
-                ax_2D.text(s_1r +0.1, s_0r + 0.1, 'after', color = 'red')
-            else:
-                ax_2D.quiver(*[0, 0], [s_1i], [s_0i], scale = 1, scale_units = 'xy', color = 'red', label = 'Phase pi/2')
-                ax_2D.text(s_1i +0.1, s_0i + 0.1, 'after', color = 'red')
-
-            if i>0:
-                if self.is_2D(i-1):
-                    (s_0r, s_0i, s_1r, s_1i) = self.round_state(3,i-1)
-                    if (s_0r != 0 or s_1r != 0):
-                        ax_2D.quiver(*[0, 0], [s_1r], [s_0r], scale = 1, scale_units = 'xy', color = 'green', label = 'Phase nul')
-                        ax_2D.text(s_1r +0.1, s_0r + 0.1, 'before', color = 'green')
-                    else:
-                        ax_2D.quiver(*[0, 0], [s_1i], [s_0i], scale = 1, scale_units = 'xy', color = 'green', label = 'Phase pi/2')
-                        ax_2D.text(s_1i +0.1, s_0i + 0.1, 'before', color = 'green')
-
-
-        plt.show()
-
diff --git a/elements.py b/elements.py
new file mode 100644
index 0000000000000000000000000000000000000000..414560ae3228bcdf1ffbf4984100c53a461fa787
--- /dev/null
+++ b/elements.py
@@ -0,0 +1,131 @@
+import numpy as np
+
+### ------ Special elements ------
+class Source:
+    def __init__(self, initial_vector=np.array([1., 0.+0j]), time = 0): ## if no input state given, starts a |0>
+        self.photon = initial_vector
+        self.next = None
+        self.time = time
+
+    def out(self, side):
+        return self.photon
+    
+    def temporality(self, side):
+        return self.time
+
+
+class Screen:
+    def __init__(self, previous=(None, None)):
+        self.previous = previous
+
+    def show(self, name = 'ψ'):
+        previous_component, previous_component_output_side = self.previous
+        ket = previous_component.out(previous_component_output_side)
+        alpha = round(ket[0],3)
+        beta = round(ket[1],3)
+        print(f"|{name}⟩ = {alpha}|0⟩ + {beta}|1⟩")
+
+    def out(self):
+        previous_component, previous_component_output_side = self.previous
+        return previous_component.out(previous_component_output_side)
+
+    def temporality(self):
+        previous_component, previous_component_output_side = self.previous
+        return previous_component.temporality(previous_component_output_side) 
+
+
+
+### ------ Classic Elements ------
+class Element:
+    def __init__(self, previous): ## previous is either a single (previous_object, output) tuple, or a 2-dictionary of tuples
+        self.previous = previous  ## (as in {"left": (Object 1, "label of chosen output of object 1"), "right": (obj2, "label2")})
+        self.next = None
+
+class Path(Element):
+    def __init__(self, previous, time):
+        Element.__init__(self, previous)
+        self.time = time
+
+    def out(self, side):
+        previous_component, previous_component_output_side = self.previous
+        return previous_component.out(previous_component_output_side)
+    
+    def temporality(self, side):
+        previous_component, previous_component_output_side = self.previous
+        return previous_component.temporality(previous_component_output_side) + self.time
+
+class Plate(Element):
+    def __init__(self, previous, jones_matrix = np.eye(2, dtype='complex')):
+        Element.__init__(self, previous)
+        self.jones_matrix = jones_matrix
+        self.output_vector = None
+
+    def out(self, side):
+        previous_element = self.previous[0]
+        previous_element_output_side = self.previous[1]
+        self.output_vector = self.jones_matrix @ previous_element.out(previous_element_output_side) # matrix-vector product : jones matrix times state vector given by the out(side) method of the previous component
+        return self.output_vector
+    
+    def temporality(self, side):
+        previous_component, previous_component_output_side = self.previous
+        return previous_component.temporality(previous_component_output_side)
+
+class BeamSplitter(Element):
+    """matrix is a 4x4 complex ndarray, two first dimensions regard left input, two next are right. So it's blockwise defined."""
+    def __init__(self, previous, matrix = np.eye(4, dtype='complex')):
+        Element.__init__(self, previous)
+        self.matrix = matrix
+        self.output_vector = None
+
+    def out(self, side):
+        if self.previous["left"] is None:
+            left_input_state = np.zeros((2,), dtype="complex") # if no input on the left, then |\psi> = 0
+        else:
+            left_input_element = self.previous["left"][0]
+            left_input_element_output_side = self.previous["left"][1]
+            left_input_state = left_input_element.out(left_input_element_output_side)
+
+        if self.previous["right"] is None:
+            right_input_state = np.zeros((2,), dtype="complex") # if no input on the right, then |\psi> = 0
+        else:
+            right_input_element = self.previous["right"][0]
+            right_input_element_output_side = self.previous["right"][1]
+            right_input_state = right_input_element.out(right_input_element_output_side)
+
+        input_state = np.concatenate((left_input_state, right_input_state))
+        self.output_vector = self.matrix @ input_state
+
+        if side == "left":
+            ans, _ = np.split(self.output_vector, 2)
+
+        elif side == "right":
+            _, ans = np.split(self.output_vector, 2)
+
+        return ans
+
+    def temporality(self, side):
+        if self.previous["left"] is None and self.previous["right"] is None:
+            return 0
+
+        elif self.previous["left"] is None:
+            right_previous_component = self.previous["right"][0]
+            right_previous_component_output_side = self.previous["right"][1]
+            right_time = right_previous_component.temporality(right_previous_component_output_side)
+            return right_time
+
+        elif self.previous["right"] is None:
+            left_previous_component = self.previous["left"][0]
+            left_previous_component_output_side = self.previous["left"][1]
+            left_time = left_previous_component.temporality(left_previous_component_output_side)
+            return left_time
+        else:
+            right_previous_component = self.previous["right"][0]
+            right_previous_component_output_side = self.previous["right"][1]
+            right_time = right_previous_component.temporality(right_previous_component_output_side)
+            left_previous_component = self.previous["left"][0]
+            left_previous_component_output_side = self.previous["left"][1]
+            left_time = left_previous_component.temporality(left_previous_component_output_side)
+            if side == "left":
+                return left_time
+            elif side == "right":
+                return right_time
\ No newline at end of file
diff --git a/graphics.py b/graphics.py
new file mode 100644
index 0000000000000000000000000000000000000000..796eb8eb44ac917281ee8cd64cc89b5994996d2e
--- /dev/null
+++ b/graphics.py
@@ -0,0 +1 @@
+import matplotlib.pyplot as plt 
diff --git a/interface.py b/interface.py
deleted file mode 100644
index aaeb2d1a1e5b0a7f1325bf891bcf4e931e103dfd..0000000000000000000000000000000000000000
--- a/interface.py
+++ /dev/null
@@ -1,209 +0,0 @@
-from tkinter import *
-from copy import deepcopy
-import pygame as pyg
-from functools import partial
-from definition import Plate, Photon
-from math import pi
-from matplotlib.figure import Figure
-
-#Define the number of plates in our system
-def number_elements(root):
-    v = StringVar(root)
-    v.set(1)
-    entry = Entry(root, textvariable = v, width=2)
-    entry.config(font = ('Bahnschrift SemiLight','25'), bg = '#cccccc', fg = '#000000' )
-    title = Label(root, text = "Choose number of elements", font = ('Bahnschrift Semibold','30'), fg = '#000000', bg = '#ffffff', relief = 'raised')
-    title.grid(column = 2, row = 2)
-    entry.grid(column = 2, row = 3)
-    return v
-
-#Define the initial states of the photon
-def initial_photon(root):
-    title = Label(root, text = "Initial coefficients", font = ('Bahnschrift Semibold','30'), fg = '#000000', bg = '#ffffff', relief = 'raised')
-    title.grid(column = 5, row = 2, columnspan = 5)
-    state0 = StringVar(root)
-    state0.set(0)
-    entry0 = Entry(root, textvariable = state0, width=2)
-    entry0.config(font = ('Bahnschrift SemiLight','25'), bg = '#cccccc', fg = '#000000' )
-    entry0.grid(column = 5, row = 3)
-    label0 = Label(root, text = "|0>", font = ('Bahnschrift SemiLight','30'), fg = '#000000', bg = '#ffffff')
-    label0.grid(column = 6, row = 3)
-    labelplus = Label(root, text = "+", font = ('Bahnschrift SemiLight','30'), fg = '#000000', bg = '#ffffff')
-    labelplus.grid(column = 7, row = 3)
-    state1 = StringVar(root)
-    state1.set(0)
-    entry1 = Entry(root, textvariable = state1, width=2)
-    entry1.config(font = ('Bahnschrift SemiLight','25'), bg = '#cccccc', fg = '#000000' )
-    entry1.grid(column = 8, row = 3)
-    label1 = Label(root, text = "|1>", font = ('Bahnschrift SemiLight','30'), fg = '#000000', bg = '#ffffff')
-    label1.grid(column = 9, row = 3)
-    return state0, state1
-
-#Show the current properties of the plates while defining them
-def display(plate,window):
-    size = len(plate)
-    number_title = Label(window, text = 'Number', font = ('Bahnschrift Semibold','25'), fg = '#ffffff', bg = '#aaaaaa', relief= 'groove')
-    number_title.grid(column = 10, row = 0)
-    type_title = Label(window, text = 'Type', font = ('Bahnschrift Semibold','25'), fg = '#ffffff', bg = '#aaaaaa', relief= 'groove')
-    type_title.grid(column = 11, row = 0)
-    theta_title = Label(window, text = 'Theta', font = ('Bahnschrift Semibold','25'), fg = '#ffffff', bg = '#aaaaaa', relief= 'groove')
-    theta_title.grid(column = 12, row = 0)
-    delta_title = Label(window, text = 'Delta', font = ('Bahnschrift Semibold','25'), fg = '#ffffff', bg = '#aaaaaa', relief= 'groove')
-    delta_title.grid(column = 13, row = 0)
-    orientation_title = Label(window, text = 'Orientation', font = ('Bahnschrift Semibold','25'), fg = '#ffffff', bg = '#aaaaaa', relief= 'groove')
-    orientation_title.grid(column = 14, row = 0)
-    for i in range(size):
-        number = Label(window, text = str(i+1), font = ('Bahnschrift SemiBold Condensed','20'), fg = '#000000', bg = '#ffffff')
-        number.grid(column = 10, row = i+1)
-        str_element = StringVar(window)
-        str_element.set('0')
-        str_element = plate[i].element
-        current_element = Label(window, text = str_element, font = ('Bahnschrift SemiBold Condensed','20'), fg = '#000000', bg = '#ffffff')
-        current_element.grid(column = 11, row = i+1)
-        str_theta = StringVar(window)
-        str_theta.set('0')
-        str_theta = str(plate[i].theta*180/pi) 
-        current_theta = Label(window, text = str_theta, font = ('Bahnschrift SemiBold Condensed','20'), fg = '#000000', bg = '#ffffff')
-        current_theta.grid(column = 12, row = i+1)
-        str_delta = StringVar(window)
-        str_delta.set('0')
-        str_delta = str(plate[i].delta*180/pi) 
-        current_delta = Label(window, text = str_delta, font = ('Bahnschrift SemiBold Condensed','20'), fg = '#000000', bg = '#ffffff')
-        current_delta.grid(column = 13, row = i+1)
-        str_orientation = StringVar(window)
-        str_orientation.set('0')
-        str_orientation = plate[i].orientation 
-        current_orientation = Label(window, text = str_orientation, font = ('Bahnschrift SemiBold Condensed','20'), fg = '#000000', bg = '#ffffff')
-        current_orientation.grid(column = 14, row = i+1)
-
-#Defining the properties of the plates
-def define_elements(window,size):
-    Number = [i for i in range(1,size+1)]
-    Type = ('Linear Retarder', 'Linear Polarizer', 'Circular Polarizer')
-    Inverse_Type = {'Linear Retarder' : 'LR', 'Linear Polarizer' : 'LP', 'Circular Polarizer' : 'CP'}
-    Orientation = ('None', 'Right', 'Left')
-    number_title = Label(window, text = 'Number', font = ('Bahnschrift Semibold','30'), fg = '#000000', bg = '#ffffff', relief= 'raised')
-    number_title.grid(column = 0, row = 0)
-    type_title = Label(window, text = 'Type', font = ('Bahnschrift Semibold','30'), fg = '#000000', bg = '#ffffff', relief= 'raised')
-    type_title.grid(column = 1, row = 0)
-    theta_title = Label(window, text = 'Theta', font = ('Bahnschrift Semibold','30'), fg = '#000000', bg = '#ffffff', relief= 'raised')
-    theta_title.grid(column = 2, row = 0)
-    delta_title = Label(window, text = 'Delta', font = ('Bahnschrift Semibold','30'), fg = '#000000', bg = '#ffffff', relief= 'raised')
-    delta_title.grid(column = 3, row = 0)
-    orientation_title = Label(window, text = 'Orientation', font = ('Bahnschrift Semibold','30'), fg = '#000000', bg = '#ffffff', relief= 'raised')
-    orientation_title.grid(column = 4, row = 0)
-    global number
-    number = StringVar(window)
-    number.set(Number[0])
-    om_number = OptionMenu(window, number, *Number)
-    om_number.config(font = ('Bahnschrift SemiLight','25'), bg = '#cccccc', fg = '#000000')
-    om_number.grid(column = 0, row = 1)
-    global type
-    type = StringVar(window)
-    type.set(Type[0])
-    om_type = OptionMenu(window, type, *Type)
-    om_type.config(font = ('Bahnschrift SemiLight','25'), bg = '#cccccc', fg = '#000000' )
-    om_type.grid(column = 1, row = 1)
-    global orientation
-    orientation = StringVar(window)
-    orientation.set(Orientation[0])
-    om_orientation = OptionMenu(window, orientation, *Orientation)
-    om_orientation.config(font = ('Bahnschrift SemiLight','25'), bg = '#cccccc', fg = '#000000' )
-    om_orientation.grid(column = 4, row = 1)
-    global theta
-    theta = StringVar(window)
-    theta.set(0)
-    entry_theta = Entry(window, textvariable = theta, width=3)
-    entry_theta.config(font = ('Bahnschrift SemiLight','25'), bg = '#cccccc', fg = '#000000' )
-    entry_theta.grid(column = 2, row = 1)
-    global delta
-    delta = StringVar(window)
-    delta.set(0)
-    entry_delta = Entry(window, textvariable = delta, width=3)
-    entry_delta.config(font = ('Bahnschrift SemiLight','25'), bg = '#cccccc', fg = '#000000')
-    entry_delta.grid(column = 3, row = 1)
-    def give_plate():
-        global number
-        global type
-        global theta
-        global delta
-        global orientation
-        global plate
-        new_plate = Plate(Inverse_Type[type.get()],float(theta.get())*pi/180, float(delta.get())*pi/180, orientation.get())
-        print(new_plate.element, new_plate.theta, new_plate.delta, new_plate.orientation)
-        plate[int(number.get())-1] = new_plate
-        display(plate,window)
-    button = Button(window, fg="#ffffff",text= "Validate", command=give_plate, font = ('Bahnschrift SemiLight','20','bold'), relief = 'raised', bg = 'black')
-    button.grid(column = 4, row = 6)
-    button_quit = Button(window, text="Finish", fg="#ffffff", command=window.destroy, font = ('Bahnschrift SemiLight','20','bold'), relief = 'raised', bg = 'black')
-    button_quit.grid(column = 0, row = 6)
-    #return (Inverse_Type[type.get()],theta.get(), delta.get(), orientation.get())
-
-class Time:
-    def __init__(self, k, photon, window):
-        self.k = k
-        self.photon = photon
-        self.window = window
-    def bloch_sphere(self):
-        (self.photon).representation(self.k)
-    def do(self):
-        state_in_k = Button(self.window, text="State", fg="#ffffff", command=self.bloch_sphere, font = ('Bahnschrift SemiLight','20','bold'), relief = 'raised', bg = 'black')
-        state_in_k.grid(column = 2*self.k, row = 0, pady = 20)
-
-#Simulation of the optical path
-def optical_path(window, photon):
-    length = len(photon.state0)
-    for i in range(length):
-        time = Time(i,photon,window)
-        time.do()
-        if i != length-1:
-            plate_number = Label(window, text = 'Plate' + str(i+1), font = ('Bahnschrift Semibold','15'), fg = '#000000', bg = '#ffffff', relief = 'raised')
-            plate_number.grid(column = 2*i+1, row = 0)
-
-#Global software
-def graphical_grid_init():
-    root = Tk()
-    root.geometry("+150+20")
-    root.config(bg = '#ffffff')
-    root.title('PhotoniCS')
-    titre = Label(root, text = "PhotoniCS", font = ('Bahnschrift SemiBold Condensed','60'), fg = '#000000', bg = '#bbbbbb', relief= 'groove')
-    titre.grid(column = 4, row = 0)
-    button = Button(root, text="Quit", fg="#ffffff", command=quit, font = ('Bahnschrift SemiLight','20','bold'), relief = 'raised', bg = 'black')
-    button.grid(column = 0, row = 6)
-    global v
-    v = number_elements(root)
-    global plate
-    global state0
-    global state1
-    state0, state1 = initial_photon(root)
-    "pyg.mixer.init()"
-    def define():
-            #bg_music= pyg.mixer.Sound("bg.wav")
-            #bg_music.play()
-            #pyg.mixer.music.set_volume(0.5)
-            window = Toplevel()
-            window.config(bg = '#ffffff')
-            window.geometry('+100+350')
-            global plate
-            plate = [Plate('LR') for i in range(int(v.get()))]
-            define_elements(window,int(v.get()))
-    button = Button(root, fg="#ffffff",text= "Define", command=define, font = ('Bahnschrift SemiLight','20','bold'), relief = 'raised', bg = 'black')
-    button.grid(column = 2, row = 6)
-    def simulate():
-            #bg_music= pyg.mixer.Sound("bg.wav")
-            #bg_music.play()
-            #pyg.mixer.music.set_volume(0.5)
-            photon = Photon(float(state0.get()),float(state1.get()))
-            global plate
-            for i in range(int(v.get())):
-                photon.gate(plate[i])
-            window = Toplevel()
-            window.config(bg = '#ffffff')
-            optical_path(window,photon)
-            window.geometry('+150+20')
-            #pyg.mixer.quit()
-    button_simulate = Button(root, fg="#ffffff",text= "Simulate", command=simulate, font = ('Bahnschrift SemiLight','20','bold'), relief = 'raised', bg = 'black')
-    button_simulate.grid(column = 7, row = 6)
-    root.mainloop()
-
-graphical_grid_init()
diff --git a/matrices.py b/matrices.py
new file mode 100644
index 0000000000000000000000000000000000000000..8fb6268d9daae5ee3d6e57b2fc6b71b95e4b6749
--- /dev/null
+++ b/matrices.py
@@ -0,0 +1,80 @@
+import numpy as np
+### ------ ONE-PHOTON MATRICES ------
+### Pauli matrices
+
+X = np.array([[0, 1], [1, 0]], dtype="complex")
+Y = np.array([[0, 0.-1j], [0+1j, 0]], dtype = "complex")
+Z = np.array([[1, 0], [0, -1]], dtype = "complex")
+
+### Other One-qubit matrices
+
+H = (1/np.sqrt(2))*np.array([[1, 1],
+                             [1, -1]], dtype='complex')
+H_REV = (1/np.sqrt(2))*np.array([[-1, 1],
+                                 [1, 1]], dtype='complex')
+S = np.array([[1, 0], [0, 0+1j]], dtype='complex')
+T = np.array([[1, 0], [0, np.exp(complex(0, np.pi/8))]], dtype="complex")
+
+
+QUARTER_WAVE_RIGHT_CIRCULAR_RETARDER = (1/np.sqrt(2))*np.array([[1j, 1j],
+                                                                [-1j, 1j]])
+QUARTER_WAVE_LEFT_CIRCULAR_RETARDER = QUARTER_WAVE_RIGHT_CIRCULAR_RETARDER.T
+
+HALF_WAVE_RIGHT_CIRCULAR_RETARDER = 1j*Y
+HALF_WAVE_LEFT_CIRCULAR_RETARDER = 1j*Y.T
+
+### Matrix functions
+
+def PHASE(angle):
+    return np.array([[1, 0], [0, np.exp(complex(0, angle))]], dtype="complex")
+
+def LINEAR_POLARIZER(theta):
+    return np.array([[np.cos(theta)*np.cos(theta), np.cos(theta)*np.sin(theta)],
+                     [np.cos(theta)*np.sin(theta), np.sin(theta)*np.sin(theta)]])
+
+def ELLIPTICAL_POLARIZER(epsilon, psi):
+    J11 = np.cos(psi)**2 + (epsilon*np.sin(psi))**2
+    J12 = 1j*epsilon + (epsilon**2 - 1)*np.cons(psi)*np.sin(psi)
+    J22 = (epsilon*np.cos(psi))**2 + np.sin(psi)**2
+    return (1/(1+epsilon**2))*np.array([[J11, J12],
+                                        [-J12, J22]])
+
+def LINEAR_RETARDER(delta, theta): ### symmetric phase convention
+    L11 = np.exp(-1j*delta/2)*(np.cos(theta)**2) + np.exp(1j*delta/2)*(np.sin(theta)**2)
+    L12 =-1j*np.sin(delta/2)*np.sin(2*theta)
+    L22 = np.exp(1j*delta/2)*(np.cos(theta)**2) + np.exp(-1j*delta/2)*(np.sin(theta)**2)
+    return np.array([[L11, L12],
+                     [L12, L22]])
+
+
+### ------ TWO-PHOTON MATRICES ------
+
+SYM_BS_AMPLITUDE_MATRIX = (1/np.sqrt(2))*np.array([[1j, 1],
+                                                   [1, 1j]])
+
+HALF_MIRROR_LEFT_PLATE_AMPLITUDE_MATRIX = H_REV
+HALF_MIRROR_RIGHT_PLATE_AMPLITUDE_MATRIX = H
+
+SYM_BS_INACTIVE = np.kron(SYM_BS_AMPLITUDE_MATRIX, np.eye(2))
+HALF_MIRROR_LEFT_INACTIVE = np.kron(HALF_MIRROR_LEFT_PLATE_AMPLITUDE_MATRIX, np.eye(2))
+HALF_MIRROR_RIGHT_INACTIVE = np.kron(HALF_MIRROR_RIGHT_PLATE_AMPLITUDE_MATRIX, np.eye(2))
+
+SYM_BS_ACTIVE = (1/np.sqrt(2))*np.array([[-1j, 0,   1,   0],
+                                           [0,  1j,   0,   1],
+                                           [1,   0, -1j,   0],
+                                           [0,   1,   0,  1j]], dtype='complex')
+
+HALF_MIRROR_LEFT_ACTIVE = (1/np.sqrt(2))*np.array([[-1,  0,   1,   0],
+                                                   [0,   1,   0,   1],
+                                                   [1,   0,   1,   0],
+                                                   [0,   1,   0,  -1]], dtype='complex')
+
+HALF_MIRROR_RIGHT_ACTIVE = (1/np.sqrt(2))*np.array([[1,   0,   1,   0],
+                                                    [0,  -1,   0,   1],
+                                                    [1,   0,  -1,   0],
+                                                    [0,   1,   0,   1]], dtype='complex')
+
+POLARISER_BS = np.array([[1,   0,   0,   0],
+                        [0,  0,   1,   0],
+                        [0,   0,  0,   1],
+                        [0,   1,   0,   0]], dtype='complex')
diff --git a/requirements.txt b/requirements.txt
new file mode 100644
index 0000000000000000000000000000000000000000..f05fbbe6ada5640b7ca46df3ac2319cc42c35ef0
--- /dev/null
+++ b/requirements.txt
@@ -0,0 +1 @@
+scipy==1.5.1