diff --git a/src/scripts/burn_in_time.py b/python_scripts/burn_in_time.py
similarity index 100%
rename from src/scripts/burn_in_time.py
rename to python_scripts/burn_in_time.py
diff --git a/src/scripts/pd_estimate.py b/python_scripts/pd_estimate.py
similarity index 100%
rename from src/scripts/pd_estimate.py
rename to python_scripts/pd_estimate.py
diff --git a/python_scripts/phase_transition.py b/python_scripts/phase_transition.py
new file mode 100644
index 0000000..7751254
--- /dev/null
+++ b/python_scripts/phase_transition.py
@@ -0,0 +1,191 @@
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy.stats import linregress
+
+
+def plot_phase_transition_alt(indir, outdir):
+    files = [
+        "size_20.txt",
+        "size_40.txt",
+        "size_60.txt",
+        "size_80.txt",
+        "size_100.txt",
+        "size_500.txt",
+    ]
+    labels = ["L = 20", "L = 40", "L = 60", "L = 80", "L = 100", "L = 500"]
+
+    figure1, ax1 = plt.subplots()
+    figure2, ax2 = plt.subplots()
+    figure3, ax3 = plt.subplots()
+    figure4, ax4 = plt.subplots()
+    figure5, ax5 = plt.subplots()
+
+    # For linear regression
+    L = []
+    Tc = []
+    size = 20
+
+    for file, label in zip(files, labels):
+        t = []
+        e = []
+        m = []
+        CV = []
+        X = []
+
+        # Append the lattice size
+        L.append(size)
+        size += 20
+
+        with open(Path(indir, file)) as f:
+            lines = f.readlines()
+            for line in lines:
+                l = line.strip().split(",")
+                t.append(float(l[0]))
+                e.append(float(l[1]))
+                m.append(float(l[2]))
+                CV.append(float(l[3]))
+                X.append(float(l[4]))
+
+        # Append the critical temp for the current lattice size
+        Tc.append(t[X.index(max(X))])
+
+        ax1.plot(t, e, label=label)
+        ax2.plot(t, m, label=label)
+        ax3.plot(t, CV, label=label)
+        ax4.plot(t, X, label=label)
+
+    inv_L = list(map(lambda x: 1 / x, L))
+    # Attempt linear regression
+    x = np.linspace(0, 1 / 20, 1001)
+    regression = linregress(inv_L, Tc)
+    f = lambda x: regression[0] * x + regression[1]
+    ax5.scatter(inv_L, Tc)
+    ax5.plot(x, f(x), label=f"m = {regression[0]}, i = {regression[1]}")
+
+    figure1.legend()
+    figure2.legend()
+    figure3.legend()
+    figure4.legend()
+    figure5.legend()
+
+    figure1.savefig(Path(outdir, "energy.pdf"))
+    figure2.savefig(Path(outdir, "magnetization.pdf"))
+    figure3.savefig(Path(outdir, "heat_capacity.pdf"))
+    figure4.savefig(Path(outdir, "susceptibility.pdf"))
+    figure5.savefig(Path(outdir, "linreg.pdf"))
+
+    plt.close(figure1)
+    plt.close(figure2)
+    plt.close(figure3)
+    plt.close(figure4)
+    plt.close(figure5)
+
+
+def plot_phase_transition(indir, outdir):
+    files = [
+        "size_20.txt",
+        "size_40.txt",
+        "size_60.txt",
+        "size_80.txt",
+        "size_100.txt",
+    ]
+    labels = [
+        "L = 20",
+        "L = 40",
+        "L = 60",
+        "L = 80",
+        "L = 100",
+    ]
+
+    figure1, ax1 = plt.subplots()
+    figure2, ax2 = plt.subplots()
+    figure3, ax3 = plt.subplots()
+    figure4, ax4 = plt.subplots()
+    figure5, ax5 = plt.subplots()
+
+    # For linear regression
+    L = []
+    Tc = []
+    size = 20
+
+    for file, label in zip(files, labels):
+        t = []
+        e = []
+        m = []
+        CV = []
+        X = []
+
+        # Append the lattice size
+        L.append(size)
+        size += 20
+
+        with open(Path(indir, file)) as f:
+            lines = f.readlines()
+            for line in lines:
+                l = line.strip().split(",")
+                t.append(float(l[0]))
+                e.append(float(l[1]))
+                m.append(float(l[2]))
+                CV.append(float(l[3]))
+                X.append(float(l[4]))
+
+        # Append the critical temp for the current lattice size
+        Tc.append(t[X.index(max(X))])
+
+        ax1.plot(t, e, label=label)
+        ax2.plot(t, m, label=label)
+        ax3.plot(t, CV, label=label)
+        ax4.plot(t, X, label=label)
+
+    inv_L = list(map(lambda x: 1 / x, L))
+    # Attempt linear regression
+    x = np.linspace(0, 1 / 20, 1001)
+    regression = linregress(inv_L, Tc)
+    f = lambda x: regression[0] * x + regression[1]
+    ax5.scatter(inv_L, Tc)
+    ax5.plot(x, f(x), label=f"m = {regression[0]}, i = {regression[1]}")
+
+    figure1.legend()
+    figure2.legend()
+    figure3.legend()
+    figure4.legend()
+    figure5.legend()
+
+    figure1.savefig(Path(outdir, "energy.pdf"))
+    figure2.savefig(Path(outdir, "magnetization.pdf"))
+    figure3.savefig(Path(outdir, "heat_capacity.pdf"))
+    figure4.savefig(Path(outdir, "susceptibility.pdf"))
+    figure5.savefig(Path(outdir, "linreg.pdf"))
+
+    plt.close(figure1)
+    plt.close(figure2)
+    plt.close(figure3)
+    plt.close(figure4)
+    plt.close(figure5)
+
+
+if __name__ == "__main__":
+    plot_phase_transition_alt(
+        "fox_output/phase_transition/wide/10M/",
+        "latex/images/phase_transition/fox/wide/10M/",
+    )
+    plot_phase_transition(
+        "fox_output/phase_transition/wide/1M/",
+        "latex/images/phase_transition/fox/wide/1M/",
+    )
+    plot_phase_transition(
+        "fox_output/phase_transition/narrow/10M/",
+        "latex/images/phase_transition/fox/narrow/10M/",
+    )
+    plot_phase_transition(
+        "output/phase_transition/", "latex/images/phase_transition/hp/"
+    )
+    plot_phase_transition(
+        "output/phase_transition/", "latex/images/phase_transition/hp/"
+    )
+    plot_phase_transition(
+        "output/phase_transition/",
+        "latex/images/phase_transition/hp/",
+    )
diff --git a/python_scripts/test_burn_in.py b/python_scripts/test_burn_in.py
new file mode 100644
index 0000000..afb6728
--- /dev/null
+++ b/python_scripts/test_burn_in.py
@@ -0,0 +1,68 @@
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy.stats import linregress
+
+
+def plot_phase_transition(indir, outdir):
+    files = [
+        "no_burn_in.txt",
+        "burn_in.txt",
+    ]
+    labels = [
+        "Without burn-in time",
+        "With burn-in time",
+    ]
+
+    figure1, ax1 = plt.subplots()
+    figure2, ax2 = plt.subplots()
+    figure3, ax3 = plt.subplots()
+    figure4, ax4 = plt.subplots()
+
+
+    for file, label in zip(files, labels):
+        t = []
+        e = []
+        m = []
+        CV = []
+        X = []
+
+        with open(Path(indir, file)) as f:
+            lines = f.readlines()
+            for line in lines:
+                l = line.strip().split(",")
+                t.append(float(l[0]))
+                e.append(float(l[1]))
+                m.append(float(l[2]))
+                CV.append(float(l[3]))
+                X.append(float(l[4]))
+
+
+        ax1.plot(t, e, label=label)
+        ax2.plot(t, m, label=label)
+        ax3.plot(t, CV, label=label)
+        ax4.plot(t, X, label=label)
+
+
+    figure1.legend()
+    figure2.legend()
+    figure3.legend()
+    figure4.legend()
+
+    figure1.savefig(Path(outdir, "energy.pdf"))
+    figure2.savefig(Path(outdir, "magnetization.pdf"))
+    figure3.savefig(Path(outdir, "heat_capacity.pdf"))
+    figure4.savefig(Path(outdir, "susceptibility.pdf"))
+
+    plt.close(figure1)
+    plt.close(figure2)
+    plt.close(figure3)
+    plt.close(figure4)
+
+
+if __name__ == "__main__":
+    plot_phase_transition(
+        "output/test_burn_in_time/",
+        "../latex/images/test_burn_in",
+    )
diff --git a/python_scripts/timing.py b/python_scripts/timing.py
new file mode 100644
index 0000000..61d17d4
--- /dev/null
+++ b/python_scripts/timing.py
@@ -0,0 +1,54 @@
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy.stats import linregress
+
+
+def plot_timing(indir, outdir):
+    files = [
+        "lattice_sizes.txt",
+        "sample_sizes.txt",
+    ]
+    labels = [
+        "Lattice sizes",
+        "Sample sizes",
+    ]
+    xlabels = [
+        "Lattice size",
+        "Sampling size"
+    ]
+    outfiles = [
+        "lattice_size.pdf",
+        "sample_sizes.pdf"
+    ]
+
+    for file, label, xlabel, outfile in zip(files, labels, xlabels, outfiles):
+        figure1, ax1 = plt.subplots()
+        x = []
+        t = []
+
+        with open(Path(indir, file)) as f:
+            lines = f.readlines()
+            for line in lines:
+                l = line.strip().split(",")
+                x.append(float(l[0]))
+                t.append(float(l[1]))
+
+
+        ax1.plot(x, t, label=label)
+        ax1.set_xlabel(xlabel)
+        ax1.set_ylabel("time (seconds)")
+
+        figure1.legend()
+
+        figure1.savefig(Path(outdir, outfile))
+
+        plt.close(figure1)
+
+
+if __name__ == "__main__":
+    plot_timing(
+        "output/timing/",
+        "../latex/images/timing",
+    )
diff --git a/src/scripts/phase_transition.py b/src/scripts/phase_transition.py
deleted file mode 100644
index 2be3b98..0000000
--- a/src/scripts/phase_transition.py
+++ /dev/null
@@ -1,112 +0,0 @@
-from pathlib import Path
-
-import matplotlib.pyplot as plt
-import numpy as np
-from scipy.stats import linregress
-
-
-def plot_phase_transition(indir, outdir):
-    files = [
-        "size_20.txt",
-        "size_40.txt",
-        "size_60.txt",
-        "size_80.txt",
-        "size_100.txt",
-    ]
-    labels = [
-        "L = 20",
-        "L = 40",
-        "L = 60",
-        "L = 80",
-        "L = 100",
-    ]
-
-    figure1, ax1 = plt.subplots()
-    figure2, ax2 = plt.subplots()
-    figure3, ax3 = plt.subplots()
-    figure4, ax4 = plt.subplots()
-    figure5, ax5 = plt.subplots()
-
-    # For linear regression
-    L = []
-    Tc = []
-    size = 20
-
-    for file, label in zip(files, labels):
-        t = []
-        e = []
-        m = []
-        CV = []
-        X = []
-
-        # Append the lattice size
-        L.append(size)
-        size += 20
-
-        with open(Path(indir, file)) as f:
-            lines = f.readlines()
-            for line in lines:
-                l = line.strip().split(",")
-                t.append(float(l[0]))
-                e.append(float(l[1]))
-                m.append(float(l[2]))
-                CV.append(float(l[3]))
-                X.append(float(l[4]))
-
-        # Append the critical temp for the current lattice size
-        Tc.append(t[X.index(max(X))])
-
-        ax1.plot(t, e, label=label)
-        ax2.plot(t, m, label=label)
-        ax3.plot(t, CV, label=label)
-        ax4.plot(t, X, label=label)
-
-    # Attempt linear regression
-    x = np.linspace(0, 100, 1001)
-    regression = linregress(L, Tc)
-    f = lambda x: regression[0] * x + regression[1]
-    ax5.scatter(L, Tc)
-    ax5.plot(x, f(x), label=f"m = {regression[0]}")
-
-    figure1.legend()
-    figure2.legend()
-    figure3.legend()
-    figure4.legend()
-    figure5.legend()
-
-    figure1.savefig(Path(outdir, "energy.pdf"))
-    figure2.savefig(Path(outdir, "magnetization.pdf"))
-    figure3.savefig(Path(outdir, "heat_capacity.pdf"))
-    figure4.savefig(Path(outdir, "susceptibility.pdf"))
-    figure5.savefig(Path(outdir, "linreg.pdf"))
-
-    plt.close(figure1)
-    plt.close(figure2)
-    plt.close(figure3)
-    plt.close(figure4)
-    plt.close(figure5)
-
-
-if __name__ == "__main__":
-    plot_phase_transition(
-        "fox_output/phase_transition/wide/10M/",
-        "../latex/images/phase_transition/fox/wide/10M/",
-    )
-    plot_phase_transition(
-        "fox_output/phase_transition/wide/1M/",
-        "../latex/images/phase_transition/fox/wide/1M/",
-    )
-    plot_phase_transition(
-        "fox_output/phase_transition/narrow/10M/",
-        "../latex/images/phase_transition/fox/narrow/10M/",
-    )
-    plot_phase_transition(
-        "output/phase_transition/", "../latex/images/phase_transition/hp/"
-    )
-    plot_phase_transition(
-        "output/phase_transition/", "../latex/images/phase_transition/hp/"
-    )
-    plot_phase_transition(
-        "output/phase_transition/",
-        "../latex/images/phase_transition/hp/",
-    )