diff --git a/hyperion/conf/conf_files.py b/hyperion/conf/conf_files.py
index 20cd1b14..bcc73bd0 100644
--- a/hyperion/conf/conf_files.py
+++ b/hyperion/conf/conf_files.py
@@ -512,7 +512,7 @@ def set_convergence(self, convergence, percentile=100., absolute=0., relative=0.
             The relative threshold below which the ratio in the percentile
             value has to be for convergence.
         '''
-        self.check_convergence = True
+        self.check_convergence = convergence
         self.convergence_percentile = percentile
         self.convergence_absolute = absolute
         self.convergence_relative = relative
diff --git a/hyperion/conf/tests/test_run_conf.py b/hyperion/conf/tests/test_run_conf.py
index 38752721..8ae3c0e8 100644
--- a/hyperion/conf/tests/test_run_conf.py
+++ b/hyperion/conf/tests/test_run_conf.py
@@ -22,3 +22,12 @@ def test_propagation_check_frequency_invalid2(value):
     with pytest.raises(ValueError) as exc:
         r.set_propagation_check_frequency(value)
     assert exc.value.args[0] == "frequency should be between 0 and 1"
+
+
+def test_set_convergence_flag():
+    r = RunConf()
+    r.set_convergence(True, percentile=99., absolute=2., relative=1.02)
+    assert r.check_convergence is True
+    r = RunConf()
+    r.set_convergence(False)
+    assert r.check_convergence is False
diff --git a/hyperion/densities/alpha_disk.py b/hyperion/densities/alpha_disk.py
index ba443d16..a1df4959 100644
--- a/hyperion/densities/alpha_disk.py
+++ b/hyperion/densities/alpha_disk.py
@@ -448,6 +448,9 @@ def density(self, grid):
         # Find density factor
         rho *= self.rho_0
 
+        if np.sum(rho * grid.volumes) == 0. and self.mass > 0:
+            raise Exception("Discretized disk mass is zero, suggesting that the grid is too coarse")
+
         norm = self.mass / np.sum(rho * grid.volumes)
 
         logger.info("Disk density is being re-scaled by a factor of %.2f to give the correct mass." % norm)
diff --git a/hyperion/densities/ambient_medium.py b/hyperion/densities/ambient_medium.py
index fcec9a18..22a54e07 100644
--- a/hyperion/densities/ambient_medium.py
+++ b/hyperion/densities/ambient_medium.py
@@ -132,8 +132,8 @@ def dust(self, value):
 
     def _check_all_set(self):
 
-        if self.density is None:
-            raise Exception("density is not set")
+        if self.rho is None:
+            raise Exception("rho is not set")
         if self.rmin is None:
             raise Exception("rmin is not set")
         if self.rmax is None:
diff --git a/hyperion/densities/tests/test_densities.py b/hyperion/densities/tests/test_densities.py
index 443a2237..4d2e11f5 100644
--- a/hyperion/densities/tests/test_densities.py
+++ b/hyperion/densities/tests/test_densities.py
@@ -835,3 +835,56 @@ def test_ambient_densities_4():
 
     expected = np.repeat(2., 9)
     assert_array_almost_equal_nulp((p1.density(g) + c.density(g) + a.density(g))[0, 0, :], expected, 10)
+
+
+def test_ulrich_midplane_cumulative_density_inner():
+    # Regression test for a factor of 2 overestimate of the midplane
+    # cumulative density inside the centrifugal radius
+    e = UlrichEnvelope()
+    e.star = Star()
+    e.star.mass = 0.5
+    e.star.radius = 0.5
+    e.rho_0 = 1.
+    e.rc = 1.
+    e.rmin = 0.01
+    e.rmax = 10.
+    r = np.logspace(np.log10(e.rmin), np.log10(0.9 * e.rc), 100000)
+    rho_mid = e.rho_0 / np.sqrt(r / e.rc) / (1. - r / e.rc) / 2.
+    expected = np.sum(0.5 * (rho_mid[1:] + rho_mid[:-1]) * np.diff(r))
+    actual = e.midplane_cumulative_density(np.array([r[-1]]))[0]
+    assert abs(actual - expected) / expected < 1.e-3
+
+
+def test_ambient_medium_missing_rho():
+    a = AmbientMedium()
+    a.rmin = 1.
+    a.rmax = 10.
+    with pytest.raises(Exception) as exc:
+        a._check_all_set()
+    assert exc.value.args[0] == "rho is not set"
+
+
+def test_alpha_disk_zero_discretized_mass():
+    # A grid too coarse to resolve the disk gives a zero discretized mass,
+    # which used to silently produce NaN densities
+    e = AlphaDisk()
+    e.star = Star()
+    e.star.mass = 0.5
+    e.star.radius = 0.5
+    e.rmin = 1.
+    e.rmax = 2.
+    e.mass = 1.
+    e.r_0 = 5.
+    e.h_0 = 1.
+    e.p = -1.
+    e.beta = 1.25
+    e.lvisc = G * (1.5 - 2. / np.sqrt(2.) + 0.5)
+
+    r = np.array([0., 0.5, 100.])
+    t = np.linspace(0., np.pi, 3)
+    p = np.linspace(0., 2. * np.pi, 3)
+    g = SphericalPolarGrid(r, t, p)
+
+    with pytest.raises(Exception) as exc:
+        e.density(g)
+    assert "zero" in exc.value.args[0]
diff --git a/hyperion/densities/ulrich_envelope.py b/hyperion/densities/ulrich_envelope.py
index 1d6318c9..3ddcc3fc 100644
--- a/hyperion/densities/ulrich_envelope.py
+++ b/hyperion/densities/ulrich_envelope.py
@@ -442,7 +442,7 @@ def midplane_cumulative_density(self, r):
             # entries are overwritten just below, so the warnings can be
             # ignored.
             with np.errstate(invalid='ignore'):
-                rho[:] = (self.rho_0 * self.rc
+                rho[:] = (0.5 * self.rho_0 * self.rc
                           * (np.log((np.sqrt(gamma_1) + 1) / (1. - np.sqrt(gamma_1)))
                              - np.log((np.sqrt(gamma_0) + 1) / (1. - np.sqrt(gamma_0)))))
 
diff --git a/hyperion/dust/dust_type.py b/hyperion/dust/dust_type.py
index c4501463..9c012f77 100644
--- a/hyperion/dust/dust_type.py
+++ b/hyperion/dust/dust_type.py
@@ -747,7 +747,7 @@ def __init__(self, model):
             if np.any(np.isnan(values)):
                 logger.warning("NaN values found inside MieX %s file - interpolating" % quantity)
                 invalid = np.isnan(values)
-                values[invalid] = interp1d_fast_loglog(wav[~invalid], values[~invalid], wav[invalid])
+                values[invalid] = interp1d_fast_loglog(wav[~invalid][::-1], values[~invalid][::-1], wav[invalid])
                 if np.any(np.isnan(values)):
                     raise Exception("Did not manage to fix NaN values in MieX %s" % quantity)
 
@@ -808,7 +808,7 @@ def __init__(self, model):
                 if np.any(np.isnan(values[:, i])):
                     logger.warning("NaN values found inside MieX %s file - interpolating" % quantity)
                     invalid = np.isnan(values[:, i])
-                    values[:, i][invalid] = interp1d_fast_loglog(wav[~invalid], values[:, i][~invalid], wav[invalid])
+                    values[:, i][invalid] = interp1d_fast_loglog(wav[~invalid][::-1], values[:, i][~invalid][::-1], wav[invalid])
                     if np.any(np.isnan(values[:, i])):
                         raise Exception("Did not manage to fix NaN values in MieX %s" % quantity)
 
diff --git a/hyperion/filter/filter.py b/hyperion/filter/filter.py
index a8c83706..91e89f7b 100644
--- a/hyperion/filter/filter.py
+++ b/hyperion/filter/filter.py
@@ -136,7 +136,7 @@ def from_hdf5_group(cls, group, name):
 
         self.alpha = group[name].attrs['alpha']
         self._beta = group[name].attrs['beta']
-        self.central_spectral_coords = group[name].attrs['nu0'] * u.Hz
+        self.central_spectral_coord = group[name].attrs['nu0'] * u.Hz
 
         return self
 
diff --git a/hyperion/filter/tests/test_filter.py b/hyperion/filter/tests/test_filter.py
index dd483bd1..3e44a12e 100644
--- a/hyperion/filter/tests/test_filter.py
+++ b/hyperion/filter/tests/test_filter.py
@@ -104,3 +104,5 @@ def test_roundtrip():
                                f2.spectral_coord.to(u.Hz).value)
     np.testing.assert_allclose(f.transmission.to(u.percent).value,
                                f2.transmission.to(u.percent).value)
+    np.testing.assert_allclose(f.central_spectral_coord.to(u.Hz).value,
+                               f2.central_spectral_coord.to(u.Hz).value)
diff --git a/hyperion/grid/amr_grid.py b/hyperion/grid/amr_grid.py
index 54053797..6f468735 100644
--- a/hyperion/grid/amr_grid.py
+++ b/hyperion/grid/amr_grid.py
@@ -1,6 +1,7 @@
 from __future__ import print_function, division
 
 import os
+import posixpath
 import struct
 import hashlib
 from copy import deepcopy
@@ -509,13 +510,13 @@ def __setitem__(self, item, value):
                     grid_ref.quantities[item] = deepcopy(grid.quantities[value.viewed_quantity])
         elif isinstance(value, h5py.ExternalLink):
             filename = value.filename
-            base_path = os.path.dirname(value.path)
-            array_name = os.path.basename(value.path)
+            base_path = posixpath.dirname(value.path)
+            array_name = posixpath.basename(value.path)
             for ilevel, level_ref in enumerate(self.levels):
                 level_path = 'level_%05i' % (ilevel + 1)
                 for igrid, grid_ref in enumerate(level_ref.grids):
-                    grid_path = 'grid_%05i' % (ilevel + 1)
-                    grid_ref.quantities[item] = h5py.ExternalLink(filename, os.path.join(base_path, level_path, grid_path, array_name))
+                    grid_path = 'grid_%05i' % (igrid + 1)
+                    grid_ref.quantities[item] = h5py.ExternalLink(filename, posixpath.join(base_path, level_path, grid_path, array_name))
         elif value == []:
             for level in self.levels:
                 for grid in level.grids:
diff --git a/hyperion/grid/tests/test_views.py b/hyperion/grid/tests/test_views.py
index 64767007..77820ec3 100644
--- a/hyperion/grid/tests/test_views.py
+++ b/hyperion/grid/tests/test_views.py
@@ -253,3 +253,22 @@ def test_dustless_array(self, grid_type):
         g2['density1'] = g1['density1']
         g2['density2'] = g1['density1'][0]
         assert g2.n_dust == 1
+
+
+def test_amr_external_link_paths():
+    # Regression test for links that used the level index instead of the
+    # grid index in the per-grid path
+    import h5py
+    amr = AMRGrid()
+    level = amr.add_level()
+    for i in range(2):
+        grid = level.add_grid()
+        grid.xmin, grid.xmax = -1., 1.
+        grid.ymin, grid.ymax = -1., 1.
+        grid.zmin, grid.zmax = -1., 1.
+        grid.nx, grid.ny, grid.nz = 4, 4, 4
+    amr['density'] = h5py.ExternalLink('other_model.rtin', 'Grid/Quantities/density')
+    link_1 = amr.levels[0].grids[0].quantities['density']
+    link_2 = amr.levels[0].grids[1].quantities['density']
+    assert link_1.path.endswith('level_00001/grid_00001/density')
+    assert link_2.path.endswith('level_00001/grid_00002/density')
diff --git a/hyperion/grid/tests/test_voronoi.py b/hyperion/grid/tests/test_voronoi.py
index fb7de2e8..5ed2ee60 100644
--- a/hyperion/grid/tests/test_voronoi.py
+++ b/hyperion/grid/tests/test_voronoi.py
@@ -134,3 +134,16 @@ def test_init_sparse(tmpdir):
     g3 = VoronoiGrid(g2)
     g3.write(f3)
 
+
+
+def test_setitem_from_view():
+    # Regression test for VoronoiGrid.__setitem__, which used to check a
+    # nonexistent 'refined' attribute when assigning a view to an empty grid
+    g1 = VoronoiGrid([1., 2., 3.], [3., 4., 5.], [2., 3., 4.])
+    g1.quantities['density'] = np.array([1., 2., 3.])
+    g2 = VoronoiGrid()
+    g2['density'] = g1['density']
+    np.testing.assert_allclose(g2['density'].array, g1['density'].array)
+    assert g2.xmin == g1.xmin and g2.xmax == g1.xmax
+    assert g2.ymin == g1.ymin and g2.ymax == g1.ymax
+    assert g2.zmin == g1.zmin and g2.zmax == g1.zmax
diff --git a/hyperion/grid/voronoi_grid.py b/hyperion/grid/voronoi_grid.py
index d15a70a3..b2cf8ab3 100644
--- a/hyperion/grid/voronoi_grid.py
+++ b/hyperion/grid/voronoi_grid.py
@@ -545,9 +545,12 @@ def __getitem__(self, item):
 
     def __setitem__(self, item, value):
         if isinstance(value, VoronoiGridView):
-            if self.refined is None:
+            if self._x is None:
                 logger.warning("No geometry in target grid - copying from original grid")
-                self.set_points(value.x, value.y, value.z)
+                self.set_points(value.x, value.y, value.z,
+                                xmin=value.xmin, xmax=value.xmax,
+                                ymin=value.ymin, ymax=value.ymax,
+                                zmin=value.zmin, zmax=value.zmax)
             self.quantities[item] = deepcopy(value.quantities[value.viewed_quantity])
         elif isinstance(value, h5py.ExternalLink):
             self.quantities[item] = value
diff --git a/hyperion/grid/voronoi_helpers.py b/hyperion/grid/voronoi_helpers.py
index 328228de..1669bee1 100644
--- a/hyperion/grid/voronoi_helpers.py
+++ b/hyperion/grid/voronoi_helpers.py
@@ -197,12 +197,12 @@ def plot(self):
             # Drop the empty vertices coordinates, signalled by NaN.
             arr = v[~np.isnan(v)]
             assert(len(arr) % 3 == 0)
-            tmp = np.split(arr, len(arr) / 3)
+            tmp = np.split(arr, len(arr) // 3)
             # Append the vertices.
             points = points + tmp
             # Append the cell description.
             cells = cells + \
-                [len(tmp)] + range(cur_cell_idx, cur_cell_idx + len(tmp))
+                [len(tmp)] + list(range(cur_cell_idx, cur_cell_idx + len(tmp)))
             cur_cell_idx += len(tmp)
             # Append the offset info.
             offset.append(cur_offset)
diff --git a/hyperion/importers/orion.py b/hyperion/importers/orion.py
index 8492217b..96b4476a 100644
--- a/hyperion/importers/orion.py
+++ b/hyperion/importers/orion.py
@@ -58,9 +58,9 @@ def read_data(self, filename, offset, quantity_indices, verbose=False):
         p8 = header.find('))', p7)
         bytes = [int(x) for x in header[p7:p8].split()]
 
-        if bytes == range(1, n_bytes + 1):
+        if bytes == list(range(1, n_bytes + 1)):
             endian = '>'
-        elif bytes == range(n_bytes, 0, -1):
+        elif bytes == list(range(n_bytes, 0, -1)):
             endian = '<'
         else:
             raise Exception("Unexpected bytes: %s" % str(bytes))
@@ -72,7 +72,7 @@ def read_data(self, filename, offset, quantity_indices, verbose=False):
         for quantity in quantity_indices:
 
             f.seek(pos + quantity_indices[quantity] * n_bytes * gridsize)
-            array = np.fromstring(f.read(n_bytes * gridsize),
+            array = np.frombuffer(f.read(n_bytes * gridsize),
                                   dtype='%sf%i' % (endian, n_bytes))
             self.quantities[quantity] = array.reshape(self.nz, self.ny, self.nx)
 
diff --git a/hyperion/model/model.py b/hyperion/model/model.py
index 7701d0d5..b5a75c7b 100644
--- a/hyperion/model/model.py
+++ b/hyperion/model/model.py
@@ -125,8 +125,8 @@ def _read_monochromatic(self, group):
         self._monochromatic = str2bool(group.attrs['monochromatic'])
         if self._monochromatic:
             self._frequencies = np.array(group['frequencies']['nu'])
-            if 'energy_threshold' in group.attrs:
-                self._monochromatic_energy_threshold = group.attrs['energy_threshold']
+            if 'monochromatic_energy_threshold' in group.attrs:
+                self._monochromatic_energy_threshold = group.attrs['monochromatic_energy_threshold']
             else:
                 self._monochromatic_energy_threshold = 1.e-10
 
@@ -472,8 +472,8 @@ def use_image_config(self, filename):
             g_image = f['/Input/Output/']
 
         # Read in binned images
-        if 'n_theta' in g_image['Binned']:
-            self.binned_output = BinnedImageConf.read(g_image['Binned'])
+        if 'group_00001' in g_image['Binned']:
+            self.binned_output = BinnedImageConf.read(g_image['Binned']['group_00001'])
 
         # Read in peeled images
         for peeled in g_image['Peeled']:
diff --git a/hyperion/model/model_output.py b/hyperion/model/model_output.py
index 54ad8008..b367f36e 100644
--- a/hyperion/model/model_output.py
+++ b/hyperion/model/model_output.py
@@ -340,6 +340,8 @@ def get_sed(self, stokes='I', group=0, technique='peeled',
         if technique == 'peeled':
             n_groups = len(self.file['Peeled'])
             if (group < 0 and -group <= n_groups) or (group >= 0 and group < n_groups):
+                if group < 0:
+                    group = n_groups + group
                 g = self.file['Peeled/group_%05i' % (group + 1)]
             else:
                 raise ValueError('File only contains %i image/SED group(s)' % n_groups)
@@ -662,6 +664,8 @@ def get_image(self, stokes='I', group=0, technique='peeled',
         if technique == 'peeled':
             n_groups = len(self.file['Peeled'])
             if (group < 0 and -group <= n_groups) or (group >= 0 and group < n_groups):
+                if group < 0:
+                    group = n_groups + group
                 g = self.file['Peeled/group_%05i' % (group + 1)]
             else:
                 raise ValueError('File only contains %i image/SED group(s)' % n_groups)
@@ -958,6 +962,8 @@ def get_available_components(self, iteration=-1):
         # If iteration is last one, find iteration number
         if iteration == -1:
             iteration = find_last_iteration(self.file)
+        else:
+            iteration = iteration + 1  # Python value is zero based
 
         # Return components
         components = list(self.file['iteration_%05i' % iteration].keys())
diff --git a/hyperion/model/tests/test_get_quantities.py b/hyperion/model/tests/test_get_quantities.py
index 52c8ed35..ef9a9bc9 100644
--- a/hyperion/model/tests/test_get_quantities.py
+++ b/hyperion/model/tests/test_get_quantities.py
@@ -18,3 +18,22 @@ def test_no_initial(tmpdir):
     mo = m.run(tmpdir.join(random_id()).strpath)
     g = mo.get_quantities()
     assert 'density' in g
+
+
+@pytest.mark.requires_hyperion_binaries
+def test_get_available_components_zero_based(tmpdir):
+    # Regression test for get_available_components using 1-based iteration
+    # numbering while get_quantities is zero-based
+    m = Model()
+    m.set_cartesian_grid([-1., 1.], [-1., 1.], [-1., 1.])
+    m.add_density_grid(np.array([[[1.]]]), get_test_dust())
+    source = m.add_point_source()
+    source.luminosity = 1.
+    source.temperature = 6000.
+    m.set_n_initial_iterations(1)
+    m.set_n_photons(initial=1, imaging=0)
+    m.write(tmpdir.join(random_id()).strpath)
+    mo = m.run(tmpdir.join(random_id()).strpath)
+    components = mo.get_available_components(iteration=0)
+    assert 'specific_energy' in components
+    assert mo.get_available_components() == components
diff --git a/hyperion/model/tests/test_model.py b/hyperion/model/tests/test_model.py
index fa39641d..967d3cb4 100644
--- a/hyperion/model/tests/test_model.py
+++ b/hyperion/model/tests/test_model.py
@@ -5,6 +5,7 @@
 import shutil
 from copy import deepcopy
 
+import h5py
 import numpy as np
 import pytest
 
@@ -396,3 +397,39 @@ def test_binned_forced_first_interaction(tmpdir):
     with pytest.raises(Exception) as exc:
         m.write(tmpdir.join(random_id()).strpath)
     assert exc.value.args[0] == "can't use binned images with forced first interaction - use set_forced_first_interaction(False) to disable"
+
+
+def test_monochromatic_energy_threshold_round_trip(tmpdir):
+    # Regression test for the threshold being read back with the wrong
+    # attribute name and silently reset to the default
+    m = Model()
+    m.set_monochromatic(True, wavelengths=[1., 2.], energy_threshold=1.e-5)
+    f = h5py.File(tmpdir.join(random_id()).strpath, 'w')
+    m._write_monochromatic(f)
+    m2 = Model()
+    m2._read_monochromatic(f)
+    f.close()
+    assert m2._monochromatic_energy_threshold == 1.e-5
+
+
+def test_use_image_config_binned(tmpdir):
+    # Regression test for binned image configuration never being restored
+    m = Model()
+    m.set_cartesian_grid([-1., 1.], [-1., 1.], [-1., 1.])
+    m.add_density_grid(np.array([[[1.]]]), get_test_dust())
+    s = m.add_point_source()
+    s.luminosity = 1.
+    s.temperature = 6000.
+    i = m.add_binned_images(sed=True, image=False)
+    i.set_wavelength_range(5, 1, 10)
+    i.set_viewing_bins(3, 2)
+    m.set_forced_first_interaction(False)
+    m.set_n_photons(initial=1, imaging=1)
+    filename = tmpdir.join(random_id()).strpath
+    m.write(filename)
+
+    m2 = Model()
+    m2.use_image_config(filename)
+    assert m2.binned_output is not None
+    assert m2.binned_output.n_theta == 3
+    assert m2.binned_output.n_phi == 2
diff --git a/hyperion/model/tests/test_sed.py b/hyperion/model/tests/test_sed.py
index 7645dd50..6e86526e 100644
--- a/hyperion/model/tests/test_sed.py
+++ b/hyperion/model/tests/test_sed.py
@@ -451,3 +451,22 @@ def test_get_sed_stokes(self, stokes):
         with pytest.raises(ValueError) as exc:
             self.m2.get_sed(stokes=stokes)
         assert exc.value.args[0] == "Only the Stokes I value was stored for this SED"
+
+
+@pytest.mark.requires_hyperion_binaries
+def test_sed_negative_group(tmpdir):
+    # Regression test for negative group indices resolving to a nonexistent
+    # HDF5 group
+    m = Model()
+    m.set_cartesian_grid([-1., 1.], [-1., 1.], [-1., 1.])
+    s = m.add_point_source()
+    s.luminosity = 1.
+    s.temperature = 6000.
+    i = m.add_peeled_images(sed=True, image=False)
+    i.set_viewing_angles([45.], [45.])
+    i.set_wavelength_range(3, 0.1, 10.)
+    m.set_n_initial_iterations(0)
+    m.set_n_photons(imaging=100)
+    m.write(tmpdir.join(random_id()).strpath)
+    mo = m.run()
+    np.testing.assert_allclose(mo.get_sed(group=-1).val, mo.get_sed(group=0).val)
diff --git a/hyperion/sources/source.py b/hyperion/sources/source.py
index 942c3f3c..21d9d6c0 100644
--- a/hyperion/sources/source.py
+++ b/hyperion/sources/source.py
@@ -225,7 +225,7 @@ def get_spectrum(self, nu_range=None):
         if self.spectrum is not None:
             nu, fnu = self.spectrum['nu'], self.spectrum['fnu']
             if nu_range is not None:
-                raise NotImplemented("nu_range not yet implemented for spectrum")
+                raise NotImplementedError("nu_range not yet implemented for spectrum")
         elif self.temperature is not None:
             if nu_range is None:
                 raise ValueError("nu_range is needed for sources with Planck spectra")
diff --git a/hyperion/sources/tests/test_source.py b/hyperion/sources/tests/test_source.py
index 379a7afd..26f675aa 100644
--- a/hyperion/sources/tests/test_source.py
+++ b/hyperion/sources/tests/test_source.py
@@ -771,3 +771,11 @@ def test_map_write_invalid2():
     with pytest.raises(Exception) as exc:
         s.write(v, 'test', None)
     assert exc.value.args[0] == 'map is zero everywhere'
+
+
+def test_spectrum_nu_range_not_implemented():
+    s = Source()
+    s.luminosity = 1.
+    s.spectrum = (np.array([1., 2., 3.]), np.array([1., 1., 1.]))
+    with pytest.raises(NotImplementedError):
+        s.get_spectrum(nu_range=(1., 2.))
diff --git a/hyperion/util/convenience.py b/hyperion/util/convenience.py
index aeb0421e..0772d50f 100644
--- a/hyperion/util/convenience.py
+++ b/hyperion/util/convenience.py
@@ -8,13 +8,15 @@ class OptThinRadius(object):
     def __init__(self, temperature, value=1., min=0.):
         self.temperature = temperature
         self.value = value
-        self.min = 0.
+        self.min = min
 
     def __mul__(self, value):
-        return OptThinRadius(self.temperature, value=self.value * value)
+        return OptThinRadius(self.temperature, value=self.value * value,
+                             min=self.min)
 
     def __rmul__(self, value):
-        return OptThinRadius(self.temperature, value=self.value * value)
+        return OptThinRadius(self.temperature, value=self.value * value,
+                             min=self.min)
 
     def __str__(self):
         return "%g times the radius at which the optically thin temperature would be %gK" % (self.value, self.temperature)
diff --git a/hyperion/util/integrate.py b/hyperion/util/integrate.py
index 757fa411..19a3b6d0 100644
--- a/hyperion/util/integrate.py
+++ b/hyperion/util/integrate.py
@@ -44,7 +44,7 @@ def integrate_subset(x, y, xmin, xmax):
         ymin = interp1d_fast(x[i1 - 1:i1 + 1], y[i1 - 1:i1 + 1], xmin)
 
     if xmax == x[-1]:
-        i2 = -2
+        i2 = -1
         ymax = y[-1]
     else:
         i2 = np.searchsorted(x, xmax)
@@ -90,7 +90,7 @@ def integrate_loglin_subset(x, y, xmin, xmax):
         ymin = interp1d_fast_loglin(x[i1 - 1:i1 + 1], y[i1 - 1:i1 + 1], xmin)
 
     if xmax == x[-1]:
-        i2 = -2
+        i2 = -1
         ymax = y[-1]
     else:
         i2 = np.searchsorted(x, xmax)
@@ -136,7 +136,7 @@ def integrate_linlog_subset(x, y, xmin, xmax):
         ymin = interp1d_fast_linlog(x[i1 - 1:i1 + 1], y[i1 - 1:i1 + 1], xmin)
 
     if xmax == x[-1]:
-        i2 = -2
+        i2 = -1
         ymax = y[-1]
     else:
         i2 = np.searchsorted(x, xmax)
@@ -181,7 +181,7 @@ def integrate_loglog_subset(x, y, xmin, xmax):
         ymin = interp1d_fast_loglog(x[i1 - 1:i1 + 1], y[i1 - 1:i1 + 1], xmin)
 
     if xmax == x[-1]:
-        i2 = -2
+        i2 = -1
         ymax = y[-1]
     else:
         i2 = np.searchsorted(x, xmax)
diff --git a/hyperion/util/tests/test_convenience.py b/hyperion/util/tests/test_convenience.py
new file mode 100644
index 00000000..4285572c
--- /dev/null
+++ b/hyperion/util/tests/test_convenience.py
@@ -0,0 +1,15 @@
+from __future__ import print_function, division
+
+from ..convenience import OptThinRadius
+
+
+def test_optthin_radius_min_stored():
+    r = OptThinRadius(1600., min=3.)
+    assert r.temperature == 1600.
+    assert r.min == 3.
+
+
+def test_optthin_radius_min_propagated():
+    r = OptThinRadius(1600., min=3.)
+    assert (r * 2.).min == 3.
+    assert (2. * r).min == 3.
diff --git a/hyperion/util/tests/test_integrate.py b/hyperion/util/tests/test_integrate.py
index 085ae393..b3f8b715 100644
--- a/hyperion/util/tests/test_integrate.py
+++ b/hyperion/util/tests/test_integrate.py
@@ -160,3 +160,16 @@ def test_linlog_not_monotonic():
         assert exc.value == 'x is not monotonically increasing'
     else:
         assert exc.value.args[0] == 'x is not monotonically increasing'
+
+
+@pytest.mark.parametrize(('subset', 'full'),
+                         [(integrate_subset, integrate),
+                          (integrate_loglin_subset, integrate_loglin),
+                          (integrate_linlog_subset, integrate_linlog),
+                          (integrate_loglog_subset, integrate_loglog)])
+def test_subset_full_range(subset, full):
+    # Regression test for a bug that caused the second-to-last tabulated
+    # point to be dropped when the upper limit was the last tabulated value
+    x = np.array([1., 2., 3., 4., 5.])
+    y = np.array([1., 3., 2., 5., 4.])
+    assert almost_equal(subset(x, y, x[0], x[-1]), full(x, y))
diff --git a/src/grid/grid_geometry_amr.f90 b/src/grid/grid_geometry_amr.f90
index 94aa4997..30eaacde 100644
--- a/src/grid/grid_geometry_amr.f90
+++ b/src/grid/grid_geometry_amr.f90
@@ -709,7 +709,7 @@ logical function in_correct_cell(p)
           in_correct_cell = in_correct_cell .and. icell_actual%i3 == p%icell%i3
        end if
     else
-       in_correct_cell = icell_actual == p%icell
+       in_correct_cell = icell_actual%i1 == p%icell%i1 .and. icell_actual%i2 == p%icell%i2 .and. icell_actual%i3 == p%icell%i3
     end if
   end function in_correct_cell
 
diff --git a/src/grid/grid_geometry_common_3d.f90 b/src/grid/grid_geometry_common_3d.f90
index 0175b98e..983b1c74 100644
--- a/src/grid/grid_geometry_common_3d.f90
+++ b/src/grid/grid_geometry_common_3d.f90
@@ -81,7 +81,7 @@ subroutine grid_sample_pdf_map2(pdf, icell, prob)
 
     do
        call random(xi)
-       ic = ceiling(xi * pdf%n)
+       ic = max(ceiling(xi * pdf%n), 1)
        prob = pdf%pdf(ic)
        if(prob > 1.e-100_dp) exit
     end do
@@ -97,7 +97,7 @@ type(grid_cell) function random_cell()
     real(dp) :: xi
     integer :: ic
     call random(xi)
-    ic = ceiling(xi * geo%n_cells)
+    ic = max(ceiling(xi * geo%n_cells), 1)
     random_cell = new_grid_cell(ic, geo)
   end function random_cell
 
@@ -107,9 +107,9 @@ type(grid_cell) function random_masked_cell()
     integer :: ic
     call random(xi)
     if (geo%masked) then
-       ic = geo%mask_map(ceiling(xi * geo%n_masked))
+       ic = geo%mask_map(max(ceiling(xi * geo%n_masked), 1))
     else
-       ic = ceiling(xi * geo%n_cells)
+       ic = max(ceiling(xi * geo%n_cells), 1)
     end if
     random_masked_cell = new_grid_cell(ic, geo)
   end function random_masked_cell
diff --git a/src/grid/grid_geometry_octree.f90 b/src/grid/grid_geometry_octree.f90
index 696dda86..a0e09709 100644
--- a/src/grid/grid_geometry_octree.f90
+++ b/src/grid/grid_geometry_octree.f90
@@ -279,7 +279,7 @@ type(grid_cell) function find_cell(p) result(icell)
     type(photon),intent(in) :: p
     integer :: ic
     if(debug) write(*,'(" [debug] find_cell")')
-    if(p%r%x<geo%xmin.or.p%r%y>geo%xmax) then
+    if(p%r%x<geo%xmin.or.p%r%x>geo%xmax) then
        call warn("find_cell","photon not in grid (in x direction)")
        icell = invalid_cell
        return
@@ -289,7 +289,7 @@ type(grid_cell) function find_cell(p) result(icell)
        icell = invalid_cell
        return
     end if
-    if(p%r%z<geo%zmin.or.p%r%y>geo%zmax) then
+    if(p%r%z<geo%zmin.or.p%r%z>geo%zmax) then
        call warn("find_cell","photon not in grid (in z direction)")
        icell = invalid_cell
        return
@@ -415,11 +415,11 @@ real(dp) function distance_to_closest_wall(p) result(d)
 
     real(dp) :: d1,d2,d3,d4,d5,d6
 
-    d1 = p%r%x - geo%cells(p%icell%ic)%x - geo%cells(p%icell%ic)%dx
+    d1 = p%r%x - geo%cells(p%icell%ic)%x + geo%cells(p%icell%ic)%dx
     d2 = geo%cells(p%icell%ic)%x + geo%cells(p%icell%ic)%dx - p%r%x
-    d3 = p%r%y - geo%cells(p%icell%ic)%y - geo%cells(p%icell%ic)%dy
+    d3 = p%r%y - geo%cells(p%icell%ic)%y + geo%cells(p%icell%ic)%dy
     d4 = geo%cells(p%icell%ic)%y + geo%cells(p%icell%ic)%dy - p%r%y
-    d5 = p%r%z - geo%cells(p%icell%ic)%z - geo%cells(p%icell%ic)%dz
+    d5 = p%r%z - geo%cells(p%icell%ic)%z + geo%cells(p%icell%ic)%dz
     d6 = geo%cells(p%icell%ic)%z + geo%cells(p%icell%ic)%dz - p%r%z
 
     ! Find the smallest of the distances
diff --git a/src/grid/grid_io_1d.f90 b/src/grid/grid_io_1d.f90
index 7388046d..0fd08086 100644
--- a/src/grid/grid_io_1d.f90
+++ b/src/grid/grid_io_1d.f90
@@ -1,4 +1,4 @@
-! MD5 of template: cd43feff40838596c8e4acdfbee7311d
+! MD5 of template: 3e49d4dccf183e62f016ef07b2c97b4e
 module grid_io
 
   use core_lib
@@ -15,7 +15,6 @@ module grid_io
   public :: write_grid_3d
   public :: write_grid_4d
   public :: write_grid_5d
-  
 
   interface read_grid_3d
      module procedure read_grid_3d_sp
@@ -44,7 +43,7 @@ module grid_io
      module procedure write_grid_4d_int
      module procedure write_grid_4d_int8
   end interface write_grid_4d
-
+  
   interface write_grid_5d
      module procedure write_grid_5d_sp
      module procedure write_grid_5d_dp
@@ -53,6 +52,7 @@ module grid_io
   end interface write_grid_5d
 
 
+
 contains
 
   logical function grid_exists(group, name)
@@ -63,6 +63,7 @@ logical function grid_exists(group, name)
   end function grid_exists
 
 
+
   subroutine read_grid_4d_int8(group, path, array, geo)
 
     implicit none
@@ -106,11 +107,10 @@ subroutine read_grid_3d_int8(group, path, array, geo)
   end subroutine read_grid_3d_int8
 
 
-
   subroutine write_grid_5d_int8(group, path, array, geo)
-    
+
     implicit none
-    
+
     integer(hid_t), intent(in) :: group
     character(len=*), intent(in) :: path
     integer(idp), intent(in) :: array(:,:,:)
@@ -151,6 +151,7 @@ subroutine write_grid_3d_int8(group, path, array, geo)
   end subroutine write_grid_3d_int8
 
 
+
   subroutine read_grid_4d_int(group, path, array, geo)
 
     implicit none
@@ -194,9 +195,6 @@ subroutine read_grid_3d_int(group, path, array, geo)
   end subroutine read_grid_3d_int
 
 
-
-
-
   subroutine write_grid_5d_int(group, path, array, geo)
 
     implicit none
@@ -240,7 +238,8 @@ subroutine write_grid_3d_int(group, path, array, geo)
 
   end subroutine write_grid_3d_int
 
- 
+
+
   subroutine read_grid_4d_dp(group, path, array, geo)
 
     implicit none
@@ -284,9 +283,8 @@ subroutine read_grid_3d_dp(group, path, array, geo)
   end subroutine read_grid_3d_dp
 
 
-
   subroutine write_grid_5d_dp(group, path, array, geo)
-    
+
     implicit none
 
     integer(hid_t), intent(in) :: group
@@ -328,7 +326,8 @@ subroutine write_grid_3d_dp(group, path, array, geo)
 
   end subroutine write_grid_3d_dp
 
-  
+
+
   subroutine read_grid_4d_sp(group, path, array, geo)
 
     implicit none
@@ -372,19 +371,18 @@ subroutine read_grid_3d_sp(group, path, array, geo)
   end subroutine read_grid_3d_sp
 
 
-
   subroutine write_grid_5d_sp(group, path, array, geo)
-    
+
     implicit none
-    
+
     integer(hid_t), intent(in) :: group
     character(len=*), intent(in) :: path
     real(sp), intent(in) :: array(:,:,:)
     type(grid_geometry_desc),intent(in) :: geo
-    
+
     call mp_write_array(group, path, array)
     call mp_write_keyword(group, path, 'geometry', geo%id)
-    
+
   end subroutine write_grid_5d_sp
 
 
diff --git a/src/grid/grid_io_1d_template.f90 b/src/grid/grid_io_1d_template.f90
index f36aa0b7..1f10e420 100644
--- a/src/grid/grid_io_1d_template.f90
+++ b/src/grid/grid_io_1d_template.f90
@@ -13,6 +13,7 @@ module grid_io
   public :: read_grid_4d
   public :: write_grid_3d
   public :: write_grid_4d
+  public :: write_grid_5d
 
   interface read_grid_3d
      module procedure read_grid_3d_sp
diff --git a/src/grid/grid_io_amr.f90 b/src/grid/grid_io_amr.f90
index b49281c4..a950e188 100644
--- a/src/grid/grid_io_amr.f90
+++ b/src/grid/grid_io_amr.f90
@@ -1,4 +1,4 @@
-! MD5 of template: a9d1b24b38764cb05db2f7aeda80e7ca
+! MD5 of template: fe55a2c227b277e3b81bf41dc7696fa9
 module grid_io
 
   use core_lib
@@ -136,7 +136,7 @@ subroutine read_grid_3d_int8(group, path, array, geo)
        level => geo%levels(ilevel)
        do igrid=1,size(level%grids)
           grid => level%grids(igrid)
-          write(full_path, '("level_", I5.5, "/grid_ ", I5.5,"/")') ilevel, igrid
+          write(full_path, '("level_", I5.5, "/grid_", I5.5,"/")') ilevel, igrid
           full_path = trim(full_path)//trim(path)
           call mp_read_array_auto(group, full_path, array3d)
           if(any(is_nan(array3d))) call error("read_grid_3d", "NaN values in 3D array")
@@ -341,7 +341,7 @@ subroutine read_grid_3d_int(group, path, array, geo)
        level => geo%levels(ilevel)
        do igrid=1,size(level%grids)
           grid => level%grids(igrid)
-          write(full_path, '("level_", I5.5, "/grid_ ", I5.5,"/")') ilevel, igrid
+          write(full_path, '("level_", I5.5, "/grid_", I5.5,"/")') ilevel, igrid
           full_path = trim(full_path)//trim(path)
           call mp_read_array_auto(group, full_path, array3d)
           if(any(is_nan(array3d))) call error("read_grid_3d", "NaN values in 3D array")
@@ -540,7 +540,7 @@ subroutine read_grid_3d_dp(group, path, array, geo)
        level => geo%levels(ilevel)
        do igrid=1,size(level%grids)
           grid => level%grids(igrid)
-          write(full_path, '("level_", I5.5, "/grid_ ", I5.5,"/")') ilevel, igrid
+          write(full_path, '("level_", I5.5, "/grid_", I5.5,"/")') ilevel, igrid
           full_path = trim(full_path)//trim(path)
           call mp_read_array_auto(group, full_path, array3d)
           if(any(is_nan(array3d))) call error("read_grid_3d", "NaN values in 3D array")
@@ -740,7 +740,7 @@ subroutine read_grid_3d_sp(group, path, array, geo)
        level => geo%levels(ilevel)
        do igrid=1,size(level%grids)
           grid => level%grids(igrid)
-          write(full_path, '("level_", I5.5, "/grid_ ", I5.5,"/")') ilevel, igrid
+          write(full_path, '("level_", I5.5, "/grid_", I5.5,"/")') ilevel, igrid
           full_path = trim(full_path)//trim(path)
           call mp_read_array_auto(group, full_path, array3d)
           if(any(is_nan(array3d))) call error("read_grid_3d", "NaN values in 3D array")
diff --git a/src/grid/grid_io_amr_template.f90 b/src/grid/grid_io_amr_template.f90
index a37d6354..8799b88c 100644
--- a/src/grid/grid_io_amr_template.f90
+++ b/src/grid/grid_io_amr_template.f90
@@ -137,7 +137,7 @@ subroutine read_grid_3d_<T>(group, path, array, geo)
        level => geo%levels(ilevel)
        do igrid=1,size(level%grids)
           grid => level%grids(igrid)
-          write(full_path, '("level_", I5.5, "/grid_ ", I5.5,"/")') ilevel, igrid
+          write(full_path, '("level_", I5.5, "/grid_", I5.5,"/")') ilevel, igrid
           full_path = trim(full_path)//trim(path)
           call mp_read_array_auto(group, full_path, array3d)
           if(any(is_nan(array3d))) call error("read_grid_3d", "NaN values in 3D array")
diff --git a/src/grid/grid_monochromatic.f90 b/src/grid/grid_monochromatic.f90
index d5846d98..8f526fe0 100644
--- a/src/grid/grid_monochromatic.f90
+++ b/src/grid/grid_monochromatic.f90
@@ -142,7 +142,7 @@ type(photon) function emit_from_monochromatic_grid_pdf(inu) result(p)
     call update_optconsts(p)
 
     call random(xi)
-    dust_id = ceiling(xi*real(n_dust, dp))
+    dust_id = max(ceiling(xi*real(n_dust, dp)), 1)
 
     if(mean_prob(dust_id) == 0._dp) then
         p%energy = 0._dp
diff --git a/src/grid/grid_pda_3d.f90 b/src/grid/grid_pda_3d.f90
index 69ec9b98..dbe237b9 100644
--- a/src/grid/grid_pda_3d.f90
+++ b/src/grid/grid_pda_3d.f90
@@ -227,7 +227,7 @@ subroutine solve_pda_indiv_exact(pda_cells, id_pda_cell)
 
           if(id_pda_cell(next%ic) > 0) then
              id_next = id_pda_cell(next%ic)
-             a(id_next, id_curr) = coefficient
+             a(id_next, id_curr) = a(id_next, id_curr) + coefficient
           else
              b(id_curr) = b(id_curr) - coefficient * e_mean(next%ic)
           end if
diff --git a/src/grid/grid_physics_3d.f90 b/src/grid/grid_physics_3d.f90
index 36184344..ffc014a4 100644
--- a/src/grid/grid_physics_3d.f90
+++ b/src/grid/grid_physics_3d.f90
@@ -352,10 +352,10 @@ subroutine sublimate_dust()
     n_nu_bins = 0
     if (compute_specific_energy_spectrum) n_nu_bins = size(specific_energy_spectrum, 3)
 
-    reset = 0
-
     do id=1,n_dust
 
+       reset = 0
+
        select case(d(id)%sublimation_mode)
        case(1)
 
@@ -447,10 +447,13 @@ subroutine update_energy_abs(scale)
        if (compute_specific_energy_spectrum) then
           do idx=1,n_nu_bins
              specific_energy_spectrum(:,id,idx) = specific_energy_sum_spectrum(:,id,idx) * scale/geo%volume
+             where(geo%volume == 0._dp)
+                specific_energy_spectrum(:,id,idx) = 0._dp
+             end where
           end do
        end if
-          
-          
+
+
        where(geo%volume == 0._dp)
           specific_energy(:,id) = 0._dp
        end where
@@ -636,7 +639,7 @@ type(photon) function emit_from_grid(inu) result(p)
     end if
 
     call random(xi)
-    p%dust_id = ceiling(xi*real(n_dust,dp))
+    p%dust_id = max(ceiling(xi*real(n_dust,dp)), 1)
 
     ! Pick random cell
     p%icell = random_masked_cell()
diff --git a/src/main/counters.f90 b/src/main/counters.f90
index cbb22899..03428df4 100644
--- a/src/main/counters.f90
+++ b/src/main/counters.f90
@@ -1,11 +1,13 @@
 module counters
 
+  use core_lib, only : idp
+
   implicit none
   save
 
-  integer :: killed_photons_geo = 0
-  integer :: killed_photons_int = 0
-  integer :: killed_photons_reabs = 0
+  integer(idp) :: killed_photons_geo = 0
+  integer(idp) :: killed_photons_int = 0
+  integer(idp) :: killed_photons_reabs = 0
 
 end module counters
 
diff --git a/src/main/main.f90 b/src/main/main.f90
index e7f451ec..f8975462 100644
--- a/src/main/main.f90
+++ b/src/main/main.f90
@@ -169,6 +169,9 @@ program main
   ! Wait for all threads
   call mp_join()
 
+  ! Initialize convergence flag in case convergence checking is disabled
+  converged = .false.
+
   ! Loop over Lucy iterations
   do iter=1,n_initial_iter
 
diff --git a/src/mpi/mpi_io.f90 b/src/mpi/mpi_io.f90
index 002efbf5..f5226285 100644
--- a/src/mpi/mpi_io.f90
+++ b/src/mpi/mpi_io.f90
@@ -1,4 +1,4 @@
-! MD5 of template: 1e1c708fb2055feec7cad222682bbe09
+! MD5 of template: 744398b07cde70659e77bf8ee3986c3f
 module mpi_hdf5_io
 
   use core_lib
@@ -369,7 +369,7 @@ subroutine mp_table_read_column_auto_1d_mpi_character(handle, path, name, array)
     n1 = size(array)
     call mpi_bcast(n1, 1, mpi_integer4, rank_main, mpi_comm_world, ierr)
     if(.not. main_process()) allocate(array(n1))
-    call mpi_bcast(array, n1, mpi_character, rank_main, mpi_comm_world, ierr)
+    call mpi_bcast(array, n1 * len(array), mpi_character, rank_main, mpi_comm_world, ierr)
   end subroutine mp_table_read_column_auto_1d_mpi_character
 
   subroutine mp_table_write_column_1d_mpi_character(handle, path, name, array)
diff --git a/src/mpi/mpi_io_template.f90 b/src/mpi/mpi_io_template.f90
index 59189739..2dd83dac 100644
--- a/src/mpi/mpi_io_template.f90
+++ b/src/mpi/mpi_io_template.f90
@@ -368,7 +368,7 @@ subroutine mp_table_read_column_auto_1d_mpi_character(handle, path, name, array)
     n1 = size(array)
     call mpi_bcast(n1, 1, mpi_integer4, rank_main, mpi_comm_world, ierr)
     if(.not. main_process()) allocate(array(n1))
-    call mpi_bcast(array, n1, mpi_character, rank_main, mpi_comm_world, ierr)
+    call mpi_bcast(array, n1 * len(array), mpi_character, rank_main, mpi_comm_world, ierr)
   end subroutine mp_table_read_column_auto_1d_mpi_character
 
   subroutine mp_table_write_column_1d_mpi_character(handle, path, name, array)
diff --git a/src/mpi/mpi_routines.f90 b/src/mpi/mpi_routines.f90
index a131a93b..92fdbb94 100644
--- a/src/mpi/mpi_routines.f90
+++ b/src/mpi/mpi_routines.f90
@@ -92,6 +92,7 @@ subroutine mp_n_photons(n_photons_tot, n_photons_curr, n_photons_stats, n_photon
 
     real(dp),save :: time1 = -1._dp
     real(dp) :: time2
+    real(dp) :: dtime_send
     real(dp),allocatable,volatile,save :: dtime(:)
 
     if(debug) write(*,'("[mpi_routines] rank ",I0," requesting number of photons")') rank
@@ -243,7 +244,8 @@ subroutine mp_n_photons(n_photons_tot, n_photons_curr, n_photons_stats, n_photon
 
        ! Receive number of photons and send acknowledgments
        call mpi_recv(n_photons, 1, mpi_integer8, rank_main, tag1, mpi_comm_world, status, ierr)
-       call mpi_isend(time2-time1, 1, mpi_real8, rank_main, tag2, mpi_comm_world, request_dum, ierr)
+       dtime_send = time2 - time1
+       call mpi_send(dtime_send, 1, mpi_real8, rank_main, tag2, mpi_comm_world, ierr)
 
        if(n_photons > n_photons_tot) stop "n_photons > n_photons_tot"