diff --git a/src/scripts/2x2_analysis.py b/src/scripts/2x2_analysis.py
deleted file mode 100644
index 07d26a6..0000000
--- a/src/scripts/2x2_analysis.py
+++ /dev/null
@@ -1,48 +0,0 @@
-import matplotlib.pyplot as plt
-
-
-def plot_from_file():
-    figure1, ax1 = plt.subplots()
-    figure2, ax2 = plt.subplots()
-    figure3, ax3 = plt.subplots()
-    figure4, ax4 = plt.subplots()
-
-    label = "something"
-    with open("output/test_2x2.txt") as f:
-        lines = f.readlines()
-
-        t = []
-        energy = []
-        magnetization = []
-        CV = []
-        X = []
-        for line in lines:
-            items = line.strip().split(",")
-            t.append(int(items[0]))
-            energy.append(float(items[1]))
-            magnetization.append(float(items[2]))
-            CV.append(float(items[3]))
-            X.append(float(items[4]))
-
-        ax1.plot(t, energy, label=f"<epsilon> {label}")
-        ax2.plot(t, magnetization, label=f"<|m|> {label}")
-        ax3.plot(t, CV, label=f"CV {label}")
-        ax4.plot(t, X, label=f"X {label}")
-
-        # ax1.set_yscale("log")
-        # ax2.set_yscale("log")
-
-    figure1.legend()
-    figure1.savefig("../latex/images/2x2_energy.pdf")
-    figure2.legend()
-    figure2.savefig("../latex/images/2x2_magnetization.pdf")
-    figure3.legend()
-    figure3.savefig("../latex/images/2x2_heat_capacity.pdf")
-    figure4.legend()
-    figure4.savefig("../latex/images/2x2_susceptibility.pdf")
-
-def main():
-    plot_from_file()
-
-if __name__ == "__main__":
-    main()
diff --git a/src/scripts/burn_in_time.py b/src/scripts/burn_in_time.py
index cd86a97..517ee01 100644
--- a/src/scripts/burn_in_time.py
+++ b/src/scripts/burn_in_time.py
@@ -2,12 +2,16 @@ import matplotlib.pyplot as plt
 
 def plot_from_file():
     files = [
-        "output/burn_in_time_1_0.txt",
-        "output/burn_in_time_2_4.txt",
+        "output/burn_in_time/unordered_1_0.txt",
+        "output/burn_in_time/ordered_1_0.txt",
+        "output/burn_in_time/unordered_2_4.txt",
+        "output/burn_in_time/ordered_2_4.txt",
     ]
     labels = [
-        "1.0",
-        "2.4"
+        "1.0, unordered",
+        "1.0, ordered",
+        "2.4, unordered",
+        "2.4, ordered"
     ]
     figure1, ax1 = plt.subplots()
     figure2, ax2 = plt.subplots()
@@ -25,8 +29,8 @@ def plot_from_file():
                 energy.append(float(items[1]))
                 magnetization.append(float(items[5]))
 
-            ax1.plot(t, energy, label=f"<epsilon> {label}")
-            ax2.plot(t, magnetization, label=f"<|m|> {label}")
+            ax1.plot(t, energy, label=fr"$\langle \epsilon \rangle$ {label}")
+            ax2.plot(t, magnetization, label=fr"$\langle | m | \rangle$ {label}")
 
     figure1.legend()
     figure1.savefig("../latex/images/burn_in_time_energy.pdf")
diff --git a/src/scripts/pd_estimate.py b/src/scripts/pd_estimate.py
index c5cd327..9ee982a 100644
--- a/src/scripts/pd_estimate.py
+++ b/src/scripts/pd_estimate.py
@@ -1,8 +1,13 @@
-import matplotlib.pyplot as plt
-import numpy as np
+import os
 from pathlib import Path
 
+import matplotlib.pyplot as plt
+import numpy as np
+
+
 def plot(infile, outfile):
+    if not (outdir := Path(outfile).parent).exists():
+        os.makedirs(outdir)
 
     figure1, ax1 = plt.subplots()
     arr = []
@@ -12,13 +17,11 @@ def plot(infile, outfile):
             vals = line.strip().split(",")
             arr.append(float(vals[0]))
 
-    ax1.hist(arr, np.arange(min(arr), max(arr) + .02, .02), density=True)
+    ax1.hist(arr, np.arange(min(arr), max(arr) + 0.02, 0.02), density=True, ec="black")
 
     figure1.savefig(outfile)
 
+
 if __name__ == "__main__":
     plot("output/pd_estimate/estimate_1_0.txt", "../latex/images/pd_estimate_1_0.pdf")
     plot("output/pd_estimate/estimate_2_4.txt", "../latex/images/pd_estimate_2_4.pdf")
-
-
-
diff --git a/src/scripts/phase_transition.py b/src/scripts/phase_transition.py
index 7782041..2be3b98 100644
--- a/src/scripts/phase_transition.py
+++ b/src/scripts/phase_transition.py
@@ -1,7 +1,10 @@
-import matplotlib.pyplot as plt
-
 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",
@@ -22,6 +25,12 @@ def plot_phase_transition(indir, outdir):
     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 = []
@@ -30,6 +39,10 @@ def plot_phase_transition(indir, outdir):
         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:
@@ -40,22 +53,60 @@ def plot_phase_transition(indir, outdir):
                 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)
+        # 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/10M/", "../latex/images/phase_transition/wide/10M/")
-    plot_phase_transition("fox_output/phase_transition/1M/", "../latex/images/phase_transition/wide/1M/")
-    plot_phase_transition("fox_output/phase_transition/narrow/10M/", "../latex/images/phase_transition/narrow/10M/")
+    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/",
+    )