diff --git a/binarycpython/utils/grid.py b/binarycpython/utils/grid.py
index 3e6b4ab4cc96bbb8f154aa2e0f8771a91d902182..ea06fc4f567a61e6c46c26e170f70364185151a4 100644
--- a/binarycpython/utils/grid.py
+++ b/binarycpython/utils/grid.py
@@ -2353,8 +2353,7 @@ class Population:
                 "sampled_values_{} = {}".format(
                     grid_variable["name"], grid_variable["samplerfunc"]
                 )
-                + "\n",
-            )
+                 + "\n")
 
             if vb:
                 self._add_code(
@@ -2387,11 +2386,46 @@ class Population:
                     + "\n"
                 )
 
+
+            # if we're on a centred grid, we have actually
+            # one fewer grid point
+            if (
+                    grid_variable["gridtype"] == "centred"
+                    or grid_variable["gridtype"] == "centre"
+                    or grid_variable["gridtype"] == "center"
+            ):
+                offset = -1
+            else:
+                offset = 0
+
+            self._add_code(
+                "# for loop for {}".format(grid_variable["name"]) + "\n",
+                "nvalues_{} = len(sampled_values_{})+{}".format(
+                    grid_variable["name"],
+                    grid_variable["name"],
+                    offset
+                )
+                 + "\n")
+
             # if grid_variable["resolution"] == 0:
             # sample from the sampled_values
+            if vb:
+                self._add_code(
+                    "print(\"var {} values \",sampled_values_{},\" len \",len(sampled_values_{})+{},\" gridtype {} offset {}\\n\")\n".format(
+                        grid_variable["name"],
+                        grid_variable["name"],
+                        grid_variable["name"],
+                        offset,
+                        grid_variable['gridtype'],
+                        offset
+                    )
+                )
+
             self._add_code(
-                "for {}_sample_number in range(max(0,len(sampled_values_{})-1)):".format(
-                    grid_variable["name"], grid_variable["name"]
+                "for {}_sample_number in range(max(0,len(sampled_values_{})+{})):".format(
+                    grid_variable["name"],
+                    grid_variable["name"],
+                    offset,
                 )
                 + "\n"
             )
@@ -2438,13 +2472,7 @@ class Population:
                 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"],
-                    )
+                    "phasevol_{name} = (sampled_values_{name}[{name}_next]-sampled_values_{name}[{name}_this]) if {name}_next < nvalues_{name} else (sampled_values_{name}[{name}_next-1]-sampled_values_{name}[{name}_this-1])".format(name=grid_variable["name"])
                     + "\n"
                 )
             else:
@@ -2458,15 +2486,7 @@ class Population:
             #   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,"=",sampled_values_{}[{}_this],", next=",{}_next,"=",sampled_values_{}[{}_next],") Skipping current sample.")'.format(
-                    grid_variable["name"],
-                    grid_variable["name"],
-                    grid_variable["name"],
-                    grid_variable["name"],
-                    grid_variable["name"],
-                    grid_variable["name"],
-                    grid_variable["name"],
-                )
+                'print("Grid generator: phasevol_{name} <= 0! (this=",{name}_this,"=",sampled_values_{name}[{name}_this],", next=",{name}_next,"=",sampled_values_{name}[{name}_next],") Skipping current sample.")'.format(name=grid_variable["name"])
                 + "\n",
                 "continue\n",
                 indent=1,
@@ -2474,12 +2494,7 @@ class Population:
 
             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"],
-                    )
+                    "print('sample {name} from ',sampled_values_{name},' at this=',{name}_this,', next=',{name}_next)".format(name=grid_variable["name"])
                     + "\n"
                 )
 
@@ -2490,22 +2505,22 @@ class Population:
                 or grid_variable["gridtype"] == "left edge"
             ):
                 self._add_code(
-                    "{} = sampled_values_{}[{}_this]".format(
-                        grid_variable["name"],
-                        grid_variable["name"],
-                        grid_variable["name"],
-                    )
+                    "{name} = sampled_values_{name}[{name}_this]".format(
+                        name=grid_variable["name"])
                     + "\n"
                 )
+                if vb:
+                    self._add_code(
+                        "print(\"sampled_values_{name}[{name}_this = \",{name}_this,\" {name}_next= \",{name}_next,\"] = \",{name},\"\\\n\")\n".format(
+                            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"],
+                    + "{name} = sampled_values_{name}[{name}_next]".format(
+                        name=grid_variable["name"],
                     )
                     + "\n"
                 )
@@ -2515,13 +2530,7 @@ class Population:
                 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"],
-                    )
+                    "{name} = 0.5 * (sampled_values_{name}[{name}_next] + sampled_values_{name}[{name}_this])".format(name=grid_variable["name"])
                     + "\n"
                 )
             else:
@@ -2540,10 +2549,10 @@ class Population:
 
             # If the grid variable has a condition, write the check and the action
             if grid_variable["condition"]:
-
                 self._add_code(
                     # Add comment
                     "# Condition for {}\n".format(grid_variable["name"]),
+
                     # Add condition check
                     "if not {}:\n".format(grid_variable["condition"]),
                     indent=0,
@@ -4254,7 +4263,7 @@ class Population:
             parameter_name="multiplicity",
             longname="multiplicity",
             valuerange=[1, max_multiplicity],
-            samplerfunc="const(1, {n}, {n})".format(n=max_multiplicity),
+            samplerfunc="const_int(1, {n}, {n})".format(n=max_multiplicity),
             precode='self.grid_options["multiplicity"] = multiplicity; self.bse_options["multiplicity"] = multiplicity; options={}'.format(
                 self.grid_options["Moe2017_options"]
             ),
diff --git a/binarycpython/utils/spacing_functions.py b/binarycpython/utils/spacing_functions.py
index ac4b7554e9354d491dec9fa92c90c49a76d8db72..d4ae9c4c5025928282f53e89ab1873f86bcabde4 100644
--- a/binarycpython/utils/spacing_functions.py
+++ b/binarycpython/utils/spacing_functions.py
@@ -19,7 +19,7 @@ def const(
     min_bound: Union[int, float], max_bound: Union[int, float], steps: int
 ) -> list:
     """
-    Samples a range linearly. Uses numpy linspace.
+    Samples a range linearly. Uses numpy linspace, and returns an array of floats. Do NOT use this for integers.
 
     Args:
         min_bound: lower bound of range
@@ -27,9 +27,33 @@ def const(
         steps: number of segments between min_bound and max_bound
 
     Returns:
-        np.linspace(min_bound, max_bound, steps+1)
+        np.linspace(min_bound, max_bound, steps)
     """
-    return np.linspace(min_bound, max_bound, steps + 1)
+    return np.linspace(min_bound, max_bound, steps)
+
+@functools.lru_cache(maxsize=16)
+def const_int(
+    min_bound: Union[int, float], max_bound: Union[int, float], steps: int
+) -> list:
+    """
+    Samples an integer range linearly. Returns a list of ints.
+
+    Args:
+        min_bound: lower bound of range
+        max_bound: upper bound of range
+        steps: number of segments between min_bound and max_bound
+
+    Returns:
+        range(min_bound,max_bound,step)
+
+        where step is int((max_bound-min_bound)/steps)
+    """
+
+    step = int((max_bound-min_bound)/(steps-1))
+    if steps <= 1:
+        return min_bound
+    else:
+        return range(min_bound,max_bound+step,step)
 
 
 ############################################################
@@ -245,6 +269,11 @@ def const_dt(
     lifetime_population.evolve()
     # print("Data table",lifetime_population.grid_results['interpolation table t->m'])
 
+    if not "interpolation table t->m" in lifetime_population.grid_results or len(lifetime_population.grid_results["interpolation table t->m"].keys()) == 0:
+        print("\n\n\nError: The t->m lifetime table is empty. One usual cause for this is that the tmax or max_evolution_time option (currently passed in to const_dt as {tmax}) is too short for there to be any entries in the table before the first timestep. Try increasing tmax and max_evolution_time, shorten the timestep or, if using log times, set tstart to be closer to 0.\n".format(tmax=tmax))
+        exit()
+
+
     # convert to nested lists for the interpolator
     #
     # make time -> mass table