diff --git a/VERSION b/VERSION
index 2003b639c40025a4216b7b765e800b872a9052cd..965065db5b84c28a60e0e932153c4f59ac44c4da 100644
--- a/VERSION
+++ b/VERSION
@@ -1 +1 @@
-0.9.2
+0.9.3
diff --git a/binarycpython/utils/distribution_functions.py b/binarycpython/utils/distribution_functions.py
index d6a99452af87c021707bd06012bae4695ca087f9..f5d6db84f0dae9943f62fbed43d3681cd788cafe 100644
--- a/binarycpython/utils/distribution_functions.py
+++ b/binarycpython/utils/distribution_functions.py
@@ -1894,7 +1894,7 @@ def Moe_di_Stefano_2017_pdf(options, verbosity=0):
         multiplicity = options["multiplicity"]
 
     # Immediately return 0 if the multiplicity modulator is 0
-    if options["multiplicity_modulator"][multiplicity - 1] == 0:
+    if options["multiplicity_modulator"][int(multiplicity) - 1] == 0:
         verbose_print(
             "\tMoe_di_Stefano_2017_pdf: returning 0 because of the multiplicity modulator being 0",
             verbosity,
@@ -1908,7 +1908,7 @@ def Moe_di_Stefano_2017_pdf(options, verbosity=0):
 
     multiplicity_probability = Moe_di_Stefano_2017_multiplicity_fractions(
         options, verbosity
-    )[multiplicity - 1]
+    )[int(multiplicity) - 1]
     prob_dict["multiplicity"] = multiplicity_probability
     verbose_print(
         "\tMoe_di_Stefano_2017_pdf: Appended multiplicity (mass1 = {}) probability ({}) to the prob dict ({})".format(
diff --git a/binarycpython/utils/grid.py b/binarycpython/utils/grid.py
index 39e0919c0cf919f574c418d8d0ce4bf78a43a0dc..690ce34a24489727fd7ef861c94ae690f481ffdf 100644
--- a/binarycpython/utils/grid.py
+++ b/binarycpython/utils/grid.py
@@ -2363,8 +2363,7 @@ class Population:
                 "sampled_values_{} = {}".format(
                     grid_variable["name"], grid_variable["samplerfunc"]
                 )
-                + "\n",
-            )
+                 + "\n")
 
             if vb:
                 self._add_code(
@@ -2397,11 +2396,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"
             )
@@ -2448,13 +2482,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:
@@ -2468,15 +2496,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,
@@ -2484,12 +2504,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"
                 )
 
@@ -2500,22 +2515,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"
                 )
@@ -2525,13 +2540,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:
@@ -2550,10 +2559,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,
@@ -4264,11 +4273,11 @@ 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"]
             ),
-            condition="({}[multiplicity-1] > 0)".format(
+            condition="({}[int(multiplicity)-1] > 0)".format(
                 str(self.grid_options["Moe2017_options"]["multiplicity_modulator"])
             ),
             gridtype="edge",
@@ -4921,7 +4930,7 @@ eccentricity3=0
 
         # use the multiplicity if given
         if "multiplicity" in system_dict:
-            nmult = system_dict["multiplicity"]
+            nmult = int(system_dict["multiplicity"])
         else:
             nmult = 4
 
diff --git a/binarycpython/utils/spacing_functions.py b/binarycpython/utils/spacing_functions.py
index ac4b7554e9354d491dec9fa92c90c49a76d8db72..c3e0be388abf1320089604e61e8831a722480141 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, must be an integer (is converted to int)
+        max_bound: upper bound of range, must be an integer (is converted to int)
+        steps: number of segments between min_bound and max_bound
+
+    Returns:
+        range(min_bound,max_bound,step)
+
+        where step is int((int(max_bound)-int(min_bound))/steps)
+    """
+
+    step = int((int(max_bound)-int(min_bound))/(steps-1))
+    if steps <= 1:
+        return int(min_bound)
+    else:
+        return range(int(min_bound),int(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
diff --git a/development_requirements.txt b/development_requirements.txt
index 231c873dc319e77b31115e774d1cad3be24869a0..4f93d477aed1b6d0d997d745f0c2334e77ea70bf 100644
--- a/development_requirements.txt
+++ b/development_requirements.txt
@@ -3,3 +3,4 @@ gitpython
 sphinx
 sphinx_autodoc_typehints
 nbsphinx
+twine
diff --git a/upload_to_pypi.sh b/upload_to_pypi.sh
new file mode 100755
index 0000000000000000000000000000000000000000..db8bb43d3ed5f702c3e78529df10800a5b891936
--- /dev/null
+++ b/upload_to_pypi.sh
@@ -0,0 +1,22 @@
+#!/bin/bash
+
+# Get current version
+VERSION_NUMBER=$(cat "VERSION")
+
+# Create dist
+echo "Creating source distribution for binarycpython-$VERSION_NUMBER"
+python setup.py sdist
+
+# Checking validity
+echo ""
+echo "Checking validity of for binarycpython-$VERSION_NUMBER source distribution:"
+twine check dist/binarycpython-$VERSION_NUMBER.tar.gz
+
+echo ""
+echo "Uploading binarcpython version $VERSION_NUMBER to pypi."
+read -p "Continue? y/n " -n 1 -r
+echo    # (optional) move to a new line
+if [[ $REPLY =~ ^[Yy]$ ]]
+then
+    twine upload dist/binarycpython-$VERSION_NUMBER.tar.gz
+fi