diff --git a/__main__.py b/__main__.py
index bc5b490ae5905afbdfbaa55d9d0b1ccf9deb8443..873ba8c5c58e76d133407faff3d0bda7d658fd9f 100644
--- a/__main__.py
+++ b/__main__.py
@@ -1,4 +1,6 @@
 import numpy as np
+import matplotlib as mpl
+import matplotlib.pyplot as plt
 import elements
 import graphics
 
diff --git a/matrices.py b/matrices.py
new file mode 100644
index 0000000000000000000000000000000000000000..3976ce0adfd851dc38fc9c77390a9b0a5d11690f
--- /dev/null
+++ b/matrices.py
@@ -0,0 +1,44 @@
+import numpy as np
+
+### 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 = 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]])