diff --git a/tests/qse/magnetic_test.py b/tests/qse/magnetic_test.py
index 1929454..7825fd5 100644
--- a/tests/qse/magnetic_test.py
+++ b/tests/qse/magnetic_test.py
@@ -9,6 +9,7 @@
 
 
 def random_state(hdim):
+    "Outputs a random normalised state"
     statevector = np.random.rand(hdim) + 1j * np.random.rand(hdim)
     statevector /= np.linalg.norm(statevector)
     return statevector
@@ -27,12 +28,44 @@ def expval(op, state):
     return (np.conj(state) @ (op @ state[:, None])).item()
 
 
+def ij_operator(n, i, j, op_i, op_j):
+    ops = [np.eye(2) for _ in range(n)]
+    ops[i] = op_i
+    ops[j] = op_j
+    return kron_list(ops)
+
+
+def get_sisj_algebraic(i, j, statevector):
+    n = int(np.log2(len(statevector)))
+    term_1 = expval(ij_operator(n, i, j, pauli_x, pauli_x), statevector)
+    term_2 = expval(ij_operator(n, i, j, pauli_y, pauli_y), statevector)
+    term_3 = expval(ij_operator(n, i, j, pauli_z, pauli_z), statevector)
+    return (1 / 3) * (
+        term_1 + term_2 + term_3
+    ).real  # Normalise it to 1 for structure factor
+
+
+@pytest.mark.parametrize("n_qubits, hsize", [(4, 2**4), (3, 2**3 - 1)])
+def test_get_basis_shape(n_qubits, hsize):
+    """Check that get_basis outputs a n_qubits * h_size shape"""
+    ibasis = qse.magnetic.get_basis(n_qubits, hsize)
+    assert ibasis.shape == (hsize, n_qubits)
+
+
+@pytest.mark.parametrize("hdim", [1, 2, 3, 4, 5])
+def test_spin_values_are_less_than_one(hdim):
+    "Test that the absolute values of the computed spins is less than 1"
+    statevector = random_state(2**hdim)
+    ibasis = qse.magnetic.get_basis(hdim, 2**hdim)
+    spins = qse.magnetic.get_spins(statevector, hdim, ibasis)
+    assert np.all(np.abs(spins) <= 1)
+
+
 def test_spin_two_qubits():
     """Check get_spins works for two qubits."""
     n = 2
     statevector = random_state(2**n)
     spin_qse = qse.magnetic.get_spins(statevector, n)
-
     op_q0 = [np.kron(p, np.eye(2)) for p in [pauli_x, pauli_y, pauli_z]]
     expval_q0 = [expval(op, statevector) for op in op_q0]
     op_q1 = [np.kron(np.eye(2), p) for p in [pauli_x, pauli_y, pauli_z]]
@@ -66,3 +99,16 @@ def test_number_operator_random(n):
 
     n_op = qse.magnetic.get_number_operator(psi, n)
     assert np.allclose(n_op, np_n_op)
+
+
+@pytest.mark.parametrize("n", [2, 3, 4])
+def test_get_sisj(n):
+    """Tests the get_sisj function."""
+    psi = random_state(2**n)
+    sisj = qse.magnetic.get_sisj(psi, n, qse.magnetic.get_basis(n, 2**n))
+    sisj_algebraic = np.eye(n)
+    for i in range(n):
+        for j in range(i + 1, n):
+            sisj_algebraic[i][j] = get_sisj_algebraic(i, j, psi)
+            sisj_algebraic[j][i] = sisj_algebraic[i][j]
+    assert np.allclose(sisj, sisj_algebraic)