diff --git a/sw/simulate_BER_curve.py b/sw/simulate_2d_BER.py
similarity index 67%
rename from sw/simulate_BER_curve.py
rename to sw/simulate_2d_BER.py
index 257705c..f2f9682 100644
--- a/sw/simulate_BER_curve.py
+++ b/sw/simulate_2d_BER.py
@@ -1,3 +1,5 @@
+import typing
+
 import numpy as np
 import pandas as pd
 import seaborn as sns
@@ -5,11 +7,12 @@ import matplotlib.pyplot as plt
 import signal
 from timeit import default_timer
 from tqdm import tqdm
+from dataclasses import dataclass
+from types import MappingProxyType
 
 from utility import codes, noise, misc
 from utility.simulation.simulators import GenericMultithreadedSimulator
 
-# from cpp_modules.cpp_decoders import ProximalDecoder
 from cpp_modules.cpp_decoders import ProximalDecoder_204_102 as ProximalDecoder
 
 
@@ -18,9 +21,18 @@ def count_bit_errors(d: np.array, d_hat: np.array) -> int:
 
 
 def task_func(params):
+    """Function called by the GenericMultithreadedSimulator instance.
+
+    Calculate the BER, FER, and DFR  for a given SNR and gamma.
+    """
     signal.signal(signal.SIGINT, signal.SIG_IGN)
 
-    decoder, max_iterations, SNR, n, k = params
+    decoder = params["decoder"]
+    max_iterations = params["max_iterations"]
+    SNR = params["SNR"]
+    n = params["n"]
+    k = params["k"]
+
     c = np.zeros(n)
     x_bpsk = c + 1
 
@@ -44,14 +56,15 @@ def task_func(params):
         if k_max == -1:
             dec_fails += 1
 
-        if total_frame_errors > 500:
+        if total_frame_errors > 100:
             break
 
     BER = total_bit_errors / (num_iterations * n)
     FER = total_frame_errors / num_iterations
     DFR = dec_fails / (num_iterations + dec_fails)
 
-    return BER, FER, DFR, num_iterations
+    return {"BER": BER, "FER": FER, "DFR": DFR,
+            "num_iterations": num_iterations}
 
 
 def simulate(H_file, SNRs, max_iterations, omega, K, gammas):
@@ -65,62 +78,45 @@ def simulate(H_file, SNRs, max_iterations, omega, K, gammas):
 
     # Define params different for each task
 
-    params = {}
+    task_params = []
     for i, SNR in enumerate(SNRs):
         for j, gamma in enumerate(gammas):
             decoder = ProximalDecoder(H=H.astype('int32'), K=K, omega=omega,
                                       gamma=gamma)
-            params[f"{i}_{j}"] = (decoder, max_iterations, SNR, n, k)
+
+            task_params.append(
+                {"decoder": decoder, "max_iterations": max_iterations,
+                 "SNR": SNR, "gamma": gamma, "n": n, "k": k})
 
     # Set up simulation
 
-    sim.task_params = params
+    sim.task_params = task_params
     sim.task_func = task_func
 
     sim.start_or_continue()
 
-    return sim.get_current_results()
-
-
-def reformat_data(results, SNRs, gammas):
-    data = {"BER": np.zeros(3 * 10), "FER": np.zeros(3 * 10),
-            "DFR": np.zeros(3 * 10), "gamma": np.zeros(3 * 10),
-            "SNR": np.zeros(3 * 10), "num_iter": np.zeros(3 * 10)}
-
-    for i, (key, (BER, FER, DFR, num_iter)) in enumerate(results.items()):
-        i_SNR, i_gamma = key.split('_')
-        data["BER"][i] = BER
-        data["FER"][i] = FER
-        data["DFR"][i] = DFR
-        data["num_iter"][i] = num_iter
-        data["SNR"][i] = SNRs[int(i_SNR)]
-        data["gamma"][i] = gammas[int(i_gamma)]
-
-    print(pd.DataFrame(data))
-    return pd.DataFrame(data)
+    return sim.current_results
 
 
 def main():
     # Set up simulation params
 
-    sim_name = "BER_FER_DFR"
+    sim_name = "2d_BER_FER_DFR"
 
+    # H_file = "BCH_7_4.alist"
+    # H_file = "BCH_31_11.alist"
+    # H_file = "BCH_31_26.alist"
     # H_file = "96.3.965.alist"
     H_file = "204.33.486.alist"
     # H_file = "204.33.484.alist"
     # H_file = "204.55.187.alist"
     # H_file = "408.33.844.alist"
-    # H_file = "BCH_7_4.alist"
-    # H_file = "BCH_31_11.alist"
-    # H_file = "BCH_31_26.alist"
-    SNRs = np.arange(1, 6, 0.5)
 
+    SNRs = np.arange(1, 6, 0.5)
     max_iterations = 20000
-    # omega = 0.005
-    # K = 60
     omega = 0.05
-    K = 60
-    gammas = [0.15, 0.01, 0.05]
+    K = 100
+    gammas = np.arange(0.0, 0.17, 0.01)
 
     # Run simulation
 
@@ -130,10 +126,13 @@ def main():
 
     print(f"duration: {end_time - start_time}")
 
-    df = reformat_data(results, SNRs, gammas)
+    df = misc.pgf_reformat_data_3d(results=results, x_param_name="SNR",
+                                   y_param_name="gamma",
+                                   z_param_names=["BER", "FER", "DFR",
+                                                  "num_iterations"])
 
-    df.to_csv(
-        f"sim_results/{sim_name}_{misc.slugify(H_file)}.csv")
+    # df.sort_values(by=["gamma", "SNR"]).to_csv(
+    #     f"sim_results/{sim_name}_{misc.slugify(H_file)}.csv", index=False)
 
     sns.set_theme()
     ax = sns.lineplot(data=df, x="SNR", y="BER", hue="gamma")