diff --git a/badges/test_coverage.svg b/badges/test_coverage.svg
index c97e50489fae4a7075aff56a73e3e1936d1f76f1..dfb99a19bbb554f31bc566f509dce7df5d8bb40d 100644
--- a/badges/test_coverage.svg
+++ b/badges/test_coverage.svg
@@ -15,7 +15,7 @@
     <g fill="#fff" text-anchor="middle" font-family="DejaVu Sans,Verdana,Geneva,sans-serif" font-size="11">
         <text x="31.5" y="15" fill="#010101" fill-opacity=".3">coverage</text>
         <text x="31.5" y="14">coverage</text>
-        <text x="80" y="15" fill="#010101" fill-opacity=".3">66%</text>
-        <text x="80" y="14">66%</text>
+        <text x="80" y="15" fill="#010101" fill-opacity=".3">67%</text>
+        <text x="80" y="14">67%</text>
     </g>
 </svg>
diff --git a/binarycpython/tests/test_custom_logging.py b/binarycpython/tests/test_custom_logging.py
index 329dd7da59fc57aa524d74de9683f32f561d4f7c..55c618c67c58cfa4079c74d1d36355014772ebe3 100644
--- a/binarycpython/tests/test_custom_logging.py
+++ b/binarycpython/tests/test_custom_logging.py
@@ -65,20 +65,21 @@ class test_binary_c_log_code(unittest.TestCase):
     def test_binary_c_log_code(self):
         with Capturing() as _:
             self._test_binary_c_log_code()
-        # print("\n".join(output))
 
     def _test_binary_c_log_code(self):
         """
         Test to see if passing a print statement to the function results in correct binary_c output
         """
 
+        # binary_c should not compile anything here
         input_1 = "None"
         output_1 = binary_c_log_code(input_1, verbosity=1)
         self.assertEqual(output_1, None, msg="Output should be None")
 
+        #
         input_2 = 'Printf("MY_STELLAR_DATA %g %g %g %g\\n",((double)stardata->model.time),((double)stardata->star[0].mass),((double)stardata->model.probability),((double)stardata->model.dt));'
         output_2 = binary_c_log_code(input_2, verbosity=1)
-        test_value_2 = '#pragma push_macro("Max")\n#pragma push_macro("Min")\n#undef Max\n#undef Min\n#include "binary_c.h"\n\n// add visibility __attribute__ ((visibility ("default"))) to it \nvoid binary_c_API_function custom_output_function(struct stardata_t * stardata);\nvoid binary_c_API_function custom_output_function(struct stardata_t * stardata)\n{\n    // struct stardata_t * stardata = (struct stardata_t *)x;\n    Printf("MY_STELLAR_DATA %g %g %g %g\\n",((double)stardata->model.time),((double)stardata->star[0].mass),((double)stardata->model.probability),((double)stardata->model.dt));;\n}\n\n#undef Max \n#undef Min\n#pragma pop_macro("Min")\n#pragma pop_macro("Max")    '
+        test_value_2 = '#pragma push_macro("Max")\n#pragma push_macro("Min")\n#undef Max\n#undef Min\n#include "binary_c.h"\n\n// add visibility __attribute__ ((visibility ("default"))) to it\nvoid binary_c_API_function custom_output_function(struct stardata_t * stardata);\nvoid binary_c_API_function custom_output_function(struct stardata_t * stardata)\n{\n    // struct stardata_t * stardata = (struct stardata_t *)x;\n    Printf("MY_STELLAR_DATA %g %g %g %g\\n",((double)stardata->model.time),((double)stardata->star[0].mass),((double)stardata->model.probability),((double)stardata->model.dt));;\n}\n\n#undef Max\n#undef Min\n#pragma pop_macro("Min")\n#pragma pop_macro("Max")    '
 
         self.assertEqual(
             output_2,
@@ -153,7 +154,7 @@ class test_from_binary_c_config(unittest.TestCase):
         output_2 = from_binary_c_config(BINARY_C_CONFIG, input_2)
         self.assertIn(
             output_2,
-            ["2.1.7", "2.2pre1", "2.2.0", "2.2.1"],
+            ["2.1.7", "2.2pre1", "2.2.0", "2.2.1", "2.2.2"],
             msg="binary_c version doesnt match",
         )
 
@@ -213,4 +214,5 @@ class test_create_and_load_logging_function(unittest.TestCase):
 
 
 if __name__ == "__main__":
-    unittest.main()
+    # unittest.main()
+    test_binary_c_log_code().test_binary_c_log_code()
diff --git a/binarycpython/tests/test_grid.py b/binarycpython/tests/test_grid.py
index 97afaec8c9032458be7f0d2eff2f9a71ceb6a40c..413f8189cfe092ac57ce41c0ae729989c7c4888e 100644
--- a/binarycpython/tests/test_grid.py
+++ b/binarycpython/tests/test_grid.py
@@ -110,11 +110,11 @@ class test__setup(unittest.TestCase):
     Unittests for _setup function
     """
 
-    def test_setup(self):
-        with Capturing() as _:
-            self._test_setup()
+    # def test_setup(self):
+    #     with Capturing() as _:
+    #         self._test_setup()
 
-    def _test_setup(self):
+    def test_setup(self):
         """
         Unittests for function _setup
         """
@@ -124,7 +124,7 @@ class test__setup(unittest.TestCase):
         self.assertTrue("metallicity" in test_pop.defaults)
         self.assertNotIn("help_all", test_pop.cleaned_up_defaults)
         self.assertEqual(test_pop.bse_options, {})
-        self.assertEqual(test_pop.custom_options, {})
+        self.assertEqual(test_pop.custom_options, {"save_snapshot": False})
         self.assertEqual(test_pop.argline_dict, {})
         self.assertEqual(test_pop.persistent_data_memory_dict, {})
         self.assertTrue(test_pop.grid_options["parse_function"] is None)
@@ -617,6 +617,7 @@ class test_resultdict(unittest.TestCase):
             dphasevol="dlnm1",
             parameter_name="M_1",
             condition="",  # Impose a condition on this grid variable. Mostly for a check for yourself
+            gridtype="left",
         )
 
         ## Executing a population
@@ -775,7 +776,7 @@ class test_grid_evolve(unittest.TestCase):
             parse_function=parse_function_test_grid_evolve_2_threads_with_custom_logging,
         )
         test_pop.set(ensemble=0)
-        resolution = {"M_1": 2}
+        resolution = {"M_1": 20}
 
         test_pop.add_grid_variable(
             name="lnm1",
diff --git a/binarycpython/utils/population_extensions/distribution_functions.py b/binarycpython/utils/population_extensions/distribution_functions.py
index 67535d8424d75f74d5636c40d4f54eec5d7e80b8..5dedd62a81e15504855b3bdda503f3f348c1aded 100644
--- a/binarycpython/utils/population_extensions/distribution_functions.py
+++ b/binarycpython/utils/population_extensions/distribution_functions.py
@@ -302,7 +302,9 @@ class distribution_functions:
             m0, m1, m2, m_max, p1, p2, p3
         )
 
-        if m < m0:
+        if (
+            m <= m0
+        ):  # TODO: this needs to be checked! Compare two choices in this with the same sampling ranges.
             prob = 0.0  # Below lower bound
         elif m <= m1:
             prob = three_part_powerlaw_constants[0] * (m**p1)  # Between m0 and m1
@@ -759,7 +761,7 @@ class distribution_functions:
         log10Pmin = max(-1.0, log10Pmin)
 
         # save mass input and limit mass used (M1 from now on) to fitted range
-        Mwas = M1
+        # Mwas = M1 TODO: why is this here?
         M1 = max(1.15, min(16.3, M1))
 
         # Calculate the normalisations
@@ -921,7 +923,7 @@ class distribution_functions:
             The multiplicity fractions of the quadruples are folded onto that of the triples
         """
 
-        if not max_multiplicity in range(1, 5):
+        if max_multiplicity not in range(1, 5):
             msg = "\tMoe and di Stefano 2017: merge_multiplicities: max_multiplicity has to be between 1 and 4. It is {} now".format(
                 max_multiplicity
             )
@@ -1422,24 +1424,26 @@ class distribution_functions:
                 )
 
                 # now we integrate and renormalise (if the table is not all zero)
-                I = self.get_integration_constant_q(
+                integration_constant_q = self.get_integration_constant_q(
                     q_interpolator, tmp_table, qdata, verbosity=verbosity
                 )
 
-                if I > 0:
+                if integration_constant_q > 0:
                     # normalise to 1.0 by dividing the data by 1.0/$I
-                    q_interpolator.multiply_table_column(1, 1.0 / I)
+                    q_interpolator.multiply_table_column(
+                        1, 1.0 / integration_constant_q
+                    )
 
                     # test this
-                    new_I = self.get_integration_constant_q(
+                    new_integration_constant_q = self.get_integration_constant_q(
                         q_interpolator, tmp_table, qdata, verbosity=verbosity
                     )
 
                     # fail if error in integral > 1e-6 (should be ~ machine precision)
-                    if abs(1.0 - new_I) > 1e-6:
+                    if abs(1.0 - new_integration_constant_q) > 1e-6:
                         verbose_print(
                             "\tMoe and di Stefano 2017: build_q_table: Error: > 1e-6 in q probability integral: {}".format(
-                                I
+                                integration_constant_q
                             ),
                             verbosity,
                             _MOE2017_VERBOSITY_LEVEL,
@@ -1520,7 +1524,7 @@ class distribution_functions:
         """
 
         dq = 1e-3  # resolution of the integration/renormalisation
-        I = 0
+        integration_constant_q = 0
 
         # integrate: note that the value of the integral is
         # meaningless to within a factor (which depends on $dq)
@@ -1538,8 +1542,8 @@ class distribution_functions:
                 )
                 raise ValueError(msg)
 
-            I += x[0] * dq
-        return I
+            integration_constant_q += x[0] * dq
+        return integration_constant_q
 
     def fill_data(self, sample_values, data_dict):
         """
@@ -1550,19 +1554,19 @@ class distribution_functions:
         """
 
         data = {}
-        I = 0
+        integration_constant = 0
 
-        dstep = float(sample_values[1]) - float(sample_values[0])
+        # dstep = float(sample_values[1]) - float(sample_values[0])
 
         # Read out the data
         for sample_value in sample_values:
             val = data_dict[sample_value]
             data[sample_value] = val
-            I += val
+            integration_constant += val
 
         # Normalise the data
         for sample_value in sample_values:
-            data[sample_value] = data[sample_value] / I
+            data[sample_value] = data[sample_value] / integration_constant
 
         return data
 
@@ -1579,6 +1583,8 @@ class distribution_functions:
         Function to calculate the e integral
 
         We need to renormalise this because min_per > 0, and not all periods should be included
+
+        TODO: create function ot actually do this in a more general way
         """
 
         global Moecache
@@ -1595,7 +1601,7 @@ class distribution_functions:
 
         # Check for cached value. If it doesn't exist: calculate
         if not Moecache[integrals_string].get(mass_period_string, None):
-            I = 0
+            integation_constant_e = 0
             decc = 1e-3
 
             for ecc in np.arange(min_ecc, max_ecc, decc):
@@ -1608,10 +1614,10 @@ class distribution_functions:
                     ]
                 )[0]
 
-                I += dp_decc * decc
+                integation_constant_e += dp_decc * decc
 
             # Set the integral value in the dict
-            Moecache[integrals_string][mass_period_string] = I
+            Moecache[integrals_string][mass_period_string] = integation_constant_e
             verbose_print(
                 "\tMoe and di Stefano 2017: calc_ecc_integral: min_ecc: {} max ecc: {} integrals_string: {} interpolator_name: {} mass_string: {} period_string: {} mass: {} period: {} I: {}".format(
                     min_ecc,
@@ -1622,7 +1628,7 @@ class distribution_functions:
                     period_string,
                     options[mass_string],
                     options[period_string],
-                    I,
+                    integation_constant_e,
                 ),
                 verbosity,
                 _MOE2017_VERBOSITY_LEVEL,
@@ -1668,7 +1674,7 @@ class distribution_functions:
 
         # Check for cached value. If it doesn't exist: calculate
         if not Moecache[integrals_string].get(options[mass_string], None):
-            I = 0
+            integation_constant_P = 0
             dlogP = 1e-3
 
             for logP in np.arange(min_logP, max_logP, dlogP):
@@ -1677,10 +1683,10 @@ class distribution_functions:
                     [np.log10(options[mass_string]), logP]
                 )[0]
 
-                I += dp_dlogP * dlogP
+                integation_constant_P += dp_dlogP * dlogP
 
             # Set the integral value in the dict
-            Moecache[integrals_string][options[mass_string]] = I
+            Moecache[integrals_string][options[mass_string]] = integation_constant_P
             verbose_print(
                 "\tMoe and di Stefano 2017: calc_P_integral: min_logP: {} integrals_string: {} interpolator_name: {} mass_string: {} mass: {} I: {}".format(
                     min_logP,
@@ -1688,7 +1694,7 @@ class distribution_functions:
                     interpolator_name,
                     mass_string,
                     options[mass_string],
-                    I,
+                    integation_constant_P,
                 ),
                 verbosity,
                 _MOE2017_VERBOSITY_LEVEL,
@@ -2142,11 +2148,11 @@ class distribution_functions:
                         # as for any other stars but Pmax must be such that
                         # sep3 < sep2 * 0.2
 
-                        # TODO: fix this here
-                        max_sep3 = 0.2 * options["sep2"] * (1.0 + options["ecc2"])
-                        max_per3 = calc_period_from_sep(
-                            options["M_1+M_2"], options["mmin"], max_sep3
-                        )
+                        # # TODO: fix this here
+                        # max_sep3 = 0.2 * options["sep2"] * (1.0 + options["ecc2"])
+                        # max_per3 = calc_period_from_sep(
+                        #     options["M_1+M_2"], options["mmin"], max_sep3
+                        # )
 
                         # Calculate P integral or use the cached value
                         # TODO: Make sure we use the correct period idea here.
diff --git a/pyproject.toml b/pyproject.toml
index 44a9414147e85013fb81672a5f4f9f0d3e192e09..e879add9c5115d905adb90543bc7e00c8c956fbe 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,5 +1,4 @@
 [tool.black]
-
 target-version = ['py39']
 extend-exclude = '''
 # A regex preceded with ^/ will apply only to files and directories
@@ -20,7 +19,10 @@ extend-exclude = '''
 max-line-length = 89
 max-complexity = 18
 
-
-[tool.coverage]
-#omit =
-#    test_*.py
+#[tool.coverage]
+#[tool.coverage.run]
+#omit =[*/.local/*,
+## omit everything in /usr
+#/usr/*,
+## omit this single file
+#utils/tirefire.py]