diff --git a/binarycpython/utils/custom_logging_functions.py b/binarycpython/utils/custom_logging_functions.py
index 0415f018ce21c318da3f0130a0d2cdf9f586f08b..e3be2e16ba91c5bfbcaae024c8fe48fdf0aa5156 100644
--- a/binarycpython/utils/custom_logging_functions.py
+++ b/binarycpython/utils/custom_logging_functions.py
@@ -132,13 +132,11 @@ def binary_c_log_code(code: str, verbose: int = 0) -> str:
 
     # Create code
     custom_logging_function_string = """\
-#pragma push_macro(\"MAX\")
-#pragma push_macro(\"MIN\")
-#undef MAX
-#undef MIN
-#include \"binary_c.h\"
-#include \"RLOF/RLOF_prototypes.h\"
-#include \"nucsyn/nucsyn_prototypes.h\"
+#pragma push_macro(\"Max\")
+#pragma push_macro(\"Min\")
+#undef Max
+#undef Min
+#include "binary_c.h"
 
 // add visibility __attribute__ ((visibility ("default"))) to it 
 void binary_c_API_function custom_output_function(struct stardata_t * stardata);
@@ -148,10 +146,10 @@ void binary_c_API_function custom_output_function(struct stardata_t * stardata)
     {};
 }}
 
-#undef MAX 
-#undef MIN
-#pragma pop_macro(\"MIN\")
-#pragma pop_macro(\"MAX\")\
+#undef Max 
+#undef Min
+#pragma pop_macro(\"Min\")
+#pragma pop_macro(\"Max\")\
     """.format(
         code
     )
@@ -257,7 +255,8 @@ def return_compilation_dict(verbose: int = 0) -> dict:
     bincincdirs = from_binary_c_config(BINARY_C_CONFIG, "incdirs")
 
     # combine
-    binclibs = " {} {} {}".format(libdirs, libbinary_c, binclibs)
+    # binclibs = " {} {} {}".format(libdirs, libbinary_c, binclibs)
+    binclibs = " {} {} {}".format(libdirs, binclibs, libbinary_c)
 
     # setup defaults:
     defaults = {
@@ -361,7 +360,7 @@ def compile_shared_lib(
             )
         )
     res = subprocess.check_output("{command}".format(command=command), shell=True)
-
+    print(res)
     if verbose > 0:
         if res:
             print("Output of compilation command:\n{}".format(res))
diff --git a/binarycpython/utils/grid.py b/binarycpython/utils/grid.py
index 9c3d25206973c06e58b4ddf33a7457ce7996d13d..b1ece9d95fef419350b007627708745a5c8183b6 100644
--- a/binarycpython/utils/grid.py
+++ b/binarycpython/utils/grid.py
@@ -356,6 +356,7 @@ class Population:
         dphasevol: Union[str, int],
         parameter_name: str,
         gridtype: str = "edge",
+        branchpoint: int = 0,
         precode: Union[str, None] = None,
         condition: Union[str, None] = None,
     ) -> None:
@@ -408,6 +409,7 @@ class Population:
         """
 
         # TODO: Add check for the gridtype input value
+        # TODO: add functionality for branchpoint
 
         # Add grid_variable
         grid_variable = {
@@ -422,6 +424,7 @@ class Population:
             "parameter_name": parameter_name,
             "condition": condition,
             "gridtype": gridtype,
+            "branchpoint": branchpoint,
             "grid_variable_number": len(self.grid_options["_grid_variables"]),
         }
 
@@ -1500,14 +1503,14 @@ class Population:
         code_string += indent * depth + "\n"
 
         # Setting up the system parameters
-        code_string += indent * depth + "m1 = None\n"
-        code_string += indent * depth + "m2 = None\n"
-        code_string += indent * depth + "m3 = None\n"
-        code_string += indent * depth + "m4 = None\n"
+        code_string += indent * depth + "M_1 = None\n"
+        code_string += indent * depth + "M_2 = None\n"
+        code_string += indent * depth + "M_3 = None\n"
+        code_string += indent * depth + "M_4 = None\n"
 
-        code_string += indent * depth + "per = None\n"
-        code_string += indent * depth + "per2 = None\n"
-        code_string += indent * depth + "per3 = None\n"
+        code_string += indent * depth + "orbital_period = None\n"
+        code_string += indent * depth + "orbital_period_triple = None\n"
+        code_string += indent * depth + "orbital_period_quadruple = None\n"
 
         code_string += indent * depth + "eccentricity = None\n"
         code_string += indent * depth + "eccentricity2 = None\n"
@@ -3048,34 +3051,34 @@ class Population:
 
         # # multiplicity_fractions = Moe_de_Stefano_2017_multiplicity_fractions(options)
 
-        import matplotlib.pyplot as plt
+        # import matplotlib.pyplot as plt
 
-        mass_range = range(1, 80)
+        # mass_range = range(1, 80)
 
-        multiplicity_dict = {}
-        for mass in mass_range:
-            options["M1"] = mass
-            multiplicity_fractions = Moe_de_Stefano_2017_multiplicity_fractions(options)
-            multiplicity_dict[mass] = multiplicity_fractions
+        # multiplicity_dict = {}
+        # for mass in mass_range:
+        #     options["M1"] = mass
+        #     multiplicity_fractions = Moe_de_Stefano_2017_multiplicity_fractions(options)
+        #     multiplicity_dict[mass] = multiplicity_fractions
 
-        single_values = [multiplicity_dict[key][0] for key in multiplicity_dict.keys()]
-        binary_values = [multiplicity_dict[key][1] for key in multiplicity_dict.keys()]
-        triple_values = [multiplicity_dict[key][2] for key in multiplicity_dict.keys()]
-        quad_values = [multiplicity_dict[key][3] for key in multiplicity_dict.keys()]
+        # single_values = [multiplicity_dict[key][0] for key in multiplicity_dict.keys()]
+        # binary_values = [multiplicity_dict[key][1] for key in multiplicity_dict.keys()]
+        # triple_values = [multiplicity_dict[key][2] for key in multiplicity_dict.keys()]
+        # quad_values = [multiplicity_dict[key][3] for key in multiplicity_dict.keys()]
 
-        print("mass = {}".format(list(mass_range)))
-        print("single_values={}".format(single_values))
-        print("binary_values={}".format(binary_values))
-        print("triple_values={}".format(triple_values))
-        print("quad_values={}".format(quad_values))
+        # print("mass = {}".format(list(mass_range)))
+        # print("single_values={}".format(single_values))
+        # print("binary_values={}".format(binary_values))
+        # print("triple_values={}".format(triple_values))
+        # print("quad_values={}".format(quad_values))
 
-        plt.plot(mass_range, single_values, label="single")
-        plt.plot(mass_range, binary_values, label="binary")
-        plt.plot(mass_range, triple_values, label="triple")
-        plt.plot(mass_range, quad_values, label="Quadruple")
-        plt.legend()
-        plt.xscale("log")
-        plt.show()
+        # plt.plot(mass_range, single_values, label="single")
+        # plt.plot(mass_range, binary_values, label="binary")
+        # plt.plot(mass_range, triple_values, label="triple")
+        # plt.plot(mass_range, quad_values, label="Quadruple")
+        # plt.legend()
+        # plt.xscale("log")
+        # plt.show()
 
         ############################################################
         # construct the grid here
@@ -3157,7 +3160,7 @@ class Population:
                 longname="log10(Orbital_Period)",
                 resolution=options["resolutions"]["logP"][0],
                 probdist=1.0,
-                condition='(self.grid_options["grid_options"]["multiplicity"] >= 2)',
+                condition='(self.grid_options["multiplicity"] >= 2)',
                 branchpoint=1,
                 gridtype="centred",
                 dphasevol="({} * dlog10per)".format(LOG_LN_CONVERTER),
@@ -3168,8 +3171,8 @@ class Population:
                     options["resolutions"]["logP"][0],
                 ),
                 precode="""orbital_period = 10.0**log10per
-    qmin={}/m1
-    qmax=self.maximum_mass_ratio_for_RLOF(m1, per)
+qmin={}/M_1
+qmax=self.maximum_mass_ratio_for_RLOF(M_1, orbital_period)
     """.format(
                     options.get("Mmin", 0.07)
                 ),
@@ -3183,7 +3186,7 @@ class Population:
                 valuerange=[
                     options["ranges"]["q"][0]
                     if options.get("ranges", {}).get("q", None)
-                    else "options.get('Mmin', 0.07)/m1",
+                    else "options.get('Mmin', 0.07)/M_1",
                     options["ranges"]["q"][1]
                     if options.get("ranges", {}).get("q", None)
                     else "qmax",
@@ -3193,13 +3196,13 @@ class Population:
                 gridtype="centred",
                 dphasevol="dq",
                 precode="""
-    M_2 = q * m1
-    sep = calc_sep_from_period(M_1, M_2, orbital_period)
+M_2 = q * M_1
+sep = calc_sep_from_period(M_1, M_2, orbital_period)
     """,
                 spacingfunc="const({}, {}, {})".format(
                     options["ranges"]["q"][0]
                     if options.get("ranges", {}).get("q", None)
-                    else "options.get('Mmin', 0.07)/m1",
+                    else "options.get('Mmin', 0.07)/M_1",
                     options["ranges"]["q"][1]
                     if options.get("ranges", {}).get("q", None)
                     else "qmax",
@@ -3238,7 +3241,7 @@ class Population:
                     longname="log10(Orbital_Period2)",
                     resolution=options["resolutions"]["logP"][1],
                     probdist=1.0,
-                    condition='(self.grid_options["grid_options"]["multiplicity"] >= 3)',
+                    condition='(self.grid_options["multiplicity"] >= 3)',
                     branchpoint=1,
                     gridtype="centred",
                     dphasevol="({} * dlog10per2)".format(LOG_LN_CONVERTER),
@@ -3252,8 +3255,8 @@ class Population:
                         options["resolutions"]["logP"][1],
                     ),
                     precode="""orbital_period_triple = 10.0**log10per2
-    q2min={}/(M_1+M_2)
-    q2max=self.maximum_mass_ratio_for_RLOF(M_1+M_2, orbital_period_triple)
+q2min={}/(M_1+M_2)
+q2max=self.maximum_mass_ratio_for_RLOF(M_1+M_2, orbital_period_triple)
     """.format(
                         options.get("Mmin", 0.07)
                     ),
@@ -3269,7 +3272,7 @@ class Population:
                     valuerange=[
                         options["ranges"]["q"][0]
                         if options.get("ranges", {}).get("q", None)
-                        else "options.get('Mmin', 0.07)/(m1+m2)",
+                        else "options.get('Mmin', 0.07)/(M_1+M_2)",
                         options["ranges"]["q"][1]
                         if options.get("ranges", {}).get("q", None)
                         else "q2max",
@@ -3279,14 +3282,14 @@ class Population:
                     gridtype="centred",
                     dphasevol="dq2",
                     precode="""
-    M_3 = q2 * (M_1 + M_2)
-    sep2 = calc_sep_from_period((M_1+M_2), M_3, orbital_period_triple)
-    eccentricity2=0
+M_3 = q2 * (M_1 + M_2)
+sep2 = calc_sep_from_period((M_1+M_2), M_3, orbital_period_triple)
+eccentricity2=0
     """,
                     spacingfunc="const({}, {}, {})".format(
                         options["ranges"]["q"][0]
                         if options.get("ranges", {}).get("q", None)
-                        else "options.get('Mmin', 0.07)/(m1+m2)",
+                        else "options.get('Mmin', 0.07)/(M_1+M_2)",
                         options["ranges"]["q"][1]
                         if options.get("ranges", {}).get("q", None)
                         else "q2max",
@@ -3324,7 +3327,7 @@ class Population:
                         longname="log10(Orbital_Period3)",
                         resolution=options["resolutions"]["logP"][2],
                         probdist=1.0,
-                        condition='(self.grid_options["grid_options"]["multiplicity"] >= 4)',
+                        condition='(self.grid_options["multiplicity"] >= 4)',
                         branchpoint=1,
                         gridtype="centred",
                         dphasevol="({} * dlog10per3)".format(LOG_LN_CONVERTER),
@@ -3338,8 +3341,8 @@ class Population:
                             options["resolutions"]["logP"][2],
                         ),
                         precode="""orbital_period_quadruple = 10.0**log10per3
-    q3min={}/(M_3)
-    q3max=self.maximum_mass_ratio_for_RLOF(M_3, orbital_period_quadruple)
+q3min={}/(M_3)
+q3max=self.maximum_mass_ratio_for_RLOF(M_3, orbital_period_quadruple)
     """.format(
                             options.get("Mmin", 0.07)
                         ),
@@ -3364,9 +3367,9 @@ class Population:
                         gridtype="centred",
                         dphasevol="dq3",
                         precode="""
-    M_4 = q3 * M_3
-    sep3 = calc_sep_from_period((M_3), M_4, orbital_period_quadruple)
-    eccentricity3=0
+M_4 = q3 * M_3
+sep3 = calc_sep_from_period((M_3), M_4, orbital_period_quadruple)
+eccentricity3=0
         """,
                         spacingfunc="const({}, {}, {})".format(
                             options["ranges"]["q"][0]