diff --git a/binarycpython/utils/grid.py b/binarycpython/utils/grid.py
index 73e73a909373dab6cc42ce2ff8379abee07a9bf6..1cc2942c82b3f45fe845dc384998df2cbd87f9cc 100644
--- a/binarycpython/utils/grid.py
+++ b/binarycpython/utils/grid.py
@@ -129,6 +129,11 @@ class Population:
         # Custom options
         self.custom_options = {}
 
+        # grid code generation
+        self.indent_depth = 0
+        self.indent_string = "    "
+        self.code_string = ""
+
         # Set the options that are passed at creation of the object
         self.set(**kwargs)
 
@@ -1930,6 +1935,38 @@ class Population:
     # a variable grid
     ###################################################
 
+    def _add_code(self,*args,indent=0,mindent=-1):
+        # add code to the code_string
+        #
+        # indent (=0) is added once at the beginning
+        # mindent (=0) is added for every line
+        #
+        # don't use both!
+
+        # if we have multiple arguments and have forgotten
+        # to set mindent = 0, do it here
+        if indent == 0 and mindent == -1 and len(args) > 0:
+            mindent = 0
+
+        # if mindent isn't set, use indent
+        if mindent == -1:
+            self.code_string += self._indent_block(indent)
+
+        # put everything in the code_string, indenting
+        # each line if required
+        for thing in args:
+            if mindent != -1 or len(args) > 1:
+                self.code_string += self._indent_block(mindent)
+            self.code_string += thing
+
+    def _indent_block(self,n=0):
+        # return an indent block, with n extra blocks in it
+        return (self.indent_depth + n) * self.indent_string
+
+    def _increment_indent_depth(self,delta):
+        # increment the indent indent_depth by delta
+        self.indent_depth += delta
+
     def _generate_grid_code(self, dry_run=False):
         """
         Function that generates the code from which the population will be made.
@@ -1955,91 +1992,86 @@ class Population:
 
         verbose_print("Generating grid code", self.grid_options["verbosity"], 1)
 
-        # Some local values
-        code_string = ""
-        depth = 0
-        indent = "    "
         total_grid_variables = len(self.grid_options["_grid_variables"])
 
-        # Import packages
-        code_string += "import math\n"
-        code_string += "import numpy as np\n"
-        code_string += "from collections import OrderedDict\n"
-        code_string += "from binarycpython.utils.distribution_functions import *\n"
-        code_string += "from binarycpython.utils.spacing_functions import *\n"
-        code_string += "from binarycpython.utils.useful_funcs import *\n"
-        code_string += "\n\n"
-
-        # Make the function
-        code_string += "def grid_code(self, print_results=True):\n"
-
-        # Increase depth
-        depth += 1
-
-        # Write some info in the function
-        code_string += (
-            indent * depth
-            + "# Grid code generated on {}\n".format(
-                datetime.datetime.now().isoformat()
-            )
-            + indent * depth
-            + "# This function generates the systems that will be evolved with binary_c\n\n"
-        )
 
-        # Set some values in the generated code:
-        code_string += indent * depth + "# Setting initial values\n"
-        code_string += indent * depth + "_total_starcount = 0\n"
-        code_string += indent * depth + "starcounts = [0 for i in range({})]\n".format(
-            total_grid_variables
-        )
-        code_string += indent * depth + "probabilities = {}\n"
-        code_string += (
-            indent * depth
-            + "probabilities_list = [0 for i in range({})]\n".format(
+        self._add_code(
+            # Import packages
+            "import math\n",
+            "import numpy as np\n",
+            "from collections import OrderedDict\n",
+            "from binarycpython.utils.distribution_functions import *\n",
+            "from binarycpython.utils.spacing_functions import *\n",
+            "from binarycpython.utils.useful_funcs import *\n",
+            "\n\n",
+            # Make the function
+            "def grid_code(self, print_results=True):\n")
+
+        # Increase indent_depth
+        self._increment_indent_depth(+1)
+
+        self._add_code(
+            # Write some info in the function
+            "# Grid code generated on {}\n".format(
+                datetime.datetime.now().isoformat()
+            ),
+            "# This function generates the systems that will be evolved with binary_c\n\n"
+            # Set some values in the generated code:
+            "# Setting initial values\n",
+            "_total_starcount = 0\n",
+            "starcounts = [0 for i in range({})]\n".format(
                 total_grid_variables
-            )
-        )
-        code_string += (
-            indent * depth
-            + "probabilities_sum = [0 for i in range({})]\n".format(
+            ),
+            "probabilities = {}\n",
+            "probabilities_list = [0 for i in range({})]\n".format(
                 total_grid_variables
-            )
+            ),
+            "probabilities_sum = [0 for i in range({})]\n".format(
+                total_grid_variables
+            ),
+            "parameter_dict = {}\n",
+            "phasevol = 1\n",
+            "\n"
+        )
+
+        # Set up the system parameters
+        self._add_code(
+            "M_1 = None\n",
+            "M_2 = None\n",
+            "M_3 = None\n",
+            "M_4 = None\n",
+
+            "orbital_period = None\n",
+            "orbital_period_triple = None\n",
+            "orbital_period_quadruple = None\n",
+
+            "eccentricity = None\n",
+            "eccentricity2 = None\n",
+            "eccentricity3 = None\n",
+            "\n",
+
+            # Prepare the probability
+            "# setting probability lists\n"
         )
-        code_string += indent * depth + "parameter_dict = {}\n"
-        code_string += indent * depth + "phasevol = 1\n"
-        code_string += indent * depth + "\n"
-
-        # Setting up the system parameters
-        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 + "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"
-        code_string += indent * depth + "eccentricity3 = None\n"
-        code_string += indent * depth + "\n"
-
-        # Prepare the probability
-        code_string += indent * depth + "# setting probability lists\n"
+
         for grid_variable_el in sorted(
             self.grid_options["_grid_variables"].items(),
             key=lambda x: x[1]["grid_variable_number"],
         ):
             # Make probabilities dict
             grid_variable = grid_variable_el[1]
-            code_string += indent * depth + 'probabilities["{}"] = 0\n'.format(
+            self._add_code('probabilities["{}"] = 0\n'.format(
                 grid_variable["parameter_name"]
-            )
+            ))
 
         #################################################################################
         # Start of code generation
         #################################################################################
-        code_string += indent * depth + "\n"
+        self._add_code("\n")
+
+        # turn vb to True to have debugging output
+        vb = False
+
         # Generate code
         print("Generating grid code")
         for loopnr, grid_variable_el in enumerate(
@@ -2054,92 +2086,165 @@ class Population:
             #########################
             # Setting up the for loop
             # Add comment for for loop
-            code_string += (
-                indent * depth
-                + "# for loop for {}".format(grid_variable["parameter_name"])
-                + "\n"
-            )
-
-            code_string += (
-                indent * depth
-                + "sampled_values_{} = {}".format(
+            self._add_code(
+                "# for loop for {}".format(grid_variable["parameter_name"])
+                + "\n",
+                "sampled_values_{} = {}".format(
                     grid_variable["name"], grid_variable["spacingfunc"]
                 )
                 + "\n"
             )
 
-            # TODO: Make clear that the phasevol only works good
-            # TODO: add option to ignore this phasevol calculation and set it to 1
-            #   if you sample linearly in that thing.
-            # if phasevol is <= 0 then we SKIP that whole loop. Its not working then.
-            if (
-                not grid_variable["dphasevol"] == -1
-            ):  # A method to turn off this calculation and allow a phasevol = 1
-                code_string += (
-                    indent * depth
-                    + "phasevol_{} = sampled_values_{}[1]-sampled_values_{}[0]".format(
-                        grid_variable["name"],
-                        grid_variable["name"],
-                        grid_variable["name"],
-                    )
-                    + "\n"
-                )
-            else:
-                code_string += indent * depth + "phasevol_{} = 1\n".format(
-                    grid_variable["name"]
+            if vb:
+                self._add_code(
+                    "print('samples','{}',':',sampled_values_{})\n".format(grid_variable["name"],grid_variable["name"])
                 )
 
             # # Some print statement
-            # code_string += (
-            #     indent * depth
-            #     + "print('phasevol_{}:', phasevol_{})".format(grid_variable["name"],
+            # self._add_code((
+            #     "print('phasevol_{}:', phasevol_{})".format(grid_variable["name"],
             #           grid_variable["name"])
             #     + "\n"
             # )
 
             # TODO: make sure this works
             # Adding for loop structure
-            # code_string += (
-            #     indent * depth
-            #     + "for {} in sampled_values_{}:".format(
+            # self._add_code((
+            #     "for {} in sampled_values_{}:".format(
             #         grid_variable["name"], grid_variable["name"]
             #     )
             #     + "\n"
             # )
 
+            if vb:
+                self._add_code(
+                    "print('sample {} from',sampled_values_{})".format(grid_variable["name"],grid_variable["name"])
+                    + "\n"
+                )
 
             if grid_variable["resolution"] == 0:
                 # generate resolution automatically
-                code_string += (
-                    indent * depth
-                    + "for {}_sample_number in range(len(sampled_values_{})):".format(
+                self._add_code(
+                    "for {}_sample_number in range(max(0,len(sampled_values_{})-1)):".format(
                         grid_variable["name"], grid_variable["name"]
                     )
                     + "\n"
                 )
             else:
                 # use resolution passed in
-                code_string += (
-                    indent * depth
-                    + "for {}_sample_number in range({}):".format(
-                        grid_variable["name"], grid_variable["resolution"]
+                self._add_code(
+                    "for {}_sample_number in range({}):".format(
+                        grid_variable["name"],
+                        grid_variable["resolution"]
+                    )
+                    + "\n"
+                )
+
+            self._increment_indent_depth(+1)
+
+            # {}_this and {}_next are this grid point's index and the next
+            # grid point's index, used in many calculations below
+            self._add_code("if {}_sample_number == 0:\n".format(grid_variable["name"]))
+            self._add_code(
+                "{}_this = 0;\n".format(grid_variable["name"]),
+                indent=1
+            )
+            self._add_code("else:\n")
+            self._add_code(
+                "{}_this = {}_sample_number; ".format(grid_variable["name"], grid_variable["name"]),
+                indent=1
+            )
+            self._add_code(
+                "\n",
+                "{}_next = {}_this + 1".format(grid_variable["name"], grid_variable["name"]),
+                "\n")
+
+            # TODO: Make clear that the phasevol only works good
+            # TODO: add option to ignore this phasevol calculation and set it to 1
+            #   if you sample linearly in that thing.
+            # if phasevol is <= 0 then we SKIP that whole loop. Its not working then.
+            if (
+                    not grid_variable["dphasevol"] == -1
+            ):  # A method to turn off this calculation and allow a phasevol = 1
+                self._add_code(
+                    "phasevol_{} = sampled_values_{}[{}_next]-sampled_values_{}[{}_this]".format(
+                        grid_variable["name"],
+                        grid_variable["name"],
+                        grid_variable["name"],
+                        grid_variable["name"],
+                        grid_variable["name"]
+                    )
+                    + "\n")
+            else:
+                self._add_code(
+                    "phasevol_{} = 1\n".format(
+                        grid_variable["name"]
                     )
+                )
+
+            ##############
+            # Add phasevol check:
+            self._add_code(
+                "if phasevol_{} <= 0:\n".format(grid_variable["name"])
+            )
+
+            # TODO: We might actually want to add the starcount and probability to the totals regardless.
+            #   n that case we need another local variable which will prevent it from being run but will track those parameters
+            # Add phasevol check action:
+            self._add_code(
+                'print("Grid generator: phasevol_{} <= 0! (this=",{}_this,", next=",{}_next,") Skipping current sample.")'.format(
+                    grid_variable["name"],
+                    grid_variable["name"],
+                    grid_variable["name"]
+                )
+                + "\n",
+                "continue" + "\n",
+                mindent=1
+            )
+
+            if vb:
+                self._add_code(
+                    "print('sample {} from ',sampled_values_{},' at this=',{}_this,', next=',{}_next)".format(
+                        grid_variable["name"],
+                        grid_variable["name"],
+                        grid_variable["name"],
+                        grid_variable["name"]
+                        )
                     + "\n"
                 )
 
-            if grid_variable['gridtype'] == 'edge':
-                code_string += (
-                    indent * (depth+1)
-                    + "{} = sampled_values_{}[0] + (sampled_values_{}[1]-sampled_values_{}[0]) * {}_sample_number".format(
-                        grid_variable["name"], grid_variable["name"], grid_variable["name"], grid_variable["name"], grid_variable["name"]
+            # select point location based on gridtype (left, centre or right)
+            if (grid_variable['gridtype'] == 'edge' or
+                grid_variable['gridtype'] == 'left' or
+                grid_variable['gridtype'] == 'left edge'):
+                self._add_code(
+                    "{} = sampled_values_{}[{}_this]".format(
+                        grid_variable["name"],
+                        grid_variable["name"],
+                        grid_variable["name"]
                     )
                     + "\n"
                 )
-            elif grid_variable['gridtype'] == 'centred':
-                code_string += (
-                    indent * (depth+1)
-                    + "{} = sampled_values_{}[0] + 0.5 * (sampled_values_{}[1]-sampled_values_{}[0]) + (sampled_values_{}[1]-sampled_values_{}[0]) * {}_sample_number".format(
-                        grid_variable["name"], grid_variable["name"], grid_variable["name"], grid_variable["name"], grid_variable["name"], grid_variable["name"], grid_variable["name"]
+            elif (grid_variable['gridtype'] == 'right' or
+                  grid_variable['gridtype'] == 'right edge'):
+                self._add_code(
+                    + "{} = sampled_values_{}[{}_next]".format(
+                        grid_variable["name"],
+                        grid_variable["name"],
+                        grid_variable["name"]
+                    )
+                    + "\n"
+                )
+            elif (grid_variable['gridtype'] == 'centred' or
+                  grid_variable['gridtype'] == 'centre' or
+                  grid_variable['gridtype'] == 'center'):
+                self._add_code(
+                    "{} = 0.5 * (sampled_values_{}[{}_next] + sampled_values_{}[{}_this])".format(
+                        grid_variable["name"],
+                        grid_variable["name"],
+                        grid_variable["name"],
+                        grid_variable["name"],
+                        grid_variable["name"]
                     )
                     + "\n"
                 )
@@ -2147,124 +2252,95 @@ class Population:
                 msg = "Unknown gridtype value. PLease choose a different one"
                 raise ValueError(msg)
 
+            if vb:
+                self._add_code(
+                    "print('hence {} = ',{})\n".format(grid_variable["name"], grid_variable["name"])
+                )
+
             #################################################################################
             # Check condition and generate for loop
 
             # If the grid variable has a condition, write the check and the action
             if grid_variable["condition"]:
-                # Add comment
-                code_string += (
-                    indent * (depth + 1)
-                    + "# Condition for {}".format(grid_variable["parameter_name"])
-                    + "\n"
-                )
 
-                # Add condition check
-                code_string += (
-                    indent * (depth + 1)
-                    + "if not {}:".format(grid_variable["condition"])
-                    + "\n"
+                self._add_code(
+                    # Add comment
+                    "# Condition for {}\n".format(grid_variable["parameter_name"]),
+                    # Add condition check
+                    "if not {}:\n".format(grid_variable["condition"]),
+                    mindent=0
                 )
 
                 # Add condition failed action:
-                code_string += (
-                    indent * (depth + 2)
-                    + 'print("Grid generator: Condition for {} not met!")'.format(
+                self._add_code(
+                    'print("Grid generator: Condition for {} not met!")'.format(
                         grid_variable["parameter_name"]
                     )
-                    + "\n"
-                )
-                code_string += indent * (depth + 2) + "continue" + "\n"
+                    + "\n",
+                    "continue" + "\n",
+                    mindent=1)
 
                 # Add some white line
-                code_string += indent * (depth + 1) + "\n"
+                self._add_code("\n")
 
             # # Add debugging:
             # if grid_variable['name']=='q':
-            #     code_string += (
-            #         indent * (depth + 1)
-            #         + 'print("sampling:", sampled_values_{}, M_1)'.format(
+            #     self._add_code((
+            #         indent=1,
+            #         'print("sampling:", sampled_values_{}, M_1)'.format(
             #             grid_variable["name"], grid_variable["name"]
             #         )
             #         + "\n"
             #     )
 
-            ##############3
-            # Add phasevol check:
-            code_string += (
-                indent * (depth + 1)
-                + "if phasevol_{} <= 0:".format(grid_variable["name"])
-                + "\n"
-            )
-
-            # TODO: We might actually want to add the starcount and probability to the totals regardless.
-            #   n that case we need another local variable which will prevent it from being run but will track those parameters
-            # Add phasevol check action:
-            code_string += (
-                indent * (depth + 2)
-                + 'print("Grid generator: phasevol_{} <= 0! Skipping current sample.")'.format(
-                    grid_variable["name"]
-                )
-                + "\n"
-            )
-            code_string += indent * (depth + 2) + "continue" + "\n"
-
-            # Add some white line
-            code_string += indent * (depth + 1) + "\n"
+            # Add some whitespace
+            self._add_code("\n")
 
             #########################
             # Setting up pre-code and value in some cases
             # Add pre-code
             if grid_variable["precode"]:
-                code_string += (
-                    indent * (depth + 1)
-                    + "{}".format(
+                self._add_code(
+                    "{}".format(
                         grid_variable["precode"].replace(
-                            "\n", "\n" + indent * (depth + 1)
+                            "\n", "\n" + self._indent_block(0)
                         )
                     )
                     + "\n"
                 )
 
             # Set phasevol
-            code_string += indent * (depth + 1) + "phasevol *= phasevol_{}\n".format(
+            self._add_code("phasevol *= phasevol_{}\n".format(
                 grid_variable["name"],
-            )
+            ))
 
             #######################
             # Probabilities
             # Calculate probability
-            code_string += indent * (depth + 1) + "\n"
-            code_string += indent * (depth + 1) + "# Setting probabilities\n"
-            code_string += (
-                indent * (depth + 1)
-                + "d{} = phasevol_{} * ({})".format(
+            self._add_code(
+                "\n",
+                "# Setting probabilities\n",
+                "d{} = phasevol_{} * ({})".format(
                     grid_variable["name"],
                     grid_variable["name"],
                     grid_variable["probdist"],
                 )
-                + "\n"
-            )
-
-            # Saving probability sum
-            code_string += (
-                indent * (depth + 1)
-                + "probabilities_sum[{}] += d{}".format(
-                    grid_variable["grid_variable_number"], grid_variable["name"]
+                + "\n",
+                # Save probability sum
+                "probabilities_sum[{}] += d{}".format(
+                grid_variable["grid_variable_number"], grid_variable["name"]
                 )
                 + "\n"
             )
 
             if grid_variable["grid_variable_number"] == 0:
-                code_string += (
-                    indent * (depth + 1)
-                    + "probabilities_list[0] = d{}".format(grid_variable["name"])
+                self._add_code(
+                    "probabilities_list[0] = d{}".format(grid_variable["name"])
                     + "\n"
                 )
             else:
-                code_string += (
-                    indent * (depth + 1)
-                    + "probabilities_list[{}] = probabilities_list[{}] * d{}".format(
+                self._add_code(
+                    "probabilities_list[{}] = probabilities_list[{}] * d{}".format(
                         grid_variable["grid_variable_number"],
                         grid_variable["grid_variable_number"] - 1,
                         grid_variable["name"],
@@ -2274,51 +2350,40 @@ class Population:
 
             #######################
             # Increment starcount for this parameter
-            code_string += "\n"
-            code_string += indent * (
-                depth + 1
-            ) + "# Increment starcount for {}\n".format(grid_variable["parameter_name"])
-
-            code_string += (
-                indent * (depth + 1)
-                + "starcounts[{}] += 1".format(
-                    grid_variable["grid_variable_number"],
-                )
-                + "\n"
-            )
-
-            # Add value to dict
-            code_string += (
-                indent * (depth + 1)
-                + 'parameter_dict["{}"] = {}'.format(
-                    grid_variable["parameter_name"], grid_variable["parameter_name"]
-                )
-                + "\n"
+            self._add_code("\n",
+                           "# Increment starcount for {}\n".format(grid_variable["parameter_name"]),
+
+                           "starcounts[{}] += 1".format(
+                               grid_variable["grid_variable_number"],
+                           )
+                           + "\n",
+
+                           # Add value to dict
+                           'parameter_dict["{}"] = {}'.format(
+                               grid_variable["parameter_name"], grid_variable["parameter_name"]
+                           )
+                           + "\n",
+                           "\n"
             )
 
-            # Add some space
-            code_string += "\n"
+            self._increment_indent_depth(-1)
 
             # The final parts of the code, where things are returned, are within the deepest loop,
             # but in some cases code from a higher loop needs to go under it again
             # SO I think its better to put an if statement here that checks
             # whether this is the last loop.
             if loopnr == len(self.grid_options["_grid_variables"]) - 1:
-
-                code_string = self._write_gridcode_system_call(
-                    code_string,
-                    indent,
-                    depth,
+                self._write_gridcode_system_call(
                     grid_variable,
                     dry_run,
                     grid_variable["branchpoint"],
                 )
 
-            # increment depth
-            depth += 1
+            # increment indent_depth
+            self._increment_indent_depth(+1)
 
-        depth -= 1
-        code_string += "\n"
+        self._increment_indent_depth(-1)
+        self._add_code("\n")
 
         # Write parts to write below the part that yield the results.
         # this has to go in a reverse order:
@@ -2332,80 +2397,74 @@ class Population:
         for loopnr, grid_variable_el in enumerate(reverse_sorted_grid_variables):
             grid_variable = grid_variable_el[1]
 
-            code_string += indent * (depth + 1) + "#" * 40 + "\n"
-            code_string += (
-                indent * (depth + 1)
-                + "# Code below is for finalising the handling of this iteration of the parameter {}\n".format(grid_variable["name"])
+            self._increment_indent_depth(+1)
+            self._add_code(
+                "#" * 40 + "\n",
+                "# Code below is for finalising the handling of this iteration of the parameter {}\n".format(grid_variable["name"])
             )
 
             # Set phasevol
             # TODO: fix. this isn't supposed to be the value that we give it here. discuss
-            code_string += indent * (depth + 1) + "phasevol /= phasevol_{}\n".format(
-                grid_variable["name"]
+            self._add_code(
+                "phasevol /= phasevol_{}\n\n".format(grid_variable["name"])
             )
-            code_string += indent * (depth + 1) + "\n"
 
-            depth -= 1
+            self._increment_indent_depth(-2)
 
             # Check the branchpoint part here. The branchpoint makes sure that we can construct
             # a grid with several multiplicities and still can make the system calls for each
             # multiplicity without reconstructing the grid each time
             if grid_variable["branchpoint"] > 0:
 
-                # Add comment
-                code_string += (
-                    indent * (depth + 1)
-                    + "# Condition for branchpoint at {}".format(
+                self._increment_indent_depth(+1)
+
+
+                self._add_code(
+                    # Add comment
+                    "# Condition for branchpoint at {}".format(
                         reverse_sorted_grid_variables[loopnr + 1][1]["parameter_name"]
                     )
-                    + "\n"
-                )
+                    + "\n",
+                    # # Add condition check
+                    #     "if not {}:".format(grid_variable["condition"])
+                    #     + "\n"
 
-                # # Add condition check
-                # code_string += (
-                #     indent * (depth + 1)
-                #     + "if not {}:".format(grid_variable["condition"])
-                #     + "\n"
-                # )
-
-                # Add branchpoint
-                code_string += (
-                    indent * (depth + 1)
-                    + "if multiplicity=={}:".format(grid_variable["branchpoint"])
+                    # Add branchpoint
+                    "if multiplicity=={}:".format(grid_variable["branchpoint"])
                     + "\n"
                 )
 
-                code_string = self._write_gridcode_system_call(
-                    code_string,
-                    indent,
-                    depth + 1,
+                self._write_gridcode_system_call(
                     reverse_sorted_grid_variables[loopnr + 1][1],
                     dry_run,
                     grid_variable["branchpoint"],
                 )
-                code_string += "\n"
+                self._increment_indent_depth(-1)
+                self._add_code("\n")
 
-        ################
+        ###############################
         # Finalising print statements
         #
-        # code_string += indent * (depth + 1) + "\n"
-        code_string += indent * (depth + 1) + "#" * 40 + "\n"
-        code_string += indent * (depth + 1) + "if print_results:\n"
-        code_string += (
-            indent * (depth + 2)
-            + "print('Grid has handled {} stars'.format(_total_starcount))\n"
+        self._increment_indent_depth(+1)
+        self._add_code(
+            "\n",
+            "#" * 40 + "\n",
+            "if print_results:\n"
         )
-        code_string += (
-            indent * (depth + 2)
-            + "print('with a total probability of {}'.format(self.grid_options['_probtot']))\n"
+        self._add_code(
+            "print('Grid has handled {} stars'.format(_total_starcount))\n",
+            "print('with a total probability of {}'.format(self.grid_options['_probtot']))\n",
+            mindent=1
         )
 
+
         ################
         # Finalising return statement for dry run.
         #
         if dry_run:
-            code_string += indent * (depth + 1) + "return _total_starcount\n"
+            self._add_code("return _total_starcount\n")
 
+        self._increment_indent_depth(-1)
         #################################################################################
         # Stop of code generation. Here the code is saved and written
 
@@ -2414,7 +2473,7 @@ class Population:
             "Saving grid code to grid_options", self.grid_options["verbosity"], 1
         )
 
-        self.grid_options["code_string"] = code_string
+        self.grid_options["code_string"] = self.code_string
 
         # Write to file
         gridcode_filename = os.path.join(
@@ -2430,125 +2489,83 @@ class Population:
         )
 
         with open(gridcode_filename, "w") as file:
-            file.write(code_string)
+            file.write(self.code_string)
 
 
     def _write_gridcode_system_call(
-        self, code_string, indent, depth, grid_variable, dry_run, branchpoint
+            self, grid_variable, dry_run, branchpoint
     ):
         #################################################################################
         # Here are the calls to the queuing or other solution. this part is for every system
         # Add comment
-        code_string += indent * (depth + 1) + "#" * 40 + "\n"
+        self._increment_indent_depth(+1)
+        self._add_code("#" * 40 + "\n")
 
         if branchpoint:
-            code_string += (
-                indent * (depth + 1)
-                + "# Code below will get evaluated for every system at this level of multiplicity (last one of that being {})\n"
-            ).format(grid_variable["name"])
+            self._add_code(
+                "# Code below will get evaluated for every system at this level of multiplicity (last one of that being {})\n".format(grid_variable["name"])
+            )
         else:
-            code_string += (
-                indent * (depth + 1)
-                + "# Code below will get evaluated for every generated system\n"
+            self._add_code(
+                "# Code below will get evaluated for every generated system\n"
             )
 
         # Factor in the custom weight input
-        code_string += (
-            indent * (depth + 1)
-            + "\n"
-        )
-        code_string += (
-            indent * (depth + 1)
-            + "# Weigh the probability by a custom weighting factor\n"
-        )
-        code_string += (
-            indent * (depth + 1)
-            + 'probability = self.grid_options["weight"] * probabilities_list[{}]'.format(
+        self._add_code(
+            "\n",
+            "# Weigh the probability by a custom weighting factor\n",
+            'probability = self.grid_options["weight"] * probabilities_list[{}]'.format(
                 grid_variable["grid_variable_number"]
             )
-            + "\n"
-        )
-
-        # Take into account the multiplicity fraction:
-        code_string += (
-            indent * (depth + 1)
-            + "\n"
-        )
-        code_string += (
-            indent * (depth + 1)
-            + "# Factor the multiplicity fraction into the probability\n"
-        )
-        code_string += (
-            indent * (depth + 1)
-            + 'probability = probability * self._calculate_multiplicity_fraction(parameter_dict)'
-            + "\n"
-        )
-
-        # Add division by amount of repeats
-        code_string += (
-            indent * (depth + 1)
-            + "\n"
-        )
-        code_string += (
-            indent * (depth + 1)
-            + "# Divide the probability by the amount of repeats\n"
-        )
-
-        code_string += (
-            indent * (depth + 1)
-            + 'probability = probability / self.grid_options["repeat"]'
-            + "\n"
-        )
-
-        # Now we yield the system self.grid_options["repeat"] times.
-        code_string += (
-            indent * (depth + 1)
-            + "\n"
-        )
-        code_string += (
-            indent * (depth + 1)
-            + "# Loop over the repeats\n"
-        )
-        code_string += (
-            indent * (depth + 1) + 'for _ in range(self.grid_options["repeat"]):' + "\n"
-        )
+            + "\n",
 
-        code_string += indent * (depth + 2) + "_total_starcount += 1\n"
 
-        # set probability and phasevol values into the system dict
-        code_string += (
-            indent * (depth + 2)
-            + 'parameter_dict["{}"] = {}'.format("probability", "probability")
-            + "\n"
+        # Take into account the multiplicity fraction:
+            "\n",
+            "# Factor the multiplicity fraction into the probability\n",
+            'probability = probability * self._calculate_multiplicity_fraction(parameter_dict)'
+            + "\n",
+
+
+            # Add division by amount of repeats
+            "\n",
+            "# Divide the probability by the amount of repeats\n",
+            'probability = probability / self.grid_options["repeat"]'
+            + "\n",
+            # Now we yield the system self.grid_options["repeat"] times.
+            "\n",
+            "# Loop over the repeats\n",
+            'for _ in range(self.grid_options["repeat"]):' + "\n"
         )
-        code_string += (
-            indent * (depth + 2)
-            + 'parameter_dict["{}"] = {}'.format("phasevol", "phasevol")
-            + "\n"
+        self._add_code(
+            "_total_starcount += 1\n",
+            # set probability and phasevol values into the system dict
+            'parameter_dict["{}"] = {}'.format("probability", "probability") + "\n",
+            'parameter_dict["{}"] = {}'.format("phasevol", "phasevol") + "\n",
+            # Increment total probability
+            "self._increment_probtot(probability)\n",
+            mindent=1
         )
 
-        # Some prints. will be removed
-        # code_string += indent * (depth + 1) + "print(probabilities)\n"
-        # code_string += (
-        #     indent * (depth + 1) + 'print("_total_starcount: ", _total_starcount)\n'
-        # )
-
-        # code_string += indent * (depth + 1) + "print(probability)\n"
-
-        # Increment total probability
-        code_string += indent * (depth + 2) + "self._increment_probtot(probability)\n"
-
         if not dry_run:
             # Handling of what is returned, or what is not.
-            code_string += indent * (depth + 2) + "yield(parameter_dict)\n"
+            self._add_code(
+                "yield(parameter_dict)\n",
+                indent=1
+            )
 
         # If its a dry run, dont do anything with it
         else:
-            code_string += indent * (depth + 2) + "pass\n"
+            self._add_code(
+                "pass\n",
+                indent=1
+            )
+
+        self._add_code("#" * 40 + "\n")
 
-        code_string += indent * (depth + 1) + "#" * 40 + "\n"
+        self._increment_indent_depth(-1)
 
-        return code_string
+        return self.code_string
 
     def _load_grid_function(self):
         """
diff --git a/binarycpython/utils/spacing_functions.py b/binarycpython/utils/spacing_functions.py
index 0f96a1c9e3ff61671fad68f206c63ddbe41ac388..e83f91048ee6e35ad54ff5927baed25d19cd9e24 100644
--- a/binarycpython/utils/spacing_functions.py
+++ b/binarycpython/utils/spacing_functions.py
@@ -6,8 +6,9 @@ Tasks:
 """
 
 from typing import Union
-import numpy as np
 import math
+import numpy as np
+import sys
 
 def const(
     min_bound: Union[int, float], max_bound: Union[int, float], steps: int