diff --git a/exercises/surface17_depolarizing.py b/exercises/surface17_depolarizing.py
new file mode 100644
index 0000000000000000000000000000000000000000..1ce7366f2cc36ced6320c9ded416b2e804ff1974
--- /dev/null
+++ b/exercises/surface17_depolarizing.py
@@ -0,0 +1,187 @@
+#!/usr/bin/env python3
+'''
+Solution for exercise 13.1f, Monte Carlo simulation of surface-17
+quantum error correction for depolarizing noise.
+This script is much more efficient if the numba package is available.
+'''
+import numpy as np
+import matplotlib.pyplot as plt
+
+# logical X and Z in symplectic notation
+xL = 0b001010100
+zL = 0b100010001
+# X stabilizer generators in symplectic notation
+sx1 = 0b000011011
+sx2 = 0b000100100
+sx3 = 0b001001000
+sx4 = 0b110110000
+# Z stabilizer generators in symplectic notation
+sz1 = 0b000000011
+sz2 = 0b000110110
+sz3 = 0b011011000
+sz4 = 0b110000000
+xstabilizerGenerators = (sx1,sx2,sx3,sx4)
+zstabilizerGenerators = (sz1,sz2,sz3,sz4)
+
+# Look-up table mapping 9 bit integers to the parity of the number of bits:
+# Definition: parity[x] = 1 if bitstring(x) has odd number of '1's, otherwise 0
+parity = bytes(bin(x).count('1') % 2 for x in range(2**9))
+
+# Z corrections detected by X stabilizers
+zcorrections = (
+        # binary data show which data qubit should be corrected:
+        # 987654321 ← data qubit number
+        0b000000000, # syndrom: -,       correction: -
+        0b000000001, # syndrom: 1,       correction: 1 (could also be 2)
+        0b000000100, # syndrom: 2,       correction: 3
+        0b000000101, # syndrom: 1,2,     correction: 1,3 (could also be 2,3 or 5,6)
+        0b001000000, # syndrom: 3,       correction: 7
+        0b000001000, # syndrom: 1,3,     correction: 4
+        0b001000100, # syndrom: 2,3,     correction: 3,7
+        0b000001100, # syndrom: 1,2,3,   correction: 3,4
+        0b100000000, # syndrom: 4,       correction: 9 (could also be 8)
+        0b000010000, # syndrom: 1,4,     correction: 5
+        0b000100000, # syndrom: 2,4,     correction: 6
+        0b000010100, # syndrom: 1,2,4,   correction: 3,5 (could also be 1,6)
+        0b101000000, # syndrom: 3,4,     correction: 7,9 (could also be 7,8 or 4,5)
+        0b001010000, # syndrom: 1,3,4,   correction: 5,7 (could also be 4,9 or 4,8)
+        0b001100000, # syndrom: 2,3,4,   correction: 6,7
+        0b000101000, # syndrom: 1,2,3,4, correction: 4,6
+        )
+
+# X corrections detected by Z stabilizers
+xcorrections = (
+        # binary data show which data qubit should be corrected:
+        # 987654321 ← data qubit number
+        0b000000000, # syndrom: -,       correction: -
+        0b000000001, # syndrom: 1,       correction: 1
+        0b000000100, # syndrom: 2,       correction: 3 (could also be 6)
+        0b000000010, # syndrom: 1,2,     correction: 2
+        0b001000000, # syndrom: 3,       correction: 7 (could also be 4)
+        0b001000001, # syndrom: 1,3,     correction: 1,7 (could also be 1,4 or 2,5)
+        0b000010000, # syndrom: 2,3,     correction: 5
+        0b000010001, # syndrom: 1,2,3,   correction: 1,5 (could also be 2,4 or 2,7)
+        0b100000000, # syndrom: 4,       correction: 9
+        0b100000001, # syndrom: 1,4,     correction: 1,9
+        0b100000100, # syndrom: 2,4,     correction: 3,9 (could also be 6,9 or 5,8)
+        0b100000010, # syndrom: 1,2,4,   correction: 2,9
+        0b010000000, # syndrom: 3,4,     correction: 8
+        0b010000001, # syndrom: 1,3,4,   correction: 1,8
+        0b100010000, # syndrom: 2,3,4,   correction: 5,9 (could also be 3,8 or 6,8)
+        0b010000010, # syndrom: 1,2,3,4, correction: 2,8
+        )
+
+
+################################################################
+# THIS BLOCK COULD BE OMITTED, but it can significantly improve the performance
+# if the numba package is available. Numba provides just-in-time compilation
+# and parallelization that speed up the for loop in montecarlo().
+try:
+    # Try to import just-in-time compilation function decorator from numba.
+    from numba import jit, prange
+    NUMBA = True
+except ImportError:
+    # prange behaves like range, but allows for parallelization. Replace it
+    # with range if numba is not available.
+    prange = range
+    # If numba is not available, jit should take a dummy value.
+    jit = lambda *args, **kwargs: (lambda x: x)
+    NUMBA = False
+# The decorator @jit causes this function to be compiled and not just
+# parsed by the python interpreter:
+@jit(nopython=True, parallel=True, cache=False)
+################################################################
+def montecarlo(p, samples):
+    '''
+    Simulate surface 17 code with error rate p/3 for X, Y and Z errors on data
+    qubits for given number of samples. Return number of logical X, Y, and Z
+    errors.
+    '''
+    xfailures = 0
+    yfailures = 0
+    zfailures = 0
+    # sample a reasonable amount of times given the probability of errors
+    for _ in prange(samples):
+
+        # start out with clean state and randomly place X, Z and XZ errors on data qubits
+        xerrors = 0
+        zerrors = 0
+        # Iterate over the 9 qubits
+        for i in range(9):
+            r = np.random.random()
+            if r < p/3:
+                xerrors |= 2**i
+            elif r < 2*p/3:
+                zerrors |= 2**i
+            elif r < p:
+                xerrors |= 2**i
+                zerrors |= 2**i
+
+        # read syndrome and apply correction operators according to look up table
+        xsyndrome = 0
+        zsyndrome = 0
+        for i in range(len(xstabilizerGenerators)):
+            # for all stabilizers: check whether the stabilizer commutes or
+            # anticommutes with the error
+            xsyndrome |= parity[xstabilizerGenerators[i] & zerrors] * 2**i
+            zsyndrome |= parity[zstabilizerGenerators[i] & xerrors] * 2**i
+        zcorrection = zcorrections[xsyndrome]
+        xcorrection = xcorrections[zsyndrome]
+        x_after_correction = xerrors ^ xcorrection
+        z_after_correction = zerrors ^ zcorrection
+
+        # check if the resulting configuration has caused a logical Z error
+        # logical Z should always anticommute with logical X
+        x_logical_error = parity[x_after_correction & zL]
+        z_logical_error = parity[z_after_correction & xL]
+        if x_logical_error and z_logical_error:
+            yfailures += 1
+        elif x_logical_error:
+            xfailures += 1
+        elif z_logical_error:
+            zfailures += 1
+    return xfailures, yfailures, zfailures
+
+
+def main(relative_samples=1000, resolution=16):
+    '''
+    Simulate and plot logical error rates of surface 17 code with depolarizing
+    noise.
+    '''
+    # set range of physical error rates to MC sample
+    pRange = np.logspace(-3, 0, resolution)
+    # define number of samples for each physical error rate
+    sampleRange = np.int_(relative_samples/pRange)
+
+    x_failure_rates = np.empty_like(pRange)
+    y_failure_rates = np.empty_like(pRange)
+    z_failure_rates = np.empty_like(pRange)
+    # run the MC simulation
+    for i in range(resolution):
+        xfailures, yfailures, zfailures = montecarlo(pRange[i], sampleRange[i])
+        x_failure_rates[i] = xfailures / sampleRange[i]
+        y_failure_rates[i] = yfailures / sampleRange[i]
+        z_failure_rates[i] = zfailures / sampleRange[i]
+
+    total_failure_rates = x_failure_rates + y_failure_rates + z_failure_rates
+
+    # Plot the results
+    fig, ax = plt.subplots()
+    ax.grid()
+    ax.set_xlabel('$p$')
+    ax.set_ylabel('$p_L$')
+    ax.set_title('Pseudothreshold of Surface-17, depolarizing noise')
+    ax.set_xscale('log')
+    ax.set_yscale('log')
+    ax.plot(pRange, x_failure_rates, '-x', label='logical X error')
+    ax.plot(pRange, y_failure_rates, '-x', label='logical Y error')
+    ax.plot(pRange, z_failure_rates, '-x', label='logical Z error')
+    ax.plot(pRange, total_failure_rates, '-x', label='any logical error')
+    ax.plot(pRange, pRange, 'k:', alpha=0.5)
+    ax.legend()
+    fig.tight_layout()
+    plt.show()
+
+
+if __name__ == '__main__':
+    main(10000 if NUMBA else 1000)
diff --git a/exercises/surface17_depolarizing_optimized.py b/exercises/surface17_depolarizing_optimized.py
new file mode 100644
index 0000000000000000000000000000000000000000..5be2ebd09b2922c69b6b91e16ed97cceaa2a4762
--- /dev/null
+++ b/exercises/surface17_depolarizing_optimized.py
@@ -0,0 +1,162 @@
+#!/usr/bin/env python3
+'''
+Solution for exercise 13.1f, Monte Carlo simulation of surface-17
+quantum error correction for depolarizing noise with some optimizations.
+This script is much more efficient if the numba package is available.
+'''
+import numpy as np
+from random import random
+import matplotlib.pyplot as plt
+
+# logical X and Z in symplectic notation
+xL = 0b001010100
+zL = 0b100010001
+# X stabilizer generators in symplectic notation
+sx1 = 0b000011011
+sx2 = 0b000100100
+sx3 = 0b001001000
+sx4 = 0b110110000
+# Z stabilizer generators in symplectic notation
+sz1 = 0b000000011
+sz2 = 0b000110110
+sz3 = 0b011011000
+sz4 = 0b110000000
+xstabilizerGenerators = (sx1,sx2,sx3,sx4)
+zstabilizerGenerators = (sz1,sz2,sz3,sz4)
+
+parity = bytes(bin(x).count('1') % 2 for x in range(1<<9))
+
+# Z corrections detected by X stabilizers
+zcorrections = (
+        0b000000000, # syndrom: -,       correction: -
+        0b000000001, # syndrom: 1,       correction: 1 (could also be 2)
+        0b000000100, # syndrom: 2,       correction: 3
+        0b000000101, # syndrom: 1,2,     correction: 1,3 (could also be 2,3 or 5,6)
+        0b001000000, # syndrom: 3,       correction: 7
+        0b000001000, # syndrom: 1,3,     correction: 4
+        0b001000100, # syndrom: 2,3,     correction: 3,7
+        0b000001100, # syndrom: 1,2,3,   correction: 3,4
+        0b100000000, # syndrom: 4,       correction: 9 (could also be 8)
+        0b000010000, # syndrom: 1,4,     correction: 5
+        0b000100000, # syndrom: 2,4,     correction: 6
+        0b000010100, # syndrom: 1,2,4,   correction: 3,5 (could also be 1,6)
+        0b101000000, # syndrom: 3,4,     correction: 7,9 (could also be 7,8 or 4,5)
+        0b001010000, # syndrom: 1,3,4,   correction: 5,7 (could also be 4,9 or 4,8)
+        0b001100000, # syndrom: 2,3,4,   correction: 6,7
+        0b000101000, # syndrom: 1,2,3,4, correction: 4,6
+        )
+
+# X corrections detected by Z stabilizers
+xcorrections = (
+        0b000000000, # syndrom: -,       correction: -
+        0b000000001, # syndrom: 1,       correction: 1
+        0b000000100, # syndrom: 2,       correction: 3 (could also be 6)
+        0b000000010, # syndrom: 1,2,     correction: 2
+        0b001000000, # syndrom: 3,       correction: 7 (could also be 4)
+        0b001000001, # syndrom: 1,3,     correction: 1,7 (could also be 1,4 or 2,5)
+        0b000010000, # syndrom: 2,3,     correction: 5
+        0b000010001, # syndrom: 1,2,3,   correction: 1,5 (could also be 2,4 or 2,7)
+        0b100000000, # syndrom: 4,       correction: 9
+        0b100000001, # syndrom: 1,4,     correction: 1,9
+        0b100000100, # syndrom: 2,4,     correction: 3,9 (could also be 6,9 or 5,8)
+        0b100000010, # syndrom: 1,2,4,   correction: 2,9
+        0b010000000, # syndrom: 3,4,     correction: 8
+        0b010000001, # syndrom: 1,3,4,   correction: 1,8
+        0b100010000, # syndrom: 2,3,4,   correction: 5,9 (could also be 3,8 or 6,8)
+        0b010000010, # syndrom: 1,2,3,4, correction: 2,8
+        )
+
+direct_xcorrections = bytes(parity[(c ^ xcorrections[sum(parity[s&c]<<i for i,s in enumerate(zstabilizerGenerators))]) & zL] for c in range(1<<9))
+direct_zcorrections = bytes(parity[(c ^ zcorrections[sum(parity[s&c]<<i for i,s in enumerate(xstabilizerGenerators))]) & xL] for c in range(1<<9))
+
+################################################################
+# THIS BLOCK COULD BE OMITTED, but it can significantly improve the performance
+# if the numba package is available. Numba provides just-in-time compilation
+# and parallelization that speed up the for loop in montecarlo().
+try:
+    # Try to import just-in-time compilation function decorator from numba.
+    from numba import jit, prange
+except ImportError:
+    # prange behaves like range, but allows for parallelization. Replace it
+    # with range if numba is not available.
+    prange = range
+    # If numba is not available, jit should take a dummy value.
+    jit = lambda *args, **kwargs: (lambda x: x)
+# The decorator @jit causes this function to be compiled and not just
+# parsed by the python interpreter:
+@jit(nopython=True, parallel=True, cache=False)
+################################################################
+def montecarlo(p, samples):
+    '''
+    Simulate surface 17 code with error rate p/3 for X, Y and Z errors on data
+    qubits for given number of samples. Return number of logical X, Y, and Z
+    errors.
+    '''
+    xfailures = 0
+    yfailures = 0
+    zfailures = 0
+    # sample a reasonable amount of times given the probability of errors
+    for _ in prange(samples):
+
+        # start out with clean state and randomly place Z-errors on data qubits
+        xerrors = 0
+        zerrors = 0
+        for i in range(9):
+            r = random()
+            #if r < p:
+            #    zerrors |= 1 << i
+            if r < p/3:
+                xerrors |= 1 << i
+            elif r < 2*p/3:
+                zerrors |= 1 << i
+            elif r < p:
+                xerrors |= 1 << i
+                zerrors |= 1 << i
+
+        # check if the resulting configuration has caused a logical Z error
+        # logical Z should always anticommute with logical X
+        x_logical_error = direct_xcorrections[xerrors]
+        z_logical_error = direct_zcorrections[zerrors]
+        if x_logical_error and z_logical_error:
+            yfailures += 1
+        elif x_logical_error:
+            xfailures += 1
+        elif z_logical_error:
+            zfailures += 1
+    return xfailures, yfailures, zfailures
+
+
+def main(relative_samples=10000, resolution=16):
+    '''
+    Simulate and plot logical error rates of surface 17 code with depolarizing
+    noise.
+    '''
+    # set range of physical error rates to MC sample
+    pRange = np.logspace(-3, 0, resolution)
+    # define number of samples for each physical error rate
+    sampleRange = np.int_(relative_samples/pRange)
+
+    # run the MC simulation
+    logical_error_counts = np.array([montecarlo(*ps) for ps in zip(pRange, sampleRange)])
+    # calculate logical failure rate from samples
+    logical_error_rates = logical_error_counts / sampleRange.reshape((resolution,1))
+
+    fig, ax = plt.subplots()
+    ax.grid()
+    ax.set_xlabel('$p$')
+    ax.set_ylabel('$p_L$')
+    ax.set_title('Pseudothreshold of Surface-17, depolarizing noise')
+    ax.set_xscale('log')
+    ax.set_yscale('log')
+    ax.plot(pRange, x_failure_rates, '-x', label='logical X error')
+    ax.plot(pRange, y_failure_rates, '-x', label='logical Y error')
+    ax.plot(pRange, z_failure_rates, '-x', label='logical Z error')
+    ax.plot(pRange, total_failure_rates, '-x', label='any logical error')
+    ax.plot(pRange, pRange, 'k:', alpha=0.5)
+    ax.legend()
+    fig.tight_layout()
+    plt.show()
+
+
+if __name__ == '__main__':
+    main()
diff --git a/exercises/surface17_exercise.py b/exercises/surface17_exercise.py
new file mode 100644
index 0000000000000000000000000000000000000000..bbbf8cc58e042f9ba01a84252d8632763b1bdab6
--- /dev/null
+++ b/exercises/surface17_exercise.py
@@ -0,0 +1,101 @@
+#!/usr/bin/env python3
+import numpy as np
+import matplotlib.pyplot as plt
+
+# logical X in symplectic notation
+xL = 0b001010100
+# X stabilizer generators in symplectic notation
+s1 = 0b000011011
+s2 = 0b000100100
+s3 = 0b001001000
+s4 = 0b110110000
+stabilizerGenerators = [s1,s2,s3,s4]
+
+def getCommutator(a,b):
+	# For two multiqubit Pauli X- or Z- operators we want to know
+	# whether they commute (0) or anticommute (1). Because single qubit Paulis
+	# anticommute we just sum up (mod 2) the mutual ones in the bitstings 
+	# that represent the multiqubit operators in symplectic notation
+	return bin(a & b)[2:].count('1')%2
+
+def getSyndrome(configuration):
+	# commutation/anticommutation with a stabilizer is marked as 0/1 in syndrome
+	syndrome = 0
+	for i in range(len(stabilizerGenerators)):
+		syndrome += getCommutator(stabilizerGenerators[i], configuration)*2**i
+	return syndrome
+
+def getCorrection(syndrome):
+	# Look Up Table: apply minimal weight correction based on syndrome information
+	correction = []
+	bin_correction = 0
+	if syndrome == 0b0001: correction = [0]
+	if syndrome == 0b1110: correction = [5,6]
+	if syndrome == 0b0010: correction = [2]
+	if syndrome == 0b1000: correction = [8]
+	if syndrome == 0b0100: correction = [6]
+	# TASK: add the correct syndromes and corrections
+	if syndrome == ...: correction = [3]
+	if syndrome == ...: correction = [2,4]
+	if syndrome == 0b1101: correction = ...
+	if syndrome == 0b1001: correction = ...
+	if syndrome == 0b1010: correction = [5]
+	if syndrome == 0b1100: correction = [6,8]
+	if syndrome == 0b0111: correction = [5,6]
+	if syndrome == 0b0011: correction = [0,2]
+	if syndrome == 0b0110: correction = [2,6]
+	if syndrome == 0b1111: correction = [3,5]
+	for i in correction:
+		bin_correction += 2**i
+	return bin_correction
+
+# set range of physical error rates to MC sample
+pRange = np.logspace(-3,0,7)
+
+fail_log = []
+logical_failure_rate = []
+for p in range(len(pRange)):
+	
+	fail_log.append([])
+	# sample a reasonable amount of times given the probability of errors
+	for _ in range(int(1000/pRange[p])):
+		
+		# start out with clean state and randomly place Z-errors on data qubits
+		configuration = 0
+		for i in range(9):
+			if np.random.random() < pRange[p]: configuration += 2**i
+		
+		# TASK: add syndrome readout and apply correction operators according to look up table
+		syndrome = ... 
+		correction = ... 
+		after_correction = configuration^correction
+
+		# TASK: check if the resulting configuration has caused a logical Z error
+		if ... :
+			fail_log[-1].append(0)
+		else:
+			fail_log[-1].append(1)
+		
+		#print('configuration', format(configuration,'09b'))
+		#print('syndrome', format(syndrome,'04b'))
+		#print('correction', format(correction,'09b'))
+		#print('result', format(after_correction,'09b'))
+		#print('logical X', getCommutator(after_correction, xL))
+	
+	# calculate logical failure rate from samples
+	logical_failure_rate.append(np.sum(fail_log[-1])/len(fail_log[-1]))
+
+#print('physical error rates', pRange)
+#print('failure', logical_failure_rate)
+
+#plt.rcParams['text.usetex'] = True
+#plt.rcParams.update({'font.size': 18})
+
+plt.grid()
+plt.xlabel('$p$')
+plt.ylabel('$p_L$')
+plt.title('Pseudothreshold of Surface-17')
+plt.loglog(pRange,logical_failure_rate,'-x')
+plt.loglog(np.logspace(-3,0,4),np.logspace(-3,0,4),'k:',alpha=0.5)
+plt.tight_layout()
+plt.show()
diff --git a/exercises/surface17_solution.py b/exercises/surface17_solution.py
new file mode 100644
index 0000000000000000000000000000000000000000..a007792474494ccb718263e3b27658a2bd27df25
--- /dev/null
+++ b/exercises/surface17_solution.py
@@ -0,0 +1,108 @@
+#!/usr/bin/env python3
+'''
+Solution for exercise 13.1f, Monte Carlo simulation of surface-17 error
+quantum error correction for dephasing noise.
+A more efficient implementation of this can be found in the solution for
+depolarizing noise.
+'''
+import numpy as np
+import matplotlib.pyplot as plt
+
+# logical X in symplectic notation
+xL = 0b001010100
+# X stabilizer generators in symplectic notation
+s1 = 0b000011011
+s2 = 0b000100100
+s3 = 0b001001000
+s4 = 0b110110000
+stabilizerGenerators = [s1,s2,s3,s4]
+
+def getCommutator(a,b):
+	# For two multiqubit Pauli X- or Z- operators we want to know
+	# whether they commute (0) or anticommute (1). Because single qubit Paulis
+	# anticommute we just sum up (mod 2) the mutual ones in the bitstings 
+	# that represent the multiqubit operators in symplectic notation
+	return bin(a & b)[2:].count('1')%2
+
+def getSyndrome(configuration):
+	# commutation/anticommutation with a stabilizer is marked as 0/1 in syndrome
+	syndrome = 0
+	for i in range(len(stabilizerGenerators)):
+		syndrome += getCommutator(stabilizerGenerators[i], configuration)*2**i
+	return syndrome
+
+def getCorrection(syndrome):
+	# Look Up Table: apply minimal weight correction based on syndrome information
+	correction = []
+	bin_correction = 0
+	if syndrome == 0b0001: correction = [0] # could also be [1]
+	if syndrome == 0b1110: correction = [5,6]
+	if syndrome == 0b0010: correction = [2]
+	if syndrome == 0b1000: correction = [8] # could also be [7]
+	if syndrome == 0b0100: correction = [6]
+	if syndrome == 0b0101: correction = [3]
+	if syndrome == 0b1011: correction = [2,4]
+	if syndrome == 0b1101: correction = [4,6] # could also be [3,8]
+	if syndrome == 0b1001: correction = [4]
+	if syndrome == 0b1010: correction = [5]
+	if syndrome == 0b1100: correction = [6,8]
+	if syndrome == 0b0111: correction = [5,6]
+	if syndrome == 0b0011: correction = [0,2]
+	if syndrome == 0b0110: correction = [2,6]
+	if syndrome == 0b1111: correction = [3,5]
+	for i in correction:
+		bin_correction += 2**i
+	return bin_correction
+
+# set range of physical error rates to MC sample
+pRange = np.logspace(-3,0,7)
+
+fail_log = []
+logical_failure_rate = []
+for p in range(len(pRange)):
+	
+	fail_log.append([])
+	# sample a reasonable amount of times given the probability of errors
+	for _ in range(int(1000/pRange[p])):
+		
+		# start out with clean state and randomly place Z-errors on data qubits
+		configuration = 0
+		for i in range(9):
+			if np.random.random() < pRange[p]: configuration += 2**i
+		
+		# read syndrome and apply correction operators according to look up table
+		syndrome = getSyndrome(configuration)
+		correction = getCorrection(syndrome)
+		after_correction = configuration^correction
+
+		# check if the resulting configuration has caused a logical Z error
+		# we notice that we induced a logical Z by applying the correction
+		# by checking if the resulting multiqubit Z-operator commutes with logical X
+		if getCommutator(after_correction, xL) == 0:
+			fail_log[-1].append(0)
+		else:
+			fail_log[-1].append(1)
+		
+		#print('configuration', format(configuration,'09b'))
+		#print('syndrome', format(syndrome,'04b'))
+		#print('correction', format(correction,'09b'))
+		#print('result', format(after_correction,'09b'))
+		#print('logical X', getCommutator(after_correction, xL))
+	
+	# calculate logical failure rate from samples
+	logical_failure_rate.append(np.sum(fail_log[-1])/len(fail_log[-1]))
+
+#print('physical error rates', pRange)
+#print('failure', logical_failure_rate)
+
+#plt.rcParams['text.usetex'] = True
+#plt.rcParams.update({'font.size': 18})
+
+plt.grid()
+plt.xlabel('$p$')
+plt.ylabel('$p_L$')
+plt.title('Pseudothreshold of Surface-17')
+plt.loglog(pRange,logical_failure_rate,'-x')
+plt.loglog(np.logspace(-3,0,4),np.logspace(-3,0,4),'k:',alpha=0.5)
+plt.tight_layout()
+plt.show()