diff --git a/scripts/continuation.py b/scripts/continuation.py
new file mode 100644
index 0000000..75801d6
--- /dev/null
+++ b/scripts/continuation.py
@@ -0,0 +1,230 @@
+import argparse, numpy as np, os
+import matplotlib.pyplot as plt
+import matplotlib.ticker as mtick
+from multiprocessing import Pool, cpu_count
+from pathlib import Path
+
+from squish import Simulation, DomainParams, Energy, ordered, Diagram
+from squish.common import OUTPUT_DIR
+from script_tools import (
+    RC_SETTINGS,
+    get_args,
+    get_data,
+    get_simulation_data,
+    get_ordered_data,
+)
+
+NAME = "Continuation"
+
+
+def proc(combo):
+    dia, min_order, min_disorder, hex_ratios, alphas, sim_alphas, sim_energies, offset, length, num_frames = (
+        combo
+    )
+    sim = dia.sim
+    out_folder = sim.path / "continuation"
+
+    plt.rcParams.update(RC_SETTINGS)
+    plt.rcParams.update({"figure.constrained_layout.use": False})
+    fig = plt.figure(figsize=(30, 15))
+
+    e_hex = ordered.e_hex(sim.domains[0])
+
+    for i in range(length):
+        gs = fig.add_gridspec(1, 2)
+        gs.update(left=0, right=0.98, top=0.93, bottom=0.12, wspace=0.08)
+        ax = fig.add_subplot(gs[1])
+
+        curr_alpha = sim.domains[i].w / sim.domains[i].h
+        curr_vee = 100 * (sim.energies[i] / sim.domains[i].n - e_hex)
+
+        ax.plot(
+            alphas,
+            100 * min_order,
+            color="C1",
+            label="Min Ordered",
+            zorder=10,
+            linestyle="dotted",
+        )
+        ax.plot(
+            alphas,
+            100 * min_disorder,
+            color="C0",
+            label="Min Disordered",
+            zorder=11,
+            linestyle="dotted",
+        )
+        ax.plot(
+            sim_alphas[:175],
+            100 * sim_energies[:175],
+            color="C2",
+            label="_nolegend_",
+            linestyle="dashed",
+        )
+        ax.plot(
+            sim_alphas[174:], 100 * sim_energies[174:], color="C2", label="Continuation"
+        )
+
+        ax.scatter(
+            hex_ratios, [0] * len(hex_ratios), color="C2", s=120, marker="H", zorder=50
+        )
+
+        ax.scatter(
+            curr_alpha,
+            curr_vee,
+            s=250,
+            facecolors="none",
+            edgecolors="C6",
+            linewidth=4,
+            zorder=100,
+        )
+
+        ax.set_xlim(0.3, 1)
+        ax.set_xticks([round(w, 2) for w in alphas[::10]])
+        ax.set_xticklabels([f"{round(w, 3):.2f}" for w in alphas[::10]], rotation=90)
+
+        # start, end = ax.get_ylim()
+        # space = np.linspace(0, end, 20)
+
+        space = np.arange(0, 10, 0.5)
+        ax.set_ylim(-0.15, 9.5)
+        ax.set_yticks(space[::-1])
+        ax.ticklabel_format(axis="y", style="sci")
+
+        ax.set_xlabel("Aspect Ratio")
+        ax.set_ylabel(r"VEE $\left[\times 10^{2}\right]$")
+
+        ax.legend(loc="upper center")
+        # ax.legend()
+        ax.grid(zorder=0)
+
+        props = dict(boxstyle="round", facecolor="white", alpha=0.8, zorder=20)
+        ax.text(
+            0.873,
+            0.97,
+            f"N={sim.domains[i].n}",
+            transform=ax.transAxes,
+            verticalalignment="top",
+            bbox=props,
+        )
+
+        ax1 = fig.add_subplot(gs[0])
+        dia.voronoi_plot(i, ax1)
+        fig.savefig(out_folder / f"img{i+offset:05}.png")
+        fig.clear()
+
+        j = len(list((sim.path / "continuation").iterdir()))
+        hashes = int(21 * j / num_frames)
+        print(
+            f'Generating frames... |{"#"*hashes}{" "*(20-hashes)}|'
+            + f" {j}/{num_frames} frames rendered.",
+            flush=True,
+            end="\r",
+        )
+
+    print(flush=True)
+
+
+def main():
+    parser = argparse.ArgumentParser(
+        description="Makes video of analytic continuation."
+    )
+    parser.add_argument(
+        "sims_path",
+        metavar="sim_dir",
+        help="folder that contains simulation data at various aspect ratios for some N",
+    )
+    parser.add_argument(
+        "shrink_sim",
+        metavar="continuation_sim",
+        help="simulation that contains the continuation data",
+    )
+    args = parser.parse_args()
+    sims_path = Path(args.sims_path)
+
+    data, n, r = get_data(
+        sims_path / "package.pkl", get_simulation_data, args=(sims_path,)
+    )
+    domain, alphas = DomainParams(n, 1, 1, r), data["all"]["alpha"]
+    ordered_data = get_data(
+        OUTPUT_DIR / "OrderedCache" / f"{n}.pkl",
+        get_ordered_data,
+        args=(domain, alphas),
+    )
+
+    min_disorder, max_disorder = [], []
+    for i, energies in enumerate(data["distinct"]["Energy"]):
+        disorder_energies = []
+        for j, energy in enumerate(energies):
+            if not data["distinct"]["Ordered"][i][j]:
+                disorder_energies.append(energy)
+        min_disorder.append(min(disorder_energies))
+        max_disorder.append(max(disorder_energies))
+
+    min_order, min_order_coer = [], []
+    for i, energies in enumerate(ordered_data["Energy"]):
+        min_order.append(energies[0])
+        min_order_coer.append(ordered_data["Coercivity"][i][0])
+    e_hex = ordered.e_hex(domain)
+    min_disorder = np.array(min_disorder) / domain.n - e_hex
+    max_disorder = np.array(max_disorder) / domain.n - e_hex
+    min_order = np.array(min_order) / domain.n - e_hex
+
+    hex_ratios = ordered.hexagon_alpha(domain.n)
+
+    sim = Simulation.from_file(Path(args.shrink_sim))
+    sim_alphas = np.array([x.w / x.h for x in sim.frames])
+    sim_energies = np.array([x.energy for x in sim.frames]) / sim.domain.n - e_hex
+    dia = Diagram(sim, ["voronoi"])
+
+    out_folder = sim.path / "continuation"
+    out_folder.mkdir(exist_ok=True)
+
+    combo_list = []
+    for i in range(cpu_count()):
+        start, end = (
+            int(i * len(sim) / cpu_count()),
+            int((i + 1) * len(sim) / cpu_count()),
+        )
+        new_dia = Diagram(None, ["voronoi"])
+        new_dia.sim = dia.sim.slice(list(range(start, end)))
+        combo_list.append(
+            (
+                new_dia,
+                min_order,
+                min_disorder,
+                hex_ratios,
+                alphas,
+                sim_alphas,
+                sim_energies,
+                start,
+                len(sim.frames[start:end]),
+                len(sim),
+            )
+        )
+
+    print("Starting figure generation...", flush=True)
+    with Pool(cpu_count()) as pool:
+        for _ in pool.imap_unordered(proc, combo_list):
+            pass
+
+    video_path = sim.path / (NAME + ".mp4")
+
+    fps = 30
+    print("Assembling MP4...", flush=True)
+    os.system(
+        f"ffmpeg -hide_banner -loglevel error -r {fps} -i"
+        + f' "{sim.path}/continuation/img%05d.png"'
+        + f" -c:v libx264 -crf 18 -preset slow -pix_fmt yuv420p -vf"
+        + f' "scale=trunc(iw/2)*2:trunc(ih/2)*2" -f mp4 "{video_path}"'
+    )
+
+    print(f"Wrote to {video_path}.")
+
+
+if __name__ == "__main__":
+    os.environ["QT_LOGGING_RULES"] = "*=false"
+    try:
+        main()
+    except KeyboardInterrupt:
+        print("Program terminated by user.")
diff --git a/scripts/cumulative_cell_vee.py b/scripts/cumulative_cell_vee.py
new file mode 100644
index 0000000..b6e07d5
--- /dev/null
+++ b/scripts/cumulative_cell_vee.py
@@ -0,0 +1,78 @@
+import numpy as np, os, csv
+import matplotlib.pyplot as plt
+import matplotlib.ticker as mtick
+
+from squish import Simulation, DomainParams, ordered
+from squish.common import OUTPUT_DIR
+from script_tools import (
+    RC_SETTINGS,
+    get_args,
+    get_data,
+    get_simulation_data,
+    get_ordered_data,
+)
+
+NAME = "Cumulative-CellVEE"
+ALPHA = 1.0
+
+
+def main():
+    sims_path, regen = get_args(
+        "Anti-cumulative distribution of VEE and percent of equilibria for fixed alpha",
+        "folders that contains various N simulations to plot",
+    )
+
+    sim, frames = Simulation.load(sims_path)
+    vees = np.linspace(-3, 3, 10000)
+    e_hex = ordered.e_hex(sim.domain)
+    energies = {"all": np.empty(0, dtype=float)}
+    counts = {}
+    for frame in frames:
+        energy = frame["stats"]["site_energies"] - e_hex
+        defects = np.count_nonzero(frame["stats"]["site_edge_count"] != 6)
+        if defects not in energies:
+            energies[defects] = np.empty(0, dtype=float)
+        energies[defects] = np.append(energies[defects], energy)
+        energies["all"] = np.append(energies["all"], energy)
+
+    all_count = None
+    for defect, energy in energies.items():
+        count = np.empty(vees.shape, dtype=float)
+        for i, vee in enumerate(vees):
+            count[i] = np.count_nonzero(energy >= vee)
+        count = 100 * count / len(energy)
+        if defect == "all":
+            all_count = count
+        else:
+            counts[defect] = count
+
+    plt.rcParams.update(RC_SETTINGS)
+
+    fig = plt.figure(figsize=(15, 15))
+    gs = fig.add_gridspec(1, 1)
+    ax = fig.add_subplot(gs[0])
+
+    for defect, count in sorted(counts.items()):
+        if defect == min(counts.keys()):
+            ax.plot(vees, count, label=fr"$\mathrm{{D}}={defect}$", color="C0")
+        elif defect == max(counts.keys()):
+            ax.plot(vees, count, label=fr"$\mathrm{{D}}={defect}$", color="C1")
+        else:
+            ax.plot(vees, count, linestyle="dotted", alpha=0.3)
+    ax.plot(vees, all_count, label="All", color="black")
+
+    ax.set_xlim(-2.5, 2)
+    ax.yaxis.set_major_formatter(mtick.PercentFormatter())
+
+    ax.set_xlabel(r"VEE")
+    ax.set_ylabel("Percent of Voronoi Regions")
+
+    ax.grid(zorder=0)
+    ax.legend()
+
+    fig.savefig(OUTPUT_DIR / (NAME + ".png"))
+    print(f"Wrote to {OUTPUT_DIR / (NAME + '.png')}")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/scripts/cumulative_vee.py b/scripts/cumulative_vee.py
index ba91d7f..16eb19a 100644
--- a/scripts/cumulative_vee.py
+++ b/scripts/cumulative_vee.py
@@ -1,4 +1,4 @@
-import numpy as np, os
+import numpy as np, os, csv
 import matplotlib.pyplot as plt
 import matplotlib.ticker as mtick
 
@@ -48,6 +48,10 @@ def main():
     gs = fig.add_gridspec(1, 1)
     ax = fig.add_subplot(gs[0])
 
+    my_cool_data = [["VEE", 61, 67, 73, 81, 84, 100]]
+    for vee in np.linspace(0, 0.06, 10000):
+        my_cool_data.append([vee])
+
     for j, package in enumerate(packages):
         data, ordered_data, domain = package
         e_hex = ordered.e_hex(domain)
@@ -66,6 +70,9 @@ def main():
             counts[i] = np.count_nonzero(energies >= vee)
         counts = 100 * counts / len(energies)
 
+        for i, count in enumerate(counts):
+            my_cool_data[i + 1].append(count)
+
         ax.plot(
             100 * vees[: index + 1],
             counts[: index + 1],
@@ -80,6 +87,10 @@ def main():
             color=f"C{j}",
         )
 
+    with open("cumulative-vee.csv", "w") as csvfile:
+        writer = csv.writer(csvfile)
+        writer.writerows(my_cool_data)
+
     ax.set_xlim(0, 6.3)
     ax.yaxis.set_major_formatter(mtick.PercentFormatter())
 
diff --git a/scripts/cumulative_vee_2.py b/scripts/cumulative_vee_2.py
new file mode 100644
index 0000000..7d59e87
--- /dev/null
+++ b/scripts/cumulative_vee_2.py
@@ -0,0 +1,213 @@
+import numpy as np, os, csv
+import matplotlib.pyplot as plt
+import matplotlib.ticker as mtick
+from scipy.optimize import curve_fit
+
+from squish import Simulation, DomainParams, ordered
+from squish.common import OUTPUT_DIR
+from script_tools import (
+    RC_SETTINGS,
+    get_args,
+    get_data,
+    get_simulation_data,
+    get_ordered_data,
+)
+
+NAME = "Cumulative-VEE-AcrossN"
+NAME2 = "Cumulative-VEE-RelError"
+NAME3 = "Cumulative-VEE-Alphas"
+
+
+def main():
+    sims_path, regen = get_args(
+        "Anti-cumulative distribution of VEE and percent of equilibria for fixed alpha",
+        "folders that contains various N simulations to plot",
+    )
+
+    vees = np.linspace(0, 0.06, 10000)
+
+    def f(path):
+        files = list(path.iterdir())
+        files = [f for f in files if f.is_dir()]
+        data = {}
+
+        for i, fol in enumerate(files):
+            sim, frames = Simulation.load(fol)
+            e_hex = ordered.e_hex(sim.domain)
+            energies = []
+            for frame in frames:
+                energies.append(frame["energy"] / sim.domain.n - e_hex)
+            energies = np.array(energies)
+
+            counts = np.empty(vees.shape, dtype=float)
+            for j, vee in enumerate(vees):
+                counts[j] = np.count_nonzero(energies >= vee)
+            counts = 100 * counts / len(energies)
+
+            data[sim.domain.n] = counts
+
+            hashes = int(20 * i / len(files))
+            print(
+                f'Loading simulations... |{"#"*hashes}{" "*(20-hashes)}|'
+                + f" {i+1}/{len(files)} simulations loaded.",
+                flush=True,
+                end="\r",
+            )
+
+        print(flush=True)
+        return data
+
+    all_data = {}
+    alpha_files = list(sims_path.iterdir())
+    alpha_files = [f for f in alpha_files if f.is_dir()]
+    for fol in alpha_files:
+        all_data[float(fol.name)] = get_data(
+            fol / "Cumulative2.pkl", f, args=(fol,), regen=regen
+        )
+
+    data = all_data[1.0]
+    points = []
+    avgs = {}
+    # print((data[106] + data[104]) / 2 - data[105])
+    for n in range(105, 396):
+        arr = np.zeros(data[100].shape)
+        for i in range(1, 6):
+            arr += data[n - i] + data[n + i]
+        arr += data[n]
+        arr /= 11
+        avgs[n] = arr
+        points.append(np.linalg.norm(data[n] - arr) / np.linalg.norm(arr))
+
+    plt.rcParams.update(RC_SETTINGS)
+
+    fig = plt.figure(figsize=(15, 15))
+    gs = fig.add_gridspec(1, 1)
+    ax = fig.add_subplot(gs[0])
+
+    ax.plot(100 * vees, avgs[395], color="black", label=r"Average", zorder=50)
+
+    for n, counts in sorted(data.items()):
+        if n in range(390, 401):
+            if n == 392:
+                ax.plot(100 * vees, counts, label="N=392", color="C0")
+            elif n == 397:
+                ax.plot(100 * vees, counts, label="N=397", color="C1")
+            elif n == 395:
+                ax.plot(100 * vees, counts, color="black", linestyle="dashed")
+            else:
+                ax.plot(
+                    100 * vees,
+                    counts,
+                    label="_nolegend_",
+                    linestyle="dashed",
+                    alpha=0.5,
+                )
+
+    ax.set_xlim(0, 6.3)
+    ax.yaxis.set_major_formatter(mtick.PercentFormatter())
+
+    ax.set_xlabel(r"VEE $\left[\times 10^{2}\right]$")
+    ax.set_ylabel("Percent of Equilibria")
+
+    ax.grid(zorder=0)
+    ax.legend()
+
+    fig.savefig(OUTPUT_DIR / (NAME + ".png"))
+    print(f"Wrote to {OUTPUT_DIR / (NAME + '.png')}")
+
+    fig = plt.figure(figsize=(30, 8))
+    gs = fig.add_gridspec(1, 1)
+
+    ax = fig.add_subplot(gs[0])
+    ns = np.array(range(105, 396))
+    ax.scatter(ns, points)
+
+    log_slope, _ = np.polyfit(np.log10(ns), np.log10(points), 1)
+    print(log_slope)
+
+    def g(x, k):
+        return k * (x ** log_slope)
+
+    # return a * np.e ** (-b * (x - c))
+
+    params, covs = curve_fit(g, ns, points, p0=[5000])
+    ax.plot(np.arange(50, 500), g(np.arange(50, 500), *params), color="C1")
+
+    # with open("cumul.csv", "w") as csvfile:
+    #    csvwriter = csv.writer(csvfile)
+    #    csvwriter.writerows(zip(ns, points))
+    print(params)
+    ax.set_xlim(95, 405)
+    ax.set_ylim(0, 0.24)
+    ax.set_xlabel("N")
+
+    ax.grid(zorder=0)
+    # ax.legend()
+
+    fig.savefig(OUTPUT_DIR / (NAME2 + ".png"))
+    print(f"Wrote to {OUTPUT_DIR / (NAME2 + '.png')}")
+
+    fig = plt.figure(figsize=(15, 15))
+    gs = fig.add_gridspec(1, 1)
+    ax = fig.add_subplot(gs[0])
+
+    alpha_avgs = {}
+    for alpha, data in all_data.items():
+        if len(data.values()) == 0 or alpha == 1.0:
+            continue
+        alpha_avgs[alpha] = sum(data.values()) / 11
+
+    for alpha, data in all_data.items():
+        if alpha == 1.0:
+            continue
+        for n, counts in data.items():
+            # ax.plot(100 * vees, counts, linestyle="dotted", alpha=0.4)
+            continue
+
+    pad = 100
+    alpha_prob_dens = {}
+    for alpha, avg in alpha_avgs.items():
+        if alpha == 1.0:
+            continue
+        mov_avgs = np.array(
+            [np.mean(avg[max(0, i - pad) : i + pad]) for i in range(len(avg))]
+        )
+        d_avgs = np.gradient(mov_avgs, vees[1] - vees[0])
+        alpha_prob_dens[alpha] = d_avgs / (np.sum(d_avgs) * (vees[1] - vees[0]))
+
+    a1_prob_dens = {}
+    for n, avg in avgs.items():
+        mov_avgs = np.array(
+            [np.mean(avg[max(0, i - pad) : i + pad]) for i in range(len(avg))]
+        )
+
+        d_avgs = np.gradient(mov_avgs, vees[1] - vees[0])
+        a1_prob_dens[n] = d_avgs / (np.sum(d_avgs) * (vees[1] - vees[0]))
+
+    for alpha, prob_dens in sorted(alpha_prob_dens.items()):
+        if alpha in [0.3, 0.5, 0.7, 0.9, 1.0]:
+            continue
+            ax.plot(100 * vees, prob_dens, label=fr"$\alpha={alpha:.1f}$")
+
+    for n, prob_dens in sorted(a1_prob_dens.items()):
+        if n in [255, 315, 335, 355, 375]:
+            ax.plot(100 * vees, prob_dens, label=f"N={n}")
+        if n == 395:
+            ax.plot(100 * vees, prob_dens, label=f"N={n}", color="black")
+            # ax.plot(100 * vees, prob_dens, label=r"$\alpha = 1.0$", zorder=50, color="black")
+    # ax.plot(100 * vees, avgs[395], zorder=50, color="black")
+
+    ax.set_xlim(0, 6.3)
+    # ax.yaxis.set_major_formatter(mtick.PercentFormatter())
+
+    ax.set_xlabel(r"VEE $\left[\times 10^{2}\right]$")
+    # ax.set_ylabel("Percent of Equilibria")
+
+    ax.grid(zorder=0)
+    ax.legend()
+    fig.savefig(OUTPUT_DIR / (NAME3 + ".png"))
+    print(f"Wrote to {OUTPUT_DIR / (NAME3 + '.png')}")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/scripts/defect_scatter.py b/scripts/defect_scatter.py
index fc7e8da..be1195c 100644
--- a/scripts/defect_scatter.py
+++ b/scripts/defect_scatter.py
@@ -27,14 +27,22 @@ def main():
     defects = np.array(defects)
     energy = 100 * np.array(energy)
 
+    counts = {}
+    for i in range(70):
+        counts[i] = np.count_nonzero(defects == i)
+    print(counts)
+    for k, v in counts.items():
+        if k % 2 == 1 and v > 0:
+            print(k, v)
+
     plt.rcParams.update(RC_SETTINGS)
-    fig = plt.figure(figsize=(18, 15))
+    fig = plt.figure(figsize=(15, 15))
     gs = fig.add_gridspec(1, 1)
     ax = fig.add_subplot(gs[0])
 
     m, b = np.polyfit(defects, energy, 1)
-    ax.scatter(defects, energy, color="C0", marker="*", s=100)
-    ax.plot(np.arange(0, 64), np.arange(0, 64) * m + b, zorder=3, color="C1")
+    ax.scatter(defects, energy, color="C0", marker="*", s=70, zorder=2)
+    ax.plot(np.arange(0, 65), np.arange(0, 65) * m + b, zorder=3, color="C1")
 
     ax.scatter(0, b, color="C2", s=500, marker="*", zorder=50)
 
@@ -49,10 +57,10 @@ def main():
     ax.set_xlim(0, 64)
     ax.set_xticks(np.arange(0, 65, 8))
 
-    ax.set_ylim(0, 3.6)
+    ax.set_ylim(0, 3.7)
     ax.set_yticks(np.arange(0, 3.6, 0.5))
 
-    ax.set_xlabel("Number of Defects")
+    ax.set_xlabel("$\mathrm{D}$")
     ax.set_ylabel(r"VEE $\left[\times 10^{2}\right]$")
     ax.grid(zorder=0)
 
diff --git a/scripts/frustration_alpha.py b/scripts/frustration_alpha.py
new file mode 100644
index 0000000..bb5485a
--- /dev/null
+++ b/scripts/frustration_alpha.py
@@ -0,0 +1,185 @@
+import numpy as np, os, math
+import matplotlib.pyplot as plt
+from multiprocessing import Pool, cpu_count
+
+from squish import Simulation, ordered, DomainParams
+from squish.common import OUTPUT_DIR
+from script_tools import RC_SETTINGS, get_args, get_data, get_ordered_data, format_data
+
+NAME = "DefectEnergy"
+NAME2 = "Frustration"
+
+
+def proc(path):
+    sim, frames = Simulation.load(path)
+    domain = sim.domain
+    e_hex = ordered.e_hex(domain)
+
+    defects, energy = [], []
+    for frame in frames:
+        # if np.var(frame["stats"]["avg_radius"]) > 1e-8:
+        defects.append(np.count_nonzero(frame["stats"]["site_edge_count"] != 6))
+        energy.append(frame["energy"] / domain.n - e_hex)
+
+    configs = ordered.configurations(domain)
+
+    avg_defects = sum(defects) / len(defects)
+    avg_vee = sum(energy) / len(energy)
+    m, b = np.polyfit(defects, energy, 1)
+
+    alpha = domain.w / domain.h
+
+    return [alpha, tuple(domain), avg_defects, avg_vee, m, b]
+
+
+def main():
+    sims_path, regen = get_args(
+        "Generates graph for Frustration at fixed N across alpha",
+        "folder that contains simulations at various alpha for fixed N",
+    )
+
+    def f():
+        data = {}
+
+        files = list(sims_path.iterdir())
+        files = [x for x in files if x.is_dir()]
+        with Pool(cpu_count()) as pool:
+            for i, res in enumerate(pool.imap_unordered(proc, files)):
+                data[res[0]] = res[1:]
+
+                hashes = int(21 * i / len(files))
+                print(
+                    f'Processed alpha={res[0]:.3f} |{"#"*hashes}{" "*(20-hashes)}|'
+                    + f" {i+1}/{len(files)} simulations processed.",
+                    flush=True,
+                    end="\r",
+                )
+
+            print(flush=True)
+
+        return format_data(
+            data,
+            key_name="alpha",
+            col_names=[
+                "Domain",
+                "Average Defects",
+                "Average VEE",
+                "Energy Per Defect",
+                "Ground VEE",
+            ],
+        )
+
+    plt.rcParams.update(RC_SETTINGS)
+    data = get_data(sims_path / (NAME + ".pkl"), f, regen=regen)
+    domain = DomainParams(*data["Domain"][0])
+    # ordered_data = get_data(
+    #     OUTPUT_DIR / "OrderedCache" / f"{domain.n}.pkl",
+    #     get_ordered_data,
+    #     args=(domain, data["alpha"]),
+    #     regen=regen,
+    # )
+    # min_order = []
+    # for i, energies in enumerate(ordered_data["Energy"]):
+    # min_order.append(energies[0])
+
+    alphas, avg_defects, avg_vees, epds, vee0s = (
+        data["alpha"],
+        data["Average Defects"],
+        100 * data["Average VEE"],
+        100 * data["Energy Per Defect"],
+        100 * data["Ground VEE"],
+        # 100 * (np.array(min_order) / domain.n - ordered.e_hex(domain)),
+    )
+
+    fig = plt.figure(figsize=(15, 15))
+    gs = fig.add_gridspec(1, 1)
+    ax = fig.add_subplot(gs[0])
+
+    m0, b0 = np.polyfit(alphas, avg_defects, 1)
+    # m1, b1 = np.polyfit(ns, avg_defects * epds, 1)
+
+    lns1 = ax.plot(
+        alphas,
+        avg_defects,
+        color="C0",
+        alpha=0.5,
+        label=r"$\langle \mathrm{D} \rangle$",
+    )
+    ax.plot(alphas, m0 * alphas + b0, color="C0", linestyle="dashed")
+
+    lns2 = ax.plot(
+        alphas, 100 * epds, color="C1", alpha=0.5, label=r"$\zeta_1 \times 10^{4}$"
+    )
+
+    lns3 = ax.plot(
+        alphas,
+        10 * avg_defects * epds,
+        color="C2",
+        alpha=0.5,
+        label=r"$\mathcal{F} \times 10^{3}$",
+    )
+    # ax3.plot(ns, 100*(m1 * ns + b1) / 10, color="C2", linestyle="dashed")
+    ax.set_xlim(0.3, 1.0)
+
+    # ax.set_ylim(3, 37)
+    # ax.set_yticks(np.arange(5, 40, 5))
+
+    # ax2.set_ylim(-3 * 0.4, 18 + 3 * 0.4)
+    # ax2.set_yticks(np.arange(0, 20, 3))
+
+    # ax3.set_ylim(-3 * 0.5, 18 + 3 * 0.4)
+    # ax3.set_yticks([])
+    # ax3.spines["right"].set_visible(False)
+    # ax3.spines.right.set_position(("axes", 1.11))
+
+    # lns = lns1 + lns2 + lns3
+    # labs = [l.get_label() for l in lns]
+    ax.grid(zorder=0)
+    # ax.legend(lns, labs, loc="lower right")
+    ax.legend()
+    ax.set_xlabel("N")
+    # ax.set_ylabel(r"$\langle \mathrm{D} \rangle$")
+    # ax2.set_ylabel("VEE")
+    # ax2.set_ylabel(r"$\zeta \left[\times 10^{4} \right]$", color="C1")
+    # ax3.set_ylabel(r"Defect Energy $\left[\times 10^{3} \right]$", color="C2")
+
+    out = OUTPUT_DIR / (NAME + f" - N{domain.n}.png")
+    fig.savefig(out)
+    print(f"Wrote to {out}")
+
+    # return
+    fig = plt.figure(figsize=(20, 15))
+    gs = fig.add_gridspec(1, 1)
+    ax = fig.add_subplot(gs[0])
+
+    ax.plot(alphas, vee0s, label="$\zeta_0$", color="C5")
+    ax.plot(alphas, avg_defects * epds, label=r"$\mathcal{F}$", color="C2")
+    ax.plot(alphas, avg_vees, label=r"$\langle \mathrm{VEE} \rangle$", color="C4")
+    # ax.plot(alphas, min_order, label="Minimum Ordered", color="C7", alpha=0.5)
+
+    print(np.linalg.norm(avg_vees - (vee0s + avg_defects * epds)))
+
+    ax.grid(zorder=0)
+    ax.legend()
+
+    ax.set_xlim(0.3, 1)
+    # ax.set_xticks([round(w, 2) for w in alphas[::10]])
+    # ax.set_xticklabels([f"{round(w, 3):.2f}" for w in alphas[::10]], rotation=90)
+
+    # ax.set_ylim(-0.1, 9.2)
+    # ax.set_yticks(np.arange(0, 9.6, 1))
+
+    ax.set_xlabel("Aspect Ratio")
+    ax.set_ylabel(r"VEE $\left[\times 10^{2}\right]$")
+
+    out = OUTPUT_DIR / (NAME2 + f" - N{domain.n}.png")
+    fig.savefig(out)
+    print(f"Wrote to {out}")
+
+
+if __name__ == "__main__":
+    os.environ["QT_log10GING_RULES"] = "*=false"
+    try:
+        main()
+    except KeyboardInterrupt:
+        print("Program terminated by user.")
diff --git a/scripts/frustration.py b/scripts/frustration_n.py
similarity index 89%
rename from scripts/frustration.py
rename to scripts/frustration_n.py
index fe030e0..c31f607 100644
--- a/scripts/frustration.py
+++ b/scripts/frustration_n.py
@@ -88,7 +88,7 @@ def main():
         100 * data["Minimum Ordered VEE"],
     )
 
-    fig = plt.figure(figsize=(18, 15))
+    fig = plt.figure(figsize=(15, 15))
     gs = fig.add_gridspec(1, 1)
     ax = fig.add_subplot(gs[0])
     # ax2 = ax.twinx()
@@ -103,17 +103,18 @@ def main():
     ax.plot(ns, m0 * ns + b0, color="C0", linestyle="dashed")
 
     lns2 = ax.plot(
-        ns, 100 * epds, color="C1", alpha=0.5, label=r"$\zeta \times 10^{4}$"
+        ns, 100 * epds, color="C1", alpha=0.5, label=r"$\zeta_1 \times 10^{4}$"
     )
 
     lns3 = ax.plot(
         ns,
-        100 * avg_defects * epds / 10,
+        10 * avg_defects * epds,
         color="C2",
         alpha=0.5,
         label=r"$\mathcal{F} \times 10^{3}$",
     )
     # ax3.plot(ns, 100*(m1 * ns + b1) / 10, color="C2", linestyle="dashed")
+    ax.set_xlim(93, 407)
 
     ax.set_ylim(3, 37)
     ax.set_yticks(np.arange(5, 40, 5))
@@ -140,23 +141,25 @@ def main():
     fig.savefig(OUTPUT_DIR / (NAME + ".png"))
     print(f"Wrote to {OUTPUT_DIR / (NAME + '.png')}")
 
-    return
-    fig = plt.figure(figsize=(30, 15))
+    # return
+    fig = plt.figure(figsize=(20, 15))
     gs = fig.add_gridspec(1, 1)
     ax = fig.add_subplot(gs[0])
 
-    ax.plot(ns, vee0s, label="$\mathrm{VEE}_0$")
-    ax.plot(ns, avg_defects * epds, label=r"$\mathcal{F}$")
-    ax.plot(ns, avg_vees, label=r"$\langle \mathrm{VEE} \rangle$")
-    ax.plot(ns, min_order, label="Minimum Ordered", color="C5", alpha=0.5)
+    ax.plot(ns, vee0s, label="$\zeta_0$", color="C5")
+    ax.plot(ns, avg_defects * epds, label=r"$\mathcal{F}$", color="C2")
+    ax.plot(ns, avg_vees, label=r"$\langle \mathrm{VEE} \rangle$", color="C4")
+    ax.scatter(ns, min_order, label="Minimum Ordered", color="C7", alpha=0.8)
 
     print(np.linalg.norm(avg_vees - (vee0s + avg_defects * epds)))
 
     ax.grid(zorder=0)
     ax.legend()
 
-    ax.set_ylim(0, 3)
-    ax.set_yticks(np.arange(0, 3.1, 0.5))
+    ax.set_xlim(95, 405)
+
+    ax.set_ylim(-0.1, 9.2)
+    ax.set_yticks(np.arange(0, 9.6, 1))
 
     ax.set_xlabel("N")
     ax.set_ylabel(r"VEE $\left[\times 10^{2}\right]$")
diff --git a/scripts/min_order_vee.py b/scripts/min_order_vee.py
new file mode 100644
index 0000000..e334cd7
--- /dev/null
+++ b/scripts/min_order_vee.py
@@ -0,0 +1,91 @@
+import numpy as np, os, math
+import matplotlib.pyplot as plt
+from multiprocessing import Pool, cpu_count
+
+from squish import Simulation, ordered, DomainParams
+from squish.common import OUTPUT_DIR
+from script_tools import RC_SETTINGS, get_args, get_data, format_data
+
+NAME = "MinOrderVEE"
+
+
+def get_min_vee(n):
+    domain = DomainParams(n, math.sqrt(n), math.sqrt(n), 4)
+    e_hex = ordered.e_hex(domain)
+
+    configs = ordered.configurations(domain)
+    min_vee, min_config = 1e10, None
+
+    for config in configs:
+        try:
+            rbar = ordered.avg_radius(domain, config)
+        except:
+            print(tuple(domain), config)
+
+        vee = (
+            2 * domain.w * domain.h
+            + 2 * math.pi * domain.n * (domain.r ** 2 - 2 * domain.r * rbar)
+        ) / domain.n - e_hex
+
+        if vee < min_vee and vee >= 0:
+            min_vee = vee
+            min_config = config
+
+    return (domain.n, min_vee, min_config)
+
+
+def main():
+    def f():
+        data = {}
+        ns = list(range(1000, 5001))
+        with Pool(cpu_count()) as pool:
+            for i, res in enumerate(pool.imap_unordered(get_min_vee, ns)):
+                data[res[0]] = res[1:]
+
+                hashes = int(21 * i / len(ns))
+                print(
+                    f'Processed N={res[0]} |{"#"*hashes}{" "*(20-hashes)}|'
+                    + f" {i+1}/{len(ns)} simulations processed.",
+                    flush=True,
+                    end="\r",
+                )
+
+            print(flush=True)
+
+        return format_data(
+            data, key_name="N", col_names=["Minimum Ordered VEE", "Config"]
+        )
+
+    plt.rcParams.update(RC_SETTINGS)
+    data = get_data(OUTPUT_DIR / "OrderedCache" / (NAME + ".pkl"), f)
+    ns, min_order, min_config = (
+        data["N"],
+        100 * data["Minimum Ordered VEE"],
+        data["Config"],
+    )
+
+    print(min_order, min_config)
+
+    fig = plt.figure(figsize=(20, 15))
+    gs = fig.add_gridspec(1, 1)
+    ax = fig.add_subplot(gs[0])
+
+    ax.scatter(ns, min_order)
+
+    ax.grid(zorder=0)
+
+    # ax.set_xlim(1000, 5000)
+
+    ax.set_xlabel("N")
+    ax.set_ylabel(r"VEE $\left[\times 10^{2}\right]$")
+
+    fig.savefig(OUTPUT_DIR / (NAME + ".png"))
+    print(f"Wrote to {OUTPUT_DIR / (NAME + '.png')}")
+
+
+if __name__ == "__main__":
+    os.environ["QT_log10GING_RULES"] = "*=false"
+    try:
+        main()
+    except KeyboardInterrupt:
+        print("Program terminated by user.")
diff --git a/scripts/ordered_rbar.py b/scripts/ordered_rbar.py
index 08edad8..9c5b66e 100644
--- a/scripts/ordered_rbar.py
+++ b/scripts/ordered_rbar.py
@@ -6,7 +6,7 @@ from squish import Simulation, ordered
 from squish.common import OUTPUT_DIR
 from script_tools import RC_SETTINGS, get_data, format_data
 
-NAME = "OrdereredRBar"
+NAME = "OrderedRBar"
 
 
 def main():
@@ -34,7 +34,7 @@ def main():
 
     plt.rcParams.update(RC_SETTINGS)
 
-    fig = plt.figure(figsize=(20, 10))
+    fig = plt.figure(figsize=(18, 14.4))
     gs = fig.add_gridspec(1, 1)
     ax = fig.add_subplot(gs[0])
 
@@ -53,6 +53,8 @@ def main():
         color="black",
     )
 
+    ax.set_ylim(0.51, 0.615)
+    ax.set_yticks(np.arange(0.52, 0.62, 0.03))
     ax.set_ylabel(r"$\bar{r}$")
     ax.grid(zorder=0)
     ax.legend()
diff --git a/scripts/ordered_scatter.py b/scripts/ordered_scatter.py
index e8b153d..095ca4a 100644
--- a/scripts/ordered_scatter.py
+++ b/scripts/ordered_scatter.py
@@ -37,9 +37,14 @@ def main():
     e_hex = ordered.e_hex(DomainParams(args.n, 1, 1, args.r))
     for alpha, energies in zip(ordered_data["alpha"], ordered_data["Energy"]):
         ax.scatter(
-            [alpha] * len(energies), 100 * (energies / args.n - e_hex), color="C0"
+            [alpha] * len(energies), 100 * (energies / args.n - e_hex), color="C0", s=25
         )
 
+    hex_ratios = ordered.hexagon_alpha(args.n)
+    ax.scatter(
+        hex_ratios, [0] * len(hex_ratios), color="C2", s=200, marker="H", zorder=10
+    )
+
     ax.set_xlim(0.3, 1)
     ax.set_xticks(np.arange(0.3, 1.01, 0.1))
 
diff --git a/scripts/perturbations.py b/scripts/perturbations.py
index a42b16c..ad83bb5 100644
--- a/scripts/perturbations.py
+++ b/scripts/perturbations.py
@@ -2,20 +2,13 @@ import argparse, numpy as np, os
 from pathlib import Path
 import matplotlib.pyplot as plt
 
-from squish import DomainParams, Simulation
+from squish import DomainParams, Simulation, ordered
 from squish.common import OUTPUT_DIR
 from script_tools import RC_SETTINGS, get_data, format_data
 
 NAME = "Perturbations"
 
 
-def toroidal_distance(domain: DomainParams, x: np.ndarray, y: np.ndarray) -> float:
-    dim = np.array([domain.w, domain.h])
-    absdist = np.abs(y - x)
-    wrap = dim * (absdist >= np.array([domain.w, domain.h]) / 2)
-    return np.linalg.norm(wrap - absdist)
-
-
 def main():
     parser = argparse.ArgumentParser(
         description="Graphs perturbation graphs for a collection of simulations."
@@ -60,8 +53,10 @@ def main():
                 )
 
                 data[delta]["norm"].append(
-                    toroidal_distance(
-                        end_sim.domain, adjusted, end_sim.frames[0].site_arr
+                    np.linalg.norm(
+                        ordered.toroidal_distance(
+                            end_sim.domain, adjusted, end_sim.frames[0].site_arr
+                        )
                     )
                 )
                 data[delta]["time"].append(sim.step_size * i)
@@ -83,7 +78,7 @@ def main():
             label=f"k = {data[delta]['k']}",
         )
 
-    ax.set_title(r"Relaxation of Perturbations")
+    # ax.set_title(r"Relaxation of Perturbations")
 
     ax.set_xlim([0, 60])
     ax.set_yscale("log")
diff --git a/scripts/probability_of_disorder.py b/scripts/probability_of_disorder.py
index e6efa01..76357f7 100644
--- a/scripts/probability_of_disorder.py
+++ b/scripts/probability_of_disorder.py
@@ -90,15 +90,15 @@ def main():
 
     ax.set_xlabel("Aspect Ratio")
     ax.set_ylabel(r"VEE $\left[\times 10^{2}\right]$", color="C2")
-    ax2.set_ylabel("Percent of Disordered Equilibria", color="C4")
+    ax2.set_ylabel("POD", color="C4")
 
     # ax.legend(loc=(0.23, 0.5))
     ax.grid(zorder=0)
 
     props = dict(boxstyle="round", facecolor="white", alpha=0.8, zorder=20)
     ax.text(
-        0.873,
-        0.97,
+        0.85,
+        0.92,
         f"N={domain.n}",
         transform=ax.transAxes,
         verticalalignment="top",
diff --git a/scripts/rbar.py b/scripts/rbar.py
new file mode 100644
index 0000000..6653896
--- /dev/null
+++ b/scripts/rbar.py
@@ -0,0 +1,85 @@
+import argparse, numpy as np, os
+from multiprocessing import Pool, cpu_count
+import matplotlib.pyplot as plt
+
+from squish import Simulation, ordered
+from squish.common import OUTPUT_DIR
+from script_tools import RC_SETTINGS, get_args, get_data, format_data
+
+NAME = "RBar"
+
+
+def main():
+    sims_path, regen = get_args(
+        "Variance of average radius for as N gets larger",
+        "folders that contains various N simulations to plot",
+    )
+
+    def f():
+        files = list(sims_path.iterdir())
+        files = [f for f in files if f.is_dir()]
+        data = {}
+
+        for i, fol in enumerate(files):
+            sim, frames = Simulation.load(fol)
+
+            variances = []
+            for frame in frames:
+                rbars = np.sort(frame["stats"]["avg_radius"])
+                variances.append(np.var(rbars))
+            avg_variance = sum(variances) / len(variances)
+
+            data[sim.domain.n] = [avg_variance]
+
+            hashes = int(20 * i / len(files))
+            print(
+                f'Loading simulations... |{"#"*hashes}{" "*(20-hashes)}|'
+                + f" {i+1}/{len(files)} simulations loaded.",
+                flush=True,
+                end="\r",
+            )
+
+        print(flush=True)
+        return format_data(data, key_name="N", col_names=["Average Variance"])
+
+    data = get_data(sims_path / (NAME + ".pkl"), f, regen=regen)
+
+    plt.rcParams.update(RC_SETTINGS)
+
+    fig = plt.figure(figsize=(18, 14.4))
+    gs = fig.add_gridspec(1, 1)
+    ax = fig.add_subplot(gs[0])
+
+    primes = [
+        i
+        for i in range(min(data["N"]), max(data["N"]))
+        if len(ordered.divisors(i)) == 2
+    ]
+
+    prime_points = []
+    for i, n in enumerate(data["N"]):
+        if n in primes:
+            prime_points.append(data["Average Variance"][i])
+
+    ax.plot(data["N"], 1e5 * data["Average Variance"])
+    ax.scatter(primes, 1e5 * np.array(prime_points), marker="o", s=60)
+    ax.set_xlim(95, 405)
+    ax.set_ylim(0, 11)
+    # ax.set_ylim(0, args.max_n)
+
+    ax.set_xlabel("N")
+    ax.set_ylabel(r"$\sigma \left[\times 10^5 \right]$")
+    ax.grid(zorder=0)
+    # ax.legend()
+
+    fig.savefig(OUTPUT_DIR / (NAME + ".png"))
+
+    print(f"Wrote to {OUTPUT_DIR / (NAME + '.png')}")
+
+
+if __name__ == "__main__":
+    os.environ["QT_log10GING_RULES"] = "*=false"
+    try:
+        main()
+    except KeyboardInterrupt:
+        print("Program terminated by user.")
diff --git a/scripts/script_tools.py b/scripts/script_tools.py
index f04042f..d439d01 100644
--- a/scripts/script_tools.py
+++ b/scripts/script_tools.py
@@ -26,7 +26,7 @@ RC_SETTINGS = {
     "font.family": "cm",
     "font.size": 40,
     "text.usetex": True,
-    "text.latex.preamble": r"\usepackage{amsmath}",
+    "text.latex.preamble": r"\usepackage{amsmath, amsfonts}",
     "figure.constrained_layout.use": True,
 }
 
diff --git a/scripts/torus.py b/scripts/torus_img.py
similarity index 67%
rename from scripts/torus.py
rename to scripts/torus_img.py
index b82c4e9..46d0a76 100644
--- a/scripts/torus.py
+++ b/scripts/torus_img.py
@@ -2,40 +2,40 @@ import numpy as np
 import matplotlib.pyplot as plt
 
 
-def torus_mesh(n, c, a):
+def torus_mesh(n, r, R):
     theta = np.linspace(0, 2 * np.pi, n)
     phi = np.linspace(0, 2 * np.pi, n)
     theta, phi = np.meshgrid(theta, phi)
-    x = (c + a * np.cos(theta)) * np.cos(phi)
-    y = (c + a * np.cos(theta)) * np.sin(phi)
-    z = a * np.sin(theta)
+    x = (R + r * np.cos(theta)) * np.cos(phi)
+    y = (R + r * np.cos(theta)) * np.sin(phi)
+    z = r * np.sin(theta)
     return x, y, z
 
 
-def torus_line(n, c, a, u, v):
+def torus_line(n, r, R, u, v):
     phi = np.linspace(0, 2 * np.pi * u, n)
     theta = (v / u) * phi
     theta, phi = theta % (2 * np.pi), phi % (2 * np.pi)
 
-    x = (c + a * np.cos(theta)) * np.cos(phi)
-    y = (c + a * np.cos(theta)) * np.sin(phi)
-    z = a * np.sin(theta)
+    x = (R + r * np.cos(theta)) * np.cos(phi)
+    y = (R + r * np.cos(theta)) * np.sin(phi)
+    z = r * np.sin(theta)
 
     line = np.column_stack([x, y, z])
-    x, y, z = c * np.cos(phi), c * np.sin(phi), 0 * phi
+    x, y, z = R * np.cos(phi), R * np.sin(phi), 0 * phi
     norm = line - np.column_stack([x, y, z])
 
     return line, norm
 
 
-def torus_points(c, a, u, v, k):
+def torus_points(r, R, u, v, k):
     phi = np.linspace(0, 2 * np.pi * u, k + 1)
     theta = (v / u) * phi
     theta, phi = theta % (2 * np.pi), phi % (2 * np.pi)
 
-    x = (c + a * np.cos(theta)) * np.cos(phi)
-    y = (c + a * np.cos(theta)) * np.sin(phi)
-    z = a * np.sin(theta)
+    x = (R + r * np.cos(theta)) * np.cos(phi)
+    y = (R + r * np.cos(theta)) * np.sin(phi)
+    z = r * np.sin(theta)
 
     return x, y, z
 
@@ -54,18 +54,24 @@ def get_obscure_segs(arr):
 def main():
     fig = plt.figure(figsize=(10, 10), constrained_layout=True)
     ax = fig.add_subplot(projection="3d")
-    ax.set_zlim(-2, 2)
-    ax.set_axis_off()
+
+    ax.set_box_aspect(aspect=(1, 1, 1))
+
+    alpha = np.sqrt(3) * 2 / 5
+    size = 5
 
     azim, elev = 0, 40
-    c, a = 2, 1
-    u, v = 4, 1
-    k = 10
+    R, r = size, size * (1 - alpha) / alpha
+    u, v = 5, 2
+    k = 20
 
+    print(alpha, r / (r + R))
+
+    ax.set_zlim(-(r + R), r + R)
     ax.view_init(elev, azim)
 
     ax.plot_surface(
-        *torus_mesh(100, c, a), rcount=200, ccount=200, color="white", alpha=0.4
+        *torus_mesh(100, r, R), rcount=200, ccount=200, color="white", alpha=0.4
     )
 
     # Calculate when surface is facing or away and apply dotted vs solid lines.
@@ -73,7 +79,7 @@ def main():
     camera_vec = np.array(
         [np.sin(el) * np.cos(az), np.sin(el) * np.sin(az), np.cos(el)]
     )
-    line, norm = torus_line(500, c, a, u, v)
+    line, norm = torus_line(500, r, R, u, v)
     cond = np.dot(norm, camera_vec) >= 0
     bounds = get_obscure_segs(cond)
 
@@ -94,9 +100,10 @@ def main():
         ax.plot(lx, ly, lz, color="blue", linewidth=2, linestyle=style)
 
     # Display points.
-    ax.scatter(*torus_points(c, a, u, v, k), color="purple", s=50)
+    ax.scatter(*torus_points(r, R, u, v, k), color="purple", s=50)
     # ax.plot(*torus_line(200, 2, 1, 2, 1), color="orange", linewidth=3)
-    plt.savefig(f"Torus - N{k}({u},{v}).png", bbox_inches="tight")
+    plt.show()
+    # plt.savefig(f"Torus - N{k}({u},{v}).png", bbox_inches="tight")
 
 
 if __name__ == "__main__":
diff --git a/scripts/torus_mesh.py b/scripts/torus_mesh.py
new file mode 100644
index 0000000..6ca1e34
--- /dev/null
+++ b/scripts/torus_mesh.py
@@ -0,0 +1,330 @@
+import argparse, numpy as np, os
+from stl.stl import ASCII
+from stl import mesh
+from squish import DomainParams, Simulation, ordered
+from squish.common import OUTPUT_DIR
+from scipy.spatial import Voronoi
+
+SUBDIVISION_AMOUNT = 8
+
+
+def not_in_order(i, j):
+    if i < j:
+        return j - i != 1
+    else:
+        return i - j == 1
+
+
+def centroid(site, verts):
+    area, c = 0, 0
+    v = verts - site
+    # print(v)
+    for i in range(len(v)):
+        x, y = v[i], v[(i + 1) % len(v)]
+        a = np.cross(x, y)
+        area += a
+        c += a * (x + y)
+
+    return c / (3 * area) + site
+
+
+def flat_sheet(k):
+    x = np.linspace(-np.pi, np.pi, 101)
+    x = np.transpose([np.tile(a, len(a)), np.repeat(a, len(a))])
+    all_verts = np.hstack((a, np.zeros((len(a), 1))))
+
+    xx = np.array([[x for x in range(101 * 100) if (x - 100) % 101 > 0]]).T
+
+    b = np.hstack((xx, xx + 1, xx + 101))
+    c = np.hstack((xx + 1, xx + 102, xx + 101))
+    all_faces = np.concatenate((b, c))
+
+    return all_verts, all_faces
+
+
+def torus_transform(v, R, r):
+    new_verts = np.empty(v.shape)
+    # Rotate by 90 degrees, and then shift to range [0, 2\pi]
+    v[:, :2] = np.matmul(np.array([[0, -1], [1, 0]]), v[:, :2].T).T
+    v[:, :2] = v[:, :2] % (2 * np.pi)
+    # v[:, 0], v[:, 1] = v[:, 1], v[:, 0]
+
+    new_verts[:, 0] = (R + r * np.cos(v[:, 1])) * np.cos(v[:, 0])
+    new_verts[:, 1] = (R + r * np.cos(v[:, 1])) * np.sin(v[:, 0])
+    new_verts[:, 2] = r * np.sin(v[:, 1])
+    return new_verts
+
+
+def flat_region(c, verts, height):
+    m = len(verts)
+
+    centers = np.hstack((np.array([c, c]), np.array([[height], [0]])))
+    v_top = np.hstack((verts, height * np.ones((m, 1))))
+    v_bot = np.hstack((verts, np.zeros((m, 1))))
+
+    v = np.concatenate((centers, v_top, v_bot))
+
+    vi = np.atleast_2d(np.arange(m, dtype=int)).T
+    cent_top, cent_bot = np.zeros((m, 1), dtype=int), np.ones((m, 1), dtype=int)
+    vert_top_m, vert_top_p = vi + 2, (vi + 1) % m + 2
+    vert_bot_m, vert_bot_p = vi + m + 2, (vi + 1) % m + m + 2
+
+    top_face = np.hstack((cent_top, vert_top_m, vert_top_p))
+    bot_face = np.hstack((cent_bot, vert_bot_p, vert_bot_m))
+
+    f = np.concatenate((top_face, bot_face))
+    return v, f, 2
+
+
+def torus_region(c, verts, height):
+    sd = SUBDIVISION_AMOUNT
+    d = verts - c
+    m = len(verts)
+    # 1 + 6 + 12 ... verts, 1 + 3 + 5 ... faces.
+    v, f = (
+        np.empty((1 + m * sd * (sd - 1) // 2, 2)),
+        np.empty((m * (sd - 1) ** 2, 3), dtype=int),
+    )
+    v[0] *= 0
+
+    # Vertices are ordered by rings going out from the center.
+    v_inds = np.arange(-1, sd)
+    v_inds = 1 + m * v_inds * (v_inds + 1) // 2
+    v_inds[0] = 0
+
+    for i in range(1, sd):
+        fj, fk = m * (i - 1) ** 2, m * i ** 2
+        region_verts = [
+            np.linspace(i * d[j] / (sd - 1), i * d[(j + 1) % m] / (sd - 1), i + 1)[:-1]
+            for j in range(m)
+        ]
+        v[v_inds[i] : v_inds[i + 1]] = np.concatenate(region_verts)
+
+        faces = []
+        for j in range(m):
+            # v_off_inds = v_inds + j * np.arange(sd + 1)
+            r1, r2 = (
+                (j * (i - 1) + np.arange(i + 1)) % (m * (i - 1)) + v_inds[i - 1],
+                (j * i + np.arange(i + 2)) % (m * i) + v_inds[i],
+            )
+            r1, r2 = np.atleast_2d(r1).T, np.atleast_2d(r2).T
+            faces.append(np.hstack((r1[:-1], r2[:-2], r2[1:-1])))
+            faces.append(np.hstack((r1[:-2], r2[1:-2], r1[1:-1])))
+
+        f[fj:fk] = np.concatenate(faces)
+
+    v += c
+    v = np.hstack((v, height * np.ones((len(v), 1))))
+    return v, f, v_inds[sd - 1]
+
+    centers = np.array([[c[0], c[1], heights[i]]])
+    v_top = np.hstack((verts, heights[i] * np.ones((len(verts), 1))))
+
+    face_verts = np.concatenate((centers, v_top))
+
+
+def main():
+    parser = argparse.ArgumentParser(description="Generates 3D model for equilibrium.")
+    parser.add_argument(
+        "sims_path", metavar="sim_dir", help="simulation to obtain equilibrium from"
+    )
+    parser.add_argument("frame_num", metavar="FRAME_NUM", type=int, help="frame number")
+    parser.add_argument("size", metavar="SIZE", type=int, help="width in mm")
+    parser.add_argument(
+        "--torus", dest="torus", action="store_true", help="generates torus model"
+    )
+    args = parser.parse_args()
+
+    torus = args.torus
+    size = args.size
+
+    # Get desired frame and load.
+    sim, frames = Simulation.load(args.sims_path)
+    frames = list(frames)
+    # num = 100
+    # nums = [838. 348. 14. 664, 725, 974]
+    # frames.sort(key=lambda x: x["energy"])
+    # frame = frames[-32]
+    # for i, f in enumerate(frames):
+    #    if i == 858:
+    frame = frames[args.frame_num]
+    n, w, h = frame["domain"][0], frame["domain"][1], frame["domain"][2]
+    frame = sim.energy.mode(*frame["domain"], frame["arr"])
+
+    # Set up size and scaling variables.
+    dim = np.array([w, h])
+    vor = frame.vor_data
+
+    alpha = w / h
+    if torus:
+        r = alpha * size / 2
+        R = size / 2 - r
+        height_bounds = (0.5 * r, 1.25 * r)
+        scale = np.array([2 * np.pi, 2 * np.pi])
+    else:
+        height_bounds = (size / 10, 0.75 * size)
+        scale = np.array([alpha * size, size])
+
+    # Rescale and translate domain.
+    points = (vor.points - dim / 2) * scale / dim
+    vertices = (vor.vertices - dim / 2) * scale / dim
+
+    # Obtain height scaling.
+    heights = frame.stats["site_energies"] - ordered.e_hex(sim.domain)
+    heights -= np.min(heights)
+    heights *= (height_bounds[1] - height_bounds[0]) / np.max(heights)
+    heights += height_bounds[0]
+
+    # Prepare oriented regions and site_vert list
+    sites = points[:n]
+    regions = []
+    for i, x in enumerate(vor.point_region[:n]):
+        region = vor.regions[x]
+        x = sites[i]
+        p, q = vertices[region[0]] - x, vertices[region[1]] - x
+        regions.append(region if np.cross(p, q) >= 0 else region[::-1])
+
+    site_verts = []
+    for region in regions:
+        site_verts.append(vertices[region])
+
+    all_verts, all_faces = [], []
+    offsets = np.empty((n,), dtype=int)
+    for i, verts in enumerate(site_verts):
+        c = centroid(sites[i], verts)
+        if torus:
+            face_verts, faces, offset = torus_region(c, verts, heights[i])
+        else:
+            face_verts, faces, offset = flat_region(c, verts, heights[i])
+
+        old_len = len(all_verts)
+        if i == 0:
+            all_verts, all_faces = face_verts, faces
+        else:
+            all_verts = np.concatenate((all_verts, face_verts))
+            all_faces = np.concatenate((all_faces, faces + old_len))
+
+        offsets[i] = old_len + offset
+
+    # Merge regions
+    for i, x in enumerate(vor.ridge_points):
+        x, y = x[0], x[1]  # Site indices.
+        no_x, no_y = x >= sim.domain.n, y >= sim.domain.n
+        if no_x and no_y:
+            continue
+
+        adj_verts = vor.ridge_vertices[i]
+        v1, v2 = vertices[adj_verts[0]], vertices[adj_verts[1]]
+
+        bot_face, top_face = None, None
+        if torus:
+            sd = SUBDIVISION_AMOUNT - 1
+            if no_x:
+                x, y = y, x
+            y = y % n
+
+            v1_x_ind = regions[x].index(adj_verts[0])
+            v2_x_ind = regions[x].index(adj_verts[1])
+            # Need to find matching vertices
+            v1, v2 = v1 % (2 * np.pi), v2 % (2 * np.pi)
+            y_verts = vertices[regions[y]] % (2 * np.pi)
+            v1_y_ind = np.argmin(np.linalg.norm(y_verts - v1, axis=1))
+            v2_y_ind = np.argmin(np.linalg.norm(y_verts - v2, axis=1))
+
+            # We want x lower than y.
+            off_x, off_y = offsets[x], offsets[y]
+            if heights[x] > heights[y]:
+                v1_x_ind, v1_y_ind = v1_y_ind, v1_x_ind
+                v2_x_ind, v2_y_ind = v2_y_ind, v2_x_ind
+                off_x, off_y = off_y, off_x
+
+            if not_in_order(v1_x_ind, v2_x_ind):
+                v1_x_ind, v2_x_ind = v2_x_ind, v1_x_ind
+
+            if not_in_order(v1_y_ind, v2_y_ind):
+                v1_y_ind, v2_y_ind = v2_y_ind, v1_y_ind
+
+            sds = np.atleast_2d(np.arange(sd)).T
+            xm, ym = sd * v1_x_ind + sds, sd * v1_y_ind + sds
+            xp, yp = xm + 1, ym + 1
+            xp[-1], yp[-1] = sd * v2_x_ind, sd * v2_y_ind
+
+            xm += off_x
+            xp += off_x
+            ym += off_y
+            yp += off_y
+
+            bot_face = np.hstack((xm, ym[::-1], xp))
+            top_face = np.hstack((xm, yp[::-1], ym[::-1]))
+
+        else:
+            if no_x ^ no_y:
+                if no_x:
+                    x = y
+                m = len(regions[x])
+                v1_ind = regions[x].index(adj_verts[0])
+                v2_ind = regions[x].index(adj_verts[1])
+                if not_in_order(v1_ind, v2_ind):
+                    v1_ind, v2_ind = v2_ind, v1_ind
+
+                v1_ind += offsets[x]
+                v2_ind += offsets[x]
+                bot_face = np.array([[v1_ind, m + v1_ind, m + v2_ind]], dtype=int)
+                top_face = np.array([[v1_ind, m + v2_ind, v2_ind]], dtype=int)
+            else:
+                v1_x_ind = regions[x].index(adj_verts[0])
+                v2_x_ind = regions[x].index(adj_verts[1])
+                v1_y_ind = regions[y].index(adj_verts[0])
+                v2_y_ind = regions[y].index(adj_verts[1])
+
+                # We want x lower than y.
+                off_x, off_y = offsets[x], offsets[y]
+                if heights[x] > heights[y]:
+                    v1_x_ind, v1_y_ind = v1_y_ind, v1_x_ind
+                    v2_x_ind, v2_y_ind = v2_y_ind, v2_x_ind
+                    off_x, off_y = off_y, off_x
+
+                if not_in_order(v1_x_ind, v2_x_ind):
+                    v1_x_ind, v2_x_ind = v2_x_ind, v1_x_ind
+
+                if not_in_order(v1_y_ind, v2_y_ind):
+                    v1_y_ind, v2_y_ind = v2_y_ind, v1_y_ind
+
+                v1_x_ind += off_x
+                v2_x_ind += off_x
+                v1_y_ind += off_y
+                v2_y_ind += off_y
+
+                bot_face = np.array([[v1_x_ind, v1_y_ind, v2_x_ind]], dtype=int)
+                top_face = np.array([[v1_x_ind, v2_y_ind, v1_y_ind]], dtype=int)
+
+        all_faces = np.concatenate((all_faces, bot_face, top_face))
+
+    if torus:
+        all_verts = torus_transform(all_verts, R, all_verts[:, 2])
+
+    # Output into mesh
+    surf = mesh.Mesh(np.zeros(all_faces.shape[0], dtype=mesh.Mesh.dtype))
+    for i, f in enumerate(all_faces):
+        for j in range(3):
+            surf.vectors[i][j] = all_verts[f[j], :]
+
+    m_type = "Torus" if torus else "Flat"
+    out = OUTPUT_DIR / f"N{n} - {m_type} - {args.frame_num}.stl"
+    # temp = f"n{n}{m_type.lower()}{num}.stl"
+    surf.save(out)
+    # os.rename(temp, out)
+    print(f"Wrote to {out}.")
+    # n, alpha = 20, np.sqrt(3) * 2 / 5
+    # u, v = 5, 2
+    # domain = DomainParams(n, alpha, 1, 1)
+    # size = 1
+
+    # R, r = size, size * (1 - alpha) / alpha
+
+    # sites = ordered.sites(domain, (u, v))
+
+
+if __name__ == "__main__":
+    np.seterr(divide="ignore")
+    main()
diff --git a/scripts/vee_diff_pod.py b/scripts/vee_diff_pod.py
index d32048c..985691b 100644
--- a/scripts/vee_diff_pod.py
+++ b/scripts/vee_diff_pod.py
@@ -43,7 +43,7 @@ def main():
 
     plt.rcParams.update(RC_SETTINGS)
 
-    fig = plt.figure(figsize=(15, 15))
+    fig = plt.figure(figsize=(20, 9))
     gs = fig.add_gridspec(1, 1)
     ax = fig.add_subplot(gs[0])
 
@@ -78,11 +78,11 @@ def main():
             100 * (min_order - min_disorder), all_disorder_count, label=f"N={domain.n}"
         )
 
-    ax.set_ylabel("Perecent of Disordered Equilibria")
+    ax.set_ylabel("POD")
     ax.set_xlabel(r"VEE Difference $\left[\times 10^{2}\right]$")
     ax.yaxis.set_major_formatter(mtick.PercentFormatter())
     ax.set_ylim(48, 102)
-    ax.set_yticks(np.arange(50, 101, 5))
+    ax.set_yticks(np.arange(50, 101, 10))
 
     ax.grid(zorder=0)
     ax.legend()
diff --git a/squish/diagram.py b/squish/diagram.py
index c320e6e..75218fa 100644
--- a/squish/diagram.py
+++ b/squish/diagram.py
@@ -9,7 +9,7 @@ from mpl_toolkits.axes_grid1.axes_divider import make_axes_locatable
 
 import numpy as np, os, math
 from matplotlib.ticker import MaxNLocator, FormatStrFormatter
-from scipy.spatial import Voronoi, voronoi_plot_2d
+from scipy.spatial import Voronoi
 from multiprocessing import Pool, cpu_count
 
 from .common import DomainParams, OUTPUT_DIR
@@ -125,37 +125,13 @@ class Diagram:
         self.diagrams = np.atleast_2d(diagrams)
         self.cumulative = cumulative
 
-    def generate_frame(self, frame: int, mode: str, fol: str, name: str = None) -> None:
+    def generate_frame(
+        self, fig: Figure, frame: int, mode: str, fol: str, name: str = None
+    ) -> None:
         if mode not in ["save", "open"]:
             raise ValueError("Not a valid mode for diagrams!")
 
         shape = self.diagrams.shape
-
-        plt.rcParams.update(
-            {
-                "axes.titlesize": 45,
-                "axes.labelsize": 45,
-                "xtick.labelsize": 40,
-                "ytick.labelsize": 40,
-                "xtick.major.width": 2,
-                "ytick.major.width": 2,
-                "xtick.major.size": 5,
-                "ytick.major.size": 5,
-                "xtick.minor.width": 1,
-                "ytick.minor.width": 1,
-                "xtick.minor.size": 3,
-                "ytick.minor.size": 3,
-                "legend.fontsize": 40,
-                "lines.linewidth": 3,
-                "font.family": "cm",
-                "font.size": 40,
-                "text.usetex": True,
-                "text.latex.preamble": r"\usepackage{amsmath}",
-                "figure.constrained_layout.use": True,
-            }
-        )
-
-        fig = plt.figure(figsize=(shape[1] * 15, shape[0] * 15))
         gs = fig.add_gridspec(shape[0], shape[1])
         # fig, axes = plt.subplots(*shape, figsize=(shape[1] * 15, shape[0] * 15))
 
@@ -178,29 +154,28 @@ class Diagram:
 
         if mode == "save":
             plt.savefig(self.sim.path / fol / name)
-            plt.close(fig)
+            fig.clear()
         elif mode == "show":
             plt.show()
 
     def voronoi_plot(self, i: int, ax: AxesSubplot) -> None:
         domain = self.sim.domains[i]
         n, w, h = domain.n, domain.w, domain.h
+        flip = w <= h
         scale = 1.5
-        area = n <= 50
 
         e_hex = ordered.e_hex(domain)
+        line_color, point_color = "white", "lightseagreen"
 
         # Make color map axis.
         divider = make_axes_locatable(ax)
         cax = divider.append_axes("right", size="5%", pad=0.05)
 
-        voronoi_plot_2d(
-            self.sim.voronois[i], ax, show_vertices=False, show_points=False
-        )
+        vor = self.sim.voronois[i]
 
         # Mark location axis with 3 ticks.
-        ax.set_xticks([0, w / 2, w])
-        ax.set_yticks([0, h / 2, h])
+        ax.set_xticks(np.round([0, w / 2, w], 2))
+        ax.set_yticks(np.round([0, h / 2, h], 2))
 
         # Obtain site energies, and make color map.
         energies = self.sim.stats[i]["site_energies"][:n]
@@ -208,74 +183,77 @@ class Diagram:
         mapper = cm.ScalarMappable(norm=norm, cmap=cm.magma)
         cbar = plt.colorbar(mapper, cax=cax, ticks=np.linspace(-2, 2, 5))
 
+        SYMM = []
+        if flip:
+            dup = int((h // w + 1) // 2 + 1)
+            for i in range(-dup, dup + 1):
+                SYMM += [[i, -1], [i, 0], [i, 1]]
+        else:
+            dup = int((w // h + 1) // 2 + 1)
+            for i in range(-dup, dup + 1):
+                SYMM += [[-1, i], [0, i], [1, i]]
+
+        SYMM = np.array(SYMM) * np.array([w, h])
+
+        regions = [vor.regions[x] for x in vor.point_region[:n]]
+        edges = np.array([len(r) for r in regions])
+
         # Fill polygons with energy colormap.
-        for j, p in enumerate(self.sim.voronois[i].point_region):
-            region = self.sim.voronois[i].regions[p]
-            if not -1 in region:
-                polygon = [self.sim.voronois[i].vertices[k] for k in region]
-                ax.fill(
-                    *zip(*polygon),
-                    color=mapper.to_rgba(energies[j % n] - e_hex),
-                    zorder=0,
+        sizes = [40 - n / 100, 200 - n / 100, 40 - n / 100, 200 - n / 100, 40 - n / 100]
+        markers = [".", "p", ".", "*", "."]
+        for j in range(5):
+            sites = vor.points[np.argwhere(edges == j + 4)]
+            if len(sites) == 0:
+                continue
+            for s in SYMM:
+                ax.scatter(
+                    *(sites + s).T,
+                    marker=markers[j],
+                    s=sizes[j],
+                    color=point_color,
+                    zorder=1,
                 )
 
-        line_color, point_color = "white", "lightseagreen"
+        polygons = []
+        for region in regions:
+            polygon = np.empty((2 * len(SYMM), len(region)), dtype=float)
+            for i, s in enumerate(SYMM):
+                polygon[i * 2 : (i + 1) * 2, :] = (vor.vertices[region] + s).T
+            polygons.append(polygon)
+
+        colors = [mapper.to_rgba(energies[i] - e_hex) for i in range(n)]
+        for col, poly_set in zip(colors, polygons):
+            ax.fill(*poly_set, color=col, edgecolor="black", zorder=0)
+
+        ax.axis("equal")
 
         # Periodic border
-        ax.plot([-w, 2 * w], [0, 0], line_color, linewidth="4")
-        ax.plot([-w, 2 * w], [h, h], line_color, linewidth="4")
-        ax.plot([0, 0], [-h, 2 * h], line_color, linewidth="4")
-        ax.plot([w, w], [-h, 2 * h], line_color, linewidth="4")
-        ax.axis("equal")
-        ax.set_xlim([(1 - 0.85 * scale) * w / 2, (1 + 0.85 * scale) * w / 2])
-        ax.set_ylim([(1 - scale) * h / 2, (1 + scale) * h / 2])
+        if flip:
+            bot, top = (1 - scale) * h / 2, (1 + scale) * h / 2
+            left, right = w / 2 - 0.88 * (top - bot) / 2, w / 2 + 0.88 * (top - bot) / 2
+            ax.plot([left, right], [0, 0], line_color, linewidth="4")
+            ax.plot([left, right], [h, h], line_color, linewidth="4")
+            for i in range(-dup, dup + 2):
+                ax.plot([i * w, i * w], [-h, 2 * h], line_color, linewidth="4")
+        else:
+            left, right = (1 - 0.85 * scale) * w / 2, (1 + 0.85 * scale) * w / 2
+            bot, top = (h / 2 - (scale * w) / 2, h / 2 + (scale * w) / 2)
+
+            ax.plot([0, 0], [bot, top], line_color, linewidth="4")
+            ax.plot([w, w], [bot, top], line_color, linewidth="4")
+            for i in range(-dup, dup + 2):
+                ax.plot([-w, 2 * w], [i * h, i * h], line_color, linewidth="4")
+
+        ax.set_xlim(left, right)
+        ax.set_ylim(bot, top)
         ax.set_title("Voronoi Tessellation")
 
-        # Marking sites and defects.
-        defects = {5: {"x": [], "y": []}, 7: {"x": [], "y": []}}
-        for j in range(n):
-            for s in SYMM:
-                vec = self.sim.voronois[i].points[j] + s * self.sim.domains[i].dim
-                if area:
-                    txt = ax.text(
-                        *vec, str(round(self.sim.stats[i]["site_areas"][j], 3))
-                    )
-                    txt.set_clip_on(True)
-
-                if self.sim.stats[i]["site_edge_count"][j] == 5:
-                    defects[5]["x"].append(vec[0])
-                    defects[5]["y"].append(vec[1])
-                elif self.sim.stats[i]["site_edge_count"][j] == 7:
-                    defects[7]["x"].append(vec[0])
-                    defects[7]["y"].append(vec[1])
-                else:
-                    ax.scatter(
-                        vec[0], vec[1], s=(40 - n / 100), color=point_color, zorder=1
-                    )
-
-        ax.scatter(
-            defects[5]["x"],
-            defects[5]["y"],
-            marker="p",
-            s=(200 - n / 100),
-            color=point_color,
-            zorder=1,
-        )
-        ax.scatter(
-            defects[7]["x"],
-            defects[7]["y"],
-            marker="*",
-            s=(200 - n / 100),
-            color=point_color,
-            zorder=1,
-        )
-
         # Make VEE text in top left.
         props = dict(boxstyle="round", facecolor="white", alpha=0.8, zorder=20)
         ax.text(
             0.065,
             0.935,
-            f"VEE: {round(sum(energies)/n - e_hex, 8)}",
+            f"VEE: {sum(energies)/n - e_hex: .8f}",
             transform=ax.transAxes,
             verticalalignment="top",
             bbox=props,
@@ -283,16 +261,28 @@ class Diagram:
 
     def energy_plot(self, i: int, ax: AxesSubplot) -> None:
         ax.set_xlim([0, len(self.sim)])
+        e_hex = ordered.e_hex(self.sim.domains[i])
 
-        energies = self.sim.energies[: (i + 1)]
-        ax.plot(list(range(i + 1)), energies)
-        ax.title.set_text("Energy vs. Time")
+        vees = np.array(self.sim.energies) / self.sim.domains[i].n - e_hex
+
+        ax.plot(list(range(len(vees))), vees)
+        ax.scatter(
+            i,
+            vees[i],
+            s=250,
+            facecolors="none",
+            edgecolors="C6",
+            linewidth=4,
+            zorder=100,
+        )
+
+        # ax.title.set_text("VEE")
         # max_value = round(self.sim[0].energy)
         # min_value = round(self.sim[-1].energy)
         # diff = max_value-min_value
         # ax.set_yticks(np.arange(int(min_value-diff/5), int(max_value+diff/5), diff/25))
-        ax.set_xlabel("Iterations")
-        ax.set_ylabel("Energy")
+        ax.set_xlabel("Frame")
+        ax.set_ylabel("VEE")
         ax.grid()
 
     def site_areas_plot(self, i: int, ax: AxesSubplot) -> None:
@@ -424,11 +414,8 @@ class Diagram:
         self.sim.path.mkdir(exist_ok=True)
         (self.sim.path / fol).mkdir(exist_ok=True)
         combo_list = []
-        for i in range(cpu_count()):
-            start, end = (
-                int(i * len(frames) / cpu_count()),
-                int((i + 1) * len(frames) / cpu_count()),
-            )
+        for i in range(1):
+            start, end = (int(i * len(frames) / 1), int((i + 1) * len(frames) / 1))
             new_dia = Diagram(None, self.diagrams, self.cumulative)
             new_dia.sim = self.sim.slice(frames[start:end])
             combo_list.append(
@@ -448,10 +435,10 @@ class Diagram:
         if mode not in ["use_all", "sample"]:
             raise ValueError("Not a valid mode for videos!")
 
+        fps = 30
         if mode == "use_all":
             frames = list(range(len(self.sim)))
         elif mode == "sample":
-            fps = 30
             if len(self.sim) < fps * time:
                 frames = list(range(len(self.sim)))
                 fps = len(self.sim) / time
@@ -481,8 +468,36 @@ class Diagram:
 
 def render_frame_range(combo: Tuple[Diagram, str, int, int, int]) -> None:
     self, fol, offset, length, num_frames = combo
+
+    plt.rcParams.update(
+        {
+            "axes.titlesize": 45,
+            "axes.labelsize": 45,
+            "xtick.labelsize": 40,
+            "ytick.labelsize": 40,
+            "xtick.major.width": 2,
+            "ytick.major.width": 2,
+            "xtick.major.size": 5,
+            "ytick.major.size": 5,
+            "xtick.minor.width": 1,
+            "ytick.minor.width": 1,
+            "xtick.minor.size": 3,
+            "ytick.minor.size": 3,
+            "legend.fontsize": 40,
+            "lines.linewidth": 3,
+            "font.family": "cm",
+            "font.size": 40,
+            "text.usetex": True,
+            "text.latex.preamble": r"\usepackage{amsmath}",
+            "figure.constrained_layout.use": True,
+        }
+    )
+
+    # Generate figure here and pass in to prevent matplotlib memory leak.
+    shape = self.diagrams.shape
+    fig = plt.figure(figsize=(shape[1] * 15, shape[0] * 15))
     for i in range(length):
-        self.generate_frame(i, "save", fol, f"img{i+offset:05}.png")
+        self.generate_frame(fig, i, "save", fol, f"img{i+offset:05}.png")
         i = len(list((self.sim.path / fol).iterdir()))
         hashes = int(21 * i / num_frames)
         print(
@@ -491,3 +506,4 @@ def render_frame_range(combo: Tuple[Diagram, str, int, int, int]) -> None:
             flush=True,
             end="\r",
         )
+    plt.close(fig)
diff --git a/squish/ordered.py b/squish/ordered.py
index 52f6900..5ed5fc6 100644
--- a/squish/ordered.py
+++ b/squish/ordered.py
@@ -14,6 +14,13 @@ def e_hex(domain: DomainParams) -> float:
     return 2 - 2 * domain.r * (2 * pi) * R_HEX + 2 * pi * domain.r ** 2
 
 
+def toroidal_distance(domain: DomainParams, x: np.ndarray, y: np.ndarray) -> float:
+    dim = np.array([domain.w, domain.h])
+    absdist = np.abs(y - x)
+    wrap = dim * (absdist >= dim / 2)
+    return np.linalg.norm(wrap - absdist, axis=1)
+
+
 def configurations(domain: DomainParams) -> List[Config]:
     n = domain.n
     coprimes, valid = [], []
@@ -41,30 +48,37 @@ def configurations(domain: DomainParams) -> List[Config]:
 def get_config_generators(
     domain: DomainParams, config: Config
 ) -> Tuple[Config, Config]:
+    # x = sites(domain, config)
+    # x = x[x[:, 1] <= (domain.h / 2)]
+    # mag_x = toroidal_distance(domain, x, np.zeros(x.shape))
+    # x = x[np.argsort(mag_x)]
+    # print(x)
+
+    # return [tuple(x[1]), tuple(x[2])]
+
     n, w, h = domain.n, domain.w, domain.h
-    q1 = sites(domain, config)[1:]
-    all_sites = np.concatenate(
-        (q1, q1 - np.array([w, 0]), q1 - np.array([0, h]), q1 - domain.dim)
+    x = sites(domain, config)[1:]  # Remove 0
+    # Concatenate quadrants 1, 2 and 4.
+    x = np.concatenate(
+        (x, x - np.array([w, 0]), x - np.array([0, h]), x - np.array([w, h]))
     )
-    # Sort sites by magnitude and smallest.
-    all_sites = sorted(list(all_sites), key=lambda x: np.linalg.norm(x))
-    v = all_sites[0]  # Smallest vector set to v.
 
-    all_sites = np.array(all_sites[1:])  # Remove v from search set.
+    # Reduce search space
+    # x = x[(np.abs(x[:, 0]) <= (domain.w / 2)) * (np.abs(x[:, 1]) <= (domain.h / 2))]
+    mag_x = np.linalg.norm(x, axis=1)
+    x = x[np.argsort(mag_x)]  # Sort by magnitude
+
+    v = x[0]
+    x = x[1:]  # Remove v from search set.
     # Checking 0 < ax + by < v*v to make the sites are within the region.
-    tol = 1e-3
-    vdot = np.matmul(all_sites, v)
-    in_box = all_sites[np.where((-tol <= vdot) & (vdot <= (v.dot(v) + tol)))[0]]
-    in_box = np.expand_dims(in_box, 0).swapaxes(
-        0, 1
-    )  # Used for the next step, getting site*site
+    tol = 1e-8
+    vdot = np.matmul(x, v)
+    x = x[(-tol <= vdot) * (vdot <= v.dot(v) + tol)]
+    mag_x = np.linalg.norm(x, axis=1)
+    x = x[np.argsort(mag_x)]
 
-    v2 = in_box[
-        np.argmin(np.squeeze(np.matmul(in_box, in_box.transpose(0, 2, 1))))
-    ].flatten()
+    v2 = x[0]
 
-    if np.all(v == v2):
-        print(v, v2, n, w, h, config)
     return tuple(v), tuple(v2)
 
 
@@ -105,7 +119,8 @@ def avg_rp(d: float, l: float) -> float:
 def circumcenter(v: numpy.ndarray, w: numpy.ndarray) -> numpy.ndarray:
     det = 1 / (2 * rot(v).dot(w))
     if rot(v).dot(w) == 0:
-        print(v, w)
+        pass
+        # print(v, w)
     v2, w2 = v.dot(v), w.dot(w)
     c = np.empty((2,))
     c[0], c[1] = w[1] * v2 - v[1] * w2, -w[0] * v2 + v[0] * w2
diff --git a/squish/simulation.py b/squish/simulation.py
index c01eb30..8d706bf 100644
--- a/squish/simulation.py
+++ b/squish/simulation.py
@@ -217,7 +217,6 @@ class Flow(Simulation):
         log: bool,
         log_steps: int,
         new_sites: Optional[numpy.ndarray] = None,
-        newton: bool = False,
     ) -> None:
         if log:
             print(f"Find - {self.domain}", flush=True)
@@ -242,11 +241,11 @@ class Flow(Simulation):
             new_frame = self.energy.mode(*self.domain, frame.add_sites(change))
             end = timer()
 
-            grad_norm = np.linalg.norm(grad)
+            grad_norm = np.linalg.norm(grad) / (self.domain.n ** 0.5)
 
             if self.adaptive:
                 error = change - grad * self.step_size
-                tol = 10 ** min(-5, 2.3 * log10(grad_norm))
+                tol = 10 ** min(-3, -2 + log10(grad_norm))
                 # tol = 10 ** -10
                 self.step_size *= (tol / np.linalg.norm(error)) ** 0.5
 
@@ -393,10 +392,7 @@ class Shrink(Simulation):
     ) -> None:
         super().__init__(domain, energy, name=name)
         self.step_size, self.thres, self.adaptive = step_size, thres, adaptive
-        self.delta, self.stop_width = (
-            self.domain.w * delta / 100,
-            self.domain.w * stop_width,
-        )
+        self.delta, self.stop_width = delta, stop_width
 
     @property
     def initial_data(self) -> Dict:
@@ -425,8 +421,19 @@ class Shrink(Simulation):
             self.save(self.initial_data, True)
 
         width = self.domain.w if len(self.frames) == 0 else self.frames[-1].w
+        if len(self.frames) != 0:
+            new_sites = self.frames[-1].site_arr
+            width = self.frames[-1].w - self.delta
+        else:
+            width = self.domain.w
         i = 0
-        while width >= self.stop_width:
+
+        if delta < 0:
+            cond = width <= self.stop_width
+        if delta > 0:
+            cond = width <= self.stop_width
+
+        while cond:
             # Get to equilibrium.
             new_domain = DomainParams(
                 self.domain.n, width, self.domain.h, self.domain.r