diff --git a/matrices.py b/matrices.py
index 3976ce0adfd851dc38fc9c77390a9b0a5d11690f..c72e1e40ee7674516c14f383bf6d9b7387a858b5 100644
--- a/matrices.py
+++ b/matrices.py
@@ -1,5 +1,5 @@
 import numpy as np
-
+### ------ ONE-PHOTON MATRICES ------
 ### Pauli matrices
 
 X = np.array([[0, 1], [1, 0]], dtype="complex")
@@ -8,7 +8,8 @@ Z = np.array([[1, 0], [0, -1]], dtype = "complex")
 
 ### Other One-qubit matrices
 
-H = np.array([[1, 1],[1, -1]], dtype='complex')
+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")
 
@@ -42,3 +43,31 @@ def LINEAR_RETARDER(delta, theta): ### symmetric phase convention
     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_INACTIVE = (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_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')