From 04409bf04a5aab6931ac5606aa4f7898097a92b5 Mon Sep 17 00:00:00 2001
From: Thomas Robitaille <thomas.robitaille@gmail.com>
Date: Sat, 1 Aug 2026 22:34:08 +0000
Subject: [PATCH 1/3] Make the test suite safe to run in parallel with
 pytest-xdist by writing temporary model files to tmpdir and making
 parametrized test collection deterministic

---
 hyperion/model/tests/test_model.py    |  8 ++++----
 hyperion/sources/tests/test_source.py | 10 +++++-----
 2 files changed, 9 insertions(+), 9 deletions(-)

diff --git a/hyperion/model/tests/test_model.py b/hyperion/model/tests/test_model.py
index 967d3cb4..9d192776 100644
--- a/hyperion/model/tests/test_model.py
+++ b/hyperion/model/tests/test_model.py
@@ -32,18 +32,18 @@ def test_noname_nofilename():
     assert e.value.args[0] == "filename= has not been specified and model has no name"
 
 
-def test_nogrid():
+def test_nogrid(tmpdir):
     m = Model()
     with pytest.raises(Exception) as e:
-        m.write('test')
+        m.write(tmpdir.join('test.rtin').strpath)
     assert e.value.args[0] == 'No coordinate grid has been set up'
 
 
-def test_nophotons():
+def test_nophotons(tmpdir):
     m = Model()
     m.set_cartesian_grid([-1., 1.], [-1., 1.], [-1., 1.])
     with pytest.raises(Exception) as e:
-        m.write('test')
+        m.write(tmpdir.join('test.rtin').strpath)
     assert e.value.args[0] == 'Photon numbers not set'
 
 
diff --git a/hyperion/sources/tests/test_source.py b/hyperion/sources/tests/test_source.py
index 26f675aa..bc42710d 100644
--- a/hyperion/sources/tests/test_source.py
+++ b/hyperion/sources/tests/test_source.py
@@ -24,13 +24,13 @@
 # SCALAR LUMINOSITY
 
 
-@pytest.mark.parametrize(('source_type'), list(set(ALL_SOURCES) - set([PointSourceCollection])))
+@pytest.mark.parametrize(('source_type'), [s for s in ALL_SOURCES if s is not PointSourceCollection])
 def test_luminosity_scalar(source_type):
     s = source_type()
     s.luminosity = 1.
 
 
-@pytest.mark.parametrize(('source_type'), list(set(ALL_SOURCES) - set([PointSourceCollection])))
+@pytest.mark.parametrize(('source_type'), [s for s in ALL_SOURCES if s is not PointSourceCollection])
 def test_luminosity_scalar_invalid2(source_type):
     s = source_type()
     with pytest.raises(ValueError) as exc:
@@ -38,7 +38,7 @@ def test_luminosity_scalar_invalid2(source_type):
     assert exc.value.args[0] == 'luminosity should be a scalar value'
 
 
-@pytest.mark.parametrize(('source_type'), list(set(ALL_SOURCES) - set([PointSourceCollection])))
+@pytest.mark.parametrize(('source_type'), [s for s in ALL_SOURCES if s is not PointSourceCollection])
 def test_luminosity_scalar_invalid3(source_type):
     s = source_type()
     with pytest.raises(ValueError) as exc:
@@ -46,7 +46,7 @@ def test_luminosity_scalar_invalid3(source_type):
     assert exc.value.args[0] == 'luminosity should be a numerical value'
 
 
-@pytest.mark.parametrize(('source_type'), list(set(ALL_SOURCES) - set([PointSourceCollection])))
+@pytest.mark.parametrize(('source_type'), [s for s in ALL_SOURCES if s is not PointSourceCollection])
 def test_luminosity_scalar_invalid4(source_type):
     s = source_type()
     with pytest.raises(ValueError) as exc:
@@ -88,7 +88,7 @@ def test_luminosity_array_invalid3(source_type):
 # TEMPERATURE
 
 
-@pytest.mark.parametrize(('source_type'), list(set(ALL_SOURCES) - set([MapSource])))
+@pytest.mark.parametrize(('source_type'), [s for s in ALL_SOURCES if s is not MapSource])
 def test_temperature(source_type):
     v = virtual_file()
     s = source_type()

From 4559b6851ecdb7aa3aac29bdf4169f7e1f61ce25 Mon Sep 17 00:00:00 2001
From: Thomas Robitaille <thomas.robitaille@gmail.com>
Date: Sat, 1 Aug 2026 22:34:08 +0000
Subject: [PATCH 2/3] Bin the frequency-resolved specific energy onto
 user-specified frequency bin edges

---
 docs/setup/setup_conf.rst                     |  44 +++---
 hyperion/conf/conf_files.py                   |  73 ++++++----
 hyperion/conf/tests/test_conf_io.py           |  19 ++-
 hyperion/grid/cartesian_grid.py               |   2 +-
 hyperion/grid/cylindrical_polar_grid.py       |   2 +-
 hyperion/grid/octree_grid.py                  |   2 +-
 hyperion/grid/spherical_polar_grid.py         |   2 +-
 hyperion/grid/voronoi_grid.py                 |   2 +-
 .../tests/test_specific_energy_spectrum.py    | 126 +++++++++++++-----
 hyperion/util/functions.py                    |   1 +
 src/dust/dust_type_4elem.f90                  |   6 +-
 src/grid/grid_generic.f90                     |   8 +-
 src/grid/grid_physics_3d.f90                  |  69 ++++------
 src/grid/grid_propagate_3d.f90                |  10 +-
 src/main/settings.f90                         |   7 +-
 src/main/setup_rt.f90                         |   9 +-
 16 files changed, 223 insertions(+), 159 deletions(-)

diff --git a/docs/setup/setup_conf.rst b/docs/setup/setup_conf.rst
index 21c7747c..046a4107 100644
--- a/docs/setup/setup_conf.rst
+++ b/docs/setup/setup_conf.rst
@@ -207,40 +207,34 @@ quantities::
                                               # 'last'  -> only final iteration
                                               # 'none'  -> not computed or saved (default)
 
-When it is not ``'none'``, two extra arrays are written out:
+When it is not ``'none'``, the frequency bins have to be set explicitly with::
+
+    import numpy as np
+    m.set_specific_energy_spectrum_bins(np.logspace(11., 16., 101))
+
+which takes the ``n + 1`` bin edges (in Hz, in increasing order) that define
+``n`` bins. Two extra datasets are then written out:
 
 * ``specific_energy_spectrum`` -- the specific energy absorbed in each cell as a
-  function of frequency, in erg/s/g.
-* ``specific_energy_spectrum_frequencies`` -- the frequencies (in Hz) corresponding to
-  the leading axis of ``specific_energy_spectrum``. These are bin centers, not edges,
-  so this array has exactly the same length as that axis (one entry per bin).
+  function of frequency, in erg/s/g, with one entry per bin.
+* ``specific_energy_spectrum_bin_edges`` -- the edges (in Hz) of the frequency
+  bins, with one more entry than the number of bins.
+
+Since the spectrum is a histogram, it is best plotted against the bin edges as
+steps (e.g. with ``matplotlib``'s ``stairs``) rather than at a single
+representative frequency per bin.
 
 This is the *absorbed* (deposited) energy spectrum, not the mean intensity: each
 contribution is weighted by the dust absorption opacity, so it is proportional
 to :math:`\kappa_\nu J_\nu`. To recover the mean intensity :math:`J_\nu` (the
 radiation field), divide by the dust absorption opacity at each frequency.
-Summed over frequency, ``specific_energy_spectrum`` recovers the total
-``specific_energy``.
 
-By default the binning uses the frequency grid of the first dust type. You can
-instead provide your own frequency grid (in Hz), independent of the dust
-properties::
+Energy absorbed from photons with frequencies outside the outermost edges is
+*not* included in the spectrum, so summed over frequency,
+``specific_energy_spectrum`` only recovers the total ``specific_energy`` if the
+edges span the full range of frequencies over which energy is absorbed.
 
-    import numpy as np
-    m.set_specific_energy_spectrum_frequencies(np.logspace(11., 16., 100))
-
-Photons are binned to the nearest of these frequencies in log space (so the
-supplied values act as bin centers). This works for all
-grid types, including AMR and Voronoi.
-
-.. warning:: The binning has no outer edges: the first and last bins collect
-             *all* photons below and above the outermost bin centers
-             respectively, however far away in frequency. This keeps the spectrum
-             energy-conserving (no photons are dropped), but it means a grid that
-             does not span the full frequency range will pile up out-of-range
-             flux in its end bins. When supplying a custom grid, make sure it
-             brackets the full range of frequencies present in the simulation
-             (the default dust-based grid already does this).
+This works for all grid types, including AMR and Voronoi.
 
 ``specific_energy_spectrum`` can be retrieved like other grid quantities, as an array
 with an extra leading frequency axis::
diff --git a/hyperion/conf/conf_files.py b/hyperion/conf/conf_files.py
index bcc73bd0..0413a18c 100644
--- a/hyperion/conf/conf_files.py
+++ b/hyperion/conf/conf_files.py
@@ -58,7 +58,7 @@ def __init__(self):
         self.set_max_reabsorptions(1000000)
         self.set_pda(False)
         self.set_mrw(False)
-        self.specific_energy_spectrum_frequencies = None
+        self.specific_energy_spectrum_bin_edges = None
 
         self.set_convergence(False)
         self.set_kill_on_absorb(False)
@@ -399,41 +399,58 @@ def _read_pda(self, group):
     def _write_pda(self, group):
         group.attrs['pda'] = bool2str(self.pda)
 
-    def _read_specific_energy_spectrum_frequencies(self, group):
-        if 'specific_energy_spectrum_frequencies' in group:
-            self.specific_energy_spectrum_frequencies = np.array(group['specific_energy_spectrum_frequencies']['nu'])
+    def _read_specific_energy_spectrum_bins(self, group):
+        if 'specific_energy_spectrum_bin_edges' in group:
+            self.specific_energy_spectrum_bin_edges = \
+                np.array(group['specific_energy_spectrum_bin_edges']['nu'])
         else:
-            self.specific_energy_spectrum_frequencies = None
+            self.specific_energy_spectrum_bin_edges = None
 
-    def _write_specific_energy_spectrum_frequencies(self, group):
-        if self.specific_energy_spectrum_frequencies is not None:
-            group.create_dataset('specific_energy_spectrum_frequencies',
-                                 data=np.array(list(zip(self.specific_energy_spectrum_frequencies)),
-                                               dtype=[('nu', float)]))
+    def _write_specific_energy_spectrum_bins(self, group):
 
-    def set_specific_energy_spectrum_frequencies(self, frequencies):
+        edges = self.specific_energy_spectrum_bin_edges
+
+        if edges is None:
+            conf = getattr(self, 'conf', None)
+            if conf is not None and conf.output.output_specific_energy_spectrum != 'none':
+                raise ValueError("output_specific_energy_spectrum is enabled but the "
+                                 "frequency bins have not been set - use "
+                                 "set_specific_energy_spectrum_bins to set them")
+            return
+
+        group.create_dataset('specific_energy_spectrum_bin_edges',
+                             data=np.array(list(zip(edges)), dtype=[('nu', float)]))
+
+    def set_specific_energy_spectrum_bins(self, edges):
 
         '''
-        Set the frequency grid onto which the frequency-resolved specific energy
-        (``specific_energy_spectrum``) is binned.
+        Set the frequency bins onto which the frequency-resolved specific
+        energy (``specific_energy_spectrum``) is binned.
 
-        This is only relevant if ``conf.output.output_specific_energy_spectrum`` is set
-        to ``'all'`` or ``'last'``. If this method is not called, the frequency
-        grid of the first dust type is used. Photons are binned to the nearest
-        frequency in log space, so the supplied values are bin centers (not
-        edges) and ``specific_energy_spectrum`` has one entry per supplied frequency.
+        This is required if ``conf.output.output_specific_energy_spectrum`` is
+        set to ``'all'`` or ``'last'``.
 
         Parameters
         ----------
-        frequencies : iterable of float
-            The frequencies (in Hz) onto which to bin ``specific_energy_spectrum``.
+        edges : iterable of float
+            The edges of the frequency bins (in Hz), in increasing order:
+            ``n + 1`` values define ``n`` bins, and
+            ``specific_energy_spectrum`` has one entry per bin. Energy
+            absorbed from photons with frequencies outside the outer edges
+            is not included in the spectrum, so the spectrum summed over the
+            bins only recovers ``specific_energy`` if the edges span the
+            frequency range over which energy is absorbed.
         '''
-        frequencies = np.asarray(frequencies, dtype=float)
-        if frequencies.ndim != 1 or frequencies.size < 1:
-            raise ValueError("frequencies should be a 1-d array of at least one value")
-        if np.any(frequencies <= 0.):
-            raise ValueError("frequencies should be positive (in Hz)")
-        self.specific_energy_spectrum_frequencies = np.sort(frequencies)
+
+        edges = np.asarray(edges, dtype=float)
+        if edges.ndim != 1 or edges.size < 2:
+            raise ValueError("edges should be a 1-d array of at least two values")
+        if np.any(edges <= 0.):
+            raise ValueError("edges should be positive (in Hz)")
+        if np.any(np.diff(edges) <= 0.):
+            raise ValueError("edges should be strictly increasing")
+
+        self.specific_energy_spectrum_bin_edges = edges
 
 
     def set_mrw(self, mrw, gamma=1.0, inter_max=1000, warn=True):
@@ -764,7 +781,7 @@ def read_run_conf(self, group):  # not a class method because inherited
         self._read_max_reabsorptions(group)
         self._read_pda(group)
         self._read_mrw(group)
-        self._read_specific_energy_spectrum_frequencies(group)
+        self._read_specific_energy_spectrum_bins(group)
         self._read_convergence(group)
         self._read_kill_on_absorb(group)
         self._read_kill_on_scatter(group)
@@ -793,7 +810,7 @@ def write_run_conf(self, group):
         self._write_max_reabsorptions(group)
         self._write_pda(group)
         self._write_mrw(group)
-        self._write_specific_energy_spectrum_frequencies(group)
+        self._write_specific_energy_spectrum_bins(group)
         self._write_convergence(group)
         self._write_kill_on_absorb(group)
         self._write_kill_on_scatter(group)
diff --git a/hyperion/conf/tests/test_conf_io.py b/hyperion/conf/tests/test_conf_io.py
index 9cfe7ed5..e4a5dc10 100644
--- a/hyperion/conf/tests/test_conf_io.py
+++ b/hyperion/conf/tests/test_conf_io.py
@@ -251,16 +251,27 @@ def test_io_run_conf_mrw(value):
     r2.read_run_conf(v)
     assert r2.mrw == r1.mrw
 
-def test_io_run_conf_specific_energy_spectrum_frequencies():
+def test_io_run_conf_specific_energy_spectrum_edges():
     r1 = RunConf()
-    r1.set_specific_energy_spectrum_frequencies(np.logspace(11., 16., 10))
+    r1.set_specific_energy_spectrum_bins(np.logspace(11., 16., 10))
     r1.set_n_photons(1, 2)
     v = virtual_file()
     r1.write_run_conf(v)
     r2 = RunConf()
     r2.read_run_conf(v)
-    np.testing.assert_allclose(r2.specific_energy_spectrum_frequencies,
-                               r1.specific_energy_spectrum_frequencies)
+    np.testing.assert_allclose(r2.specific_energy_spectrum_bin_edges,
+                               r1.specific_energy_spectrum_bin_edges)
+
+def test_io_run_conf_specific_energy_spectrum_invalid():
+    r1 = RunConf()
+    with pytest.raises(TypeError):
+        r1.set_specific_energy_spectrum_bins()
+    with pytest.raises(ValueError, match='at least two'):
+        r1.set_specific_energy_spectrum_bins([1.])
+    with pytest.raises(ValueError, match='strictly increasing'):
+        r1.set_specific_energy_spectrum_bins([1., 3., 2.])
+    with pytest.raises(ValueError, match='positive'):
+        r1.set_specific_energy_spectrum_bins([-1., 2.])
 
 @pytest.mark.parametrize(('value'), [False, True])
 def test_io_run_conf_convergence(value):
diff --git a/hyperion/grid/cartesian_grid.py b/hyperion/grid/cartesian_grid.py
index b1f00bde..d614bb82 100644
--- a/hyperion/grid/cartesian_grid.py
+++ b/hyperion/grid/cartesian_grid.py
@@ -279,7 +279,7 @@ def read_quantities(self, group, quantities='all'):
         # Read in physical quantities
         if quantities is not None:
             for quantity in group:
-                if quantity == 'specific_energy_spectrum_frequencies':
+                if quantity == 'specific_energy_spectrum_bin_edges':
                     continue  # per-frequency metadata, not a per-cell grid quantity
                 if quantities == 'all' or quantity in quantities:
                     array = np.array(group[quantity])
diff --git a/hyperion/grid/cylindrical_polar_grid.py b/hyperion/grid/cylindrical_polar_grid.py
index 2dc5b149..9b20af18 100644
--- a/hyperion/grid/cylindrical_polar_grid.py
+++ b/hyperion/grid/cylindrical_polar_grid.py
@@ -307,7 +307,7 @@ def read_quantities(self, group, quantities='all'):
         # Read in physical quantities
         if quantities is not None:
             for quantity in group:
-                if quantity == 'specific_energy_spectrum_frequencies':
+                if quantity == 'specific_energy_spectrum_bin_edges':
                     continue  # per-frequency metadata, not a per-cell grid quantity
                 if quantities == 'all' or quantity in quantities:
                     array = np.array(group[quantity])
diff --git a/hyperion/grid/octree_grid.py b/hyperion/grid/octree_grid.py
index 38485b16..565a1d2b 100644
--- a/hyperion/grid/octree_grid.py
+++ b/hyperion/grid/octree_grid.py
@@ -372,7 +372,7 @@ def read_quantities(self, group, quantities='all'):
         # Read in physical quantities
         if quantities is not None:
             for quantity in group:
-                if quantity == 'specific_energy_spectrum_frequencies':
+                if quantity == 'specific_energy_spectrum_bin_edges':
                     continue  # per-frequency metadata, not a per-cell grid quantity
                 if quantities == 'all' or quantity in quantities:
                     array = np.array(group[quantity])
diff --git a/hyperion/grid/spherical_polar_grid.py b/hyperion/grid/spherical_polar_grid.py
index e6d925d4..f4a71523 100644
--- a/hyperion/grid/spherical_polar_grid.py
+++ b/hyperion/grid/spherical_polar_grid.py
@@ -317,7 +317,7 @@ def read_quantities(self, group, quantities='all'):
         # Read in physical quantities
         if quantities is not None:
             for quantity in group:
-                if quantity == 'specific_energy_spectrum_frequencies':
+                if quantity == 'specific_energy_spectrum_bin_edges':
                     continue  # per-frequency metadata, not a per-cell grid quantity
                 if quantities == 'all' or quantity in quantities:
                     array = np.array(group[quantity])
diff --git a/hyperion/grid/voronoi_grid.py b/hyperion/grid/voronoi_grid.py
index cc12d4c7..24f42b1e 100644
--- a/hyperion/grid/voronoi_grid.py
+++ b/hyperion/grid/voronoi_grid.py
@@ -397,7 +397,7 @@ def read_quantities(self, group, quantities='all'):
         # Read in physical quantities
         if quantities is not None:
             for quantity in group:
-                if quantity == 'specific_energy_spectrum_frequencies':
+                if quantity == 'specific_energy_spectrum_bin_edges':
                     continue  # per-frequency metadata, not a per-cell grid quantity
                 if quantities == 'all' or quantity in quantities:
                     array = np.array(group[quantity])
diff --git a/hyperion/model/tests/test_specific_energy_spectrum.py b/hyperion/model/tests/test_specific_energy_spectrum.py
index 145fc545..132a22fe 100644
--- a/hyperion/model/tests/test_specific_energy_spectrum.py
+++ b/hyperion/model/tests/test_specific_energy_spectrum.py
@@ -25,6 +25,17 @@ def _read_dataset(filename, suffix):
         return np.asarray(f[matches[-1]][()], dtype=float)
 
 
+# Default bin edges for the tests, spanning the full frequency range of both
+# test dusts so that the spectrum captures all of the absorbed energy.
+_DEFAULT_EDGES = np.logspace(6., 18., 13)
+
+
+def _read_output_bin_edges(filename):
+    """Return the bin edges dataset written to the output (the last matching
+    dataset, which is the output one rather than the copied input)."""
+    return _read_dataset(filename, '/specific_energy_spectrum_bin_edges')
+
+
 def _assert_spectrum_sums_to_specific_energy(filename, exclude_empty_spectrum=False):
     # Summed over frequency, specific_energy_spectrum must reproduce specific_energy.
     # Cells with no absorption are clamped up to the minimum specific energy
@@ -46,7 +57,7 @@ def _assert_spectrum_sums_to_specific_energy(filename, exclude_empty_spectrum=Fa
 
 
 def _cartesian_model(output_specific_energy_spectrum='last', n_cells=2, n_photons=100000,
-                     density=1.e-18, frequencies=None):
+                     density=1.e-18, bin_edges=None):
     m = Model()
     edges = np.linspace(-1., 1., n_cells + 1)
     m.set_cartesian_grid(edges, edges, edges)
@@ -60,8 +71,8 @@ def _cartesian_model(output_specific_energy_spectrum='last', n_cells=2, n_photon
     m.set_seed(-12345)
     m.conf.output.output_specific_energy = 'last'
     m.conf.output.output_specific_energy_spectrum = output_specific_energy_spectrum
-    if frequencies is not None:
-        m.set_specific_energy_spectrum_frequencies(frequencies)
+    if output_specific_energy_spectrum != 'none':
+        m.set_specific_energy_spectrum_bins(_DEFAULT_EDGES if bin_edges is None else bin_edges)
     return m
 
 
@@ -104,16 +115,17 @@ def test_specific_energy_spectrum_sums_to_specific_energy_with_pda(tmpdir):
 
 
 @pytest.mark.requires_hyperion_binaries
-def test_specific_energy_spectrum_frequencies_written(tmpdir):
-    # specific_energy_spectrum_frequencies is a per-frequency (1-D) quantity, not
-    # per-cell. This guards against the regression where it was written as a grid
-    # array, which segfaulted whenever n_nu was smaller than n_cells.
-    m = _cartesian_model('last', n_cells=8)  # 512 cells, dust has 2 frequencies
+def test_specific_energy_spectrum_bin_edges_written(tmpdir):
+    # The bin edges are a per-frequency (1-D) quantity, not per-cell (this
+    # guards against the regression where the frequency metadata was written
+    # as a grid array, which segfaulted whenever n_nu was smaller than
+    # n_cells). The edges written to the output must match the ones given.
+    m = _cartesian_model('last', n_cells=8)  # 512 cells
     m.write(tmpdir.join(random_id()).strpath)
     out = m.run(tmpdir.join(random_id()).strpath)
-    bins = _read_dataset(out.filename, '/specific_energy_spectrum_frequencies')
-    assert bins is not None
-    assert bins.shape == (2,)  # == number of dust frequencies, not n_cells
+    edges = _read_output_bin_edges(out.filename)
+    assert edges is not None
+    np.testing.assert_allclose(edges, _DEFAULT_EDGES, rtol=1.e-6)
 
 
 @pytest.mark.requires_hyperion_binaries
@@ -159,27 +171,51 @@ def test_specific_energy_spectrum_amr(tmpdir):
     m.set_seed(-9)
     m.conf.output.output_specific_energy = 'last'
     m.conf.output.output_specific_energy_spectrum = 'last'
+    m.set_specific_energy_spectrum_bins(_DEFAULT_EDGES)
     m.write(tmpdir.join(random_id()).strpath)
     out = m.run(tmpdir.join(random_id()).strpath)
     _assert_spectrum_sums_to_specific_energy(out.filename)
 
 
 @pytest.mark.requires_hyperion_binaries
-def test_specific_energy_spectrum_custom_frequency_grid(tmpdir):
-    # A user-specified frequency grid (independent of the dust frequency grid)
-    # should be used for the output bins, and energy must still be conserved
-    # when summing over frequency.
-    frequencies = np.logspace(10., 16., 12)
-    m = _cartesian_model('last', density=1.e-16, frequencies=frequencies)
+def test_specific_energy_spectrum_custom_bin_edges(tmpdir):
+    # A user-specified frequency grid given as explicit bin edges: n + 1 edges
+    # should produce n bins, the edges should be written to the output, and
+    # energy must be conserved when summing over frequency as long as the edges
+    # span the frequency range over which energy is absorbed.
+    bin_edges = np.logspace(6., 18., 13)
+    m = _cartesian_model('last', density=1.e-16, bin_edges=bin_edges)
     m.write(tmpdir.join(random_id()).strpath)
     out = m.run(tmpdir.join(random_id()).strpath)
-    bins = _read_dataset(out.filename, '/specific_energy_spectrum_frequencies')
-    np.testing.assert_allclose(bins, frequencies, rtol=1.e-6)
+    edges_out = _read_output_bin_edges(out.filename)
+    np.testing.assert_allclose(edges_out, bin_edges, rtol=1.e-6)
     se_nu = _read_dataset(out.filename, '/specific_energy_spectrum')
-    assert se_nu.shape[0] == len(frequencies)
+    assert se_nu.shape[0] == len(bin_edges) - 1
     _assert_spectrum_sums_to_specific_energy(out.filename)
 
 
+@pytest.mark.requires_hyperion_binaries
+def test_specific_energy_spectrum_windowed_bin_edges(tmpdir):
+    # With explicit bin edges, energy absorbed from photons with frequencies
+    # outside the outer edges is not counted, so with edges covering only part
+    # of the frequency range over which energy is absorbed, the sum over the
+    # bins should recover only part of the scalar specific energy.
+    m_full = _cartesian_model('last', density=1.e-16, bin_edges=np.logspace(5., 20., 16))
+    m_window = _cartesian_model('last', density=1.e-16, bin_edges=np.logspace(14.8, 15.2, 5))
+    se_sum = {}
+    for key, m in (('full', m_full), ('window', m_window)):
+        m.write(tmpdir.join(random_id()).strpath)
+        out = m.run(tmpdir.join(random_id()).strpath)
+        se = _read_dataset(out.filename, '/specific_energy')
+        se_nu = _read_dataset(out.filename, '/specific_energy_spectrum')
+        se_sum[key] = se_nu.sum(axis=0)
+        # The spectrum can never contain more energy than the scalar total
+        assert np.all(se_sum[key] <= se * (1. + 1.e-6))
+    # The window contains part, but only part, of the absorbed energy (for a
+    # 6000 K source, roughly a quarter of the energy falls in this window)
+    assert 0. < se_sum['window'].sum() < 0.9 * se_sum['full'].sum()
+
+
 @pytest.mark.requires_hyperion_binaries
 def test_specific_energy_spectrum_with_mrw(tmpdir):
     # Energy deposited during modified-random-walk steps must appear in the
@@ -202,6 +238,7 @@ def test_specific_energy_spectrum_with_mrw(tmpdir):
     m.set_max_interactions(1000000000)
     m.conf.output.output_specific_energy = 'last'
     m.conf.output.output_specific_energy_spectrum = 'last'
+    m.set_specific_energy_spectrum_bins(_DEFAULT_EDGES)
     m.write(tmpdir.join(random_id()).strpath)
     out = m.run(tmpdir.join(random_id()).strpath)
     _assert_spectrum_sums_to_specific_energy(out.filename)
@@ -226,33 +263,49 @@ def test_specific_energy_spectrum_with_sublimation_cap(tmpdir):
     m.set_seed(-12345)
     m.conf.output.output_specific_energy = 'last'
     m.conf.output.output_specific_energy_spectrum = 'last'
+    m.set_specific_energy_spectrum_bins(_DEFAULT_EDGES)
     m.write(tmpdir.join(random_id()).strpath)
     out = m.run(tmpdir.join(random_id()).strpath)
     _assert_spectrum_sums_to_specific_energy(out.filename)
 
 
-def test_specific_energy_spectrum_frequencies_roundtrip(tmpdir):
-    # The custom frequency grid should survive a write/read round-trip.
+def test_specific_energy_spectrum_bins_roundtrip(tmpdir):
+    # The frequency bin edges should survive a write/read round-trip.
     from ...conf.conf_files import RunConf
-    frequencies = np.logspace(10., 16., 8)
+    edges = np.logspace(10., 16., 8)
     conf = RunConf()
-    conf.set_specific_energy_spectrum_frequencies(frequencies)
+    conf.set_specific_energy_spectrum_bins(edges)
     with h5py.File(tmpdir.join('conf.h5').strpath, 'w') as f:
-        conf._write_specific_energy_spectrum_frequencies(f)
+        conf._write_specific_energy_spectrum_bins(f)
         conf2 = RunConf()
-        conf2._read_specific_energy_spectrum_frequencies(f)
-    np.testing.assert_allclose(conf2.specific_energy_spectrum_frequencies, frequencies)
+        conf2._read_specific_energy_spectrum_bins(f)
+    np.testing.assert_allclose(conf2.specific_energy_spectrum_bin_edges, edges)
 
 
-def test_specific_energy_spectrum_frequencies_validation():
+def test_specific_energy_spectrum_bins_validation():
     from ...conf.conf_files import RunConf
     conf = RunConf()
-    with pytest.raises(ValueError):
-        conf.set_specific_energy_spectrum_frequencies([[1.e10, 1.e12], [1.e14, 1.e16]])
-    with pytest.raises(ValueError):
-        conf.set_specific_energy_spectrum_frequencies([1.e10, -1.e12])
-    # Default: no custom grid set
-    assert conf.specific_energy_spectrum_frequencies is None
+    with pytest.raises(TypeError):
+        conf.set_specific_energy_spectrum_bins()
+    with pytest.raises(ValueError, match='at least two'):
+        conf.set_specific_energy_spectrum_bins([1.e10])
+    with pytest.raises(ValueError, match='positive'):
+        conf.set_specific_energy_spectrum_bins([-1.e10, 1.e12])
+    with pytest.raises(ValueError, match='strictly increasing'):
+        conf.set_specific_energy_spectrum_bins([1.e10, 1.e16, 1.e12])
+    with pytest.raises(ValueError, match='1-d'):
+        conf.set_specific_energy_spectrum_bins([[1.e10, 1.e12], [1.e14, 1.e16]])
+    # Default: no bins set
+    assert conf.specific_energy_spectrum_bin_edges is None
+
+
+def test_specific_energy_spectrum_requires_bins(tmpdir):
+    # Enabling the frequency-resolved specific energy output without setting
+    # the bins must raise a clear error at write time.
+    m = _cartesian_model('none')
+    m.conf.output.output_specific_energy_spectrum = 'last'
+    with pytest.raises(ValueError, match='frequency bins have not been set'):
+        m.write(tmpdir.join(random_id()).strpath)
 
 
 @pytest.mark.requires_hyperion_binaries
@@ -275,6 +328,7 @@ def test_specific_energy_spectrum_voronoi(tmpdir):
     m.set_seed(-12345)
     m.conf.output.output_specific_energy = 'last'
     m.conf.output.output_specific_energy_spectrum = 'last'
+    m.set_specific_energy_spectrum_bins(_DEFAULT_EDGES)
     m.write(tmpdir.join(random_id()).strpath)
     out = m.run(tmpdir.join(random_id()).strpath)
     _assert_spectrum_sums_to_specific_energy(out.filename)
@@ -290,10 +344,10 @@ def test_specific_energy_spectrum_get_quantities(tmpdir):
     out = m.run(tmpdir.join(random_id()).strpath)
     g = out.get_quantities()
     assert 'specific_energy_spectrum' in g.quantities
-    assert 'specific_energy_spectrum_frequencies' not in g.quantities
+    assert 'specific_energy_spectrum_bin_edges' not in g.quantities
     se = np.array(g.quantities['specific_energy'])       # (dust, nz, ny, nx)
     se_nu = np.array(g.quantities['specific_energy_spectrum'])  # (nu, dust, nz, ny, nx)
-    assert se_nu.shape[0] == 2  # number of dust frequencies
+    assert se_nu.shape[0] == len(_DEFAULT_EDGES) - 1  # number of bins
     nu_sum = se_nu.sum(axis=0)
     heated = se > se.min() * (1. + 1.e-6)
     np.testing.assert_allclose(nu_sum[heated], se[heated], rtol=1.e-6)
diff --git a/hyperion/util/functions.py b/hyperion/util/functions.py
index 8712c8d9..c5b684d8 100644
--- a/hyperion/util/functions.py
+++ b/hyperion/util/functions.py
@@ -274,3 +274,4 @@ def asstr(s):
             return str(s)
     else:
         asstr = str
+
diff --git a/src/dust/dust_type_4elem.f90 b/src/dust/dust_type_4elem.f90
index ff90ce89..561c598b 100644
--- a/src/dust/dust_type_4elem.f90
+++ b/src/dust/dust_type_4elem.f90
@@ -753,7 +753,9 @@ function get_j_nu_bin_fractions(d, jnu_var_id, nu_edges) result(frac)
 
     ! Fraction of the emissivity falling into each of the bins whose edges
     ! are given by nu_edges (in ascending order). The fractions are
-    ! normalized so that they add up to exactly one.
+    ! normalized by the integral of the emissivity over all frequencies, so
+    ! they add up to the fraction of the emissivity that falls inside the
+    ! bins (which is one if the bins cover the full emissivity range).
 
     implicit none
 
@@ -770,7 +772,7 @@ function get_j_nu_bin_fractions(d, jnu_var_id, nu_edges) result(frac)
             &                      nu_edges(inu), nu_edges(inu+1))
     end do
 
-    norm = sum(frac)
+    norm = integral_loglog(d%j_nu(jnu_var_id)%x, d%j_nu(jnu_var_id)%pdf)
     if(norm > 0._dp) frac = frac / norm
 
   end function get_j_nu_bin_fractions
diff --git a/src/grid/grid_generic.f90 b/src/grid/grid_generic.f90
index c7ec2641..35bcdcd0 100644
--- a/src/grid/grid_generic.f90
+++ b/src/grid/grid_generic.f90
@@ -6,7 +6,7 @@ module grid_generic
   
   use grid_io, only : write_grid_3d, write_grid_4d, write_grid_5d
   use grid_geometry, only : geo
-  use grid_physics, only : n_photons, last_photon_id, specific_energy_sum, specific_energy_sum_spectrum, specific_energy, specific_energy_spectrum, nu_bins, density, density_original
+  use grid_physics, only : n_photons, last_photon_id, specific_energy_sum, specific_energy_sum_spectrum, specific_energy, specific_energy_spectrum, nu_bin_edges, density, density_original
   use settings, only : output_n_photons, output_specific_energy, output_specific_energy_spectrum, output_density, output_density_diff, physics_io_type, compute_specific_energy_spectrum
 
   implicit none
@@ -69,9 +69,9 @@ subroutine output_grid(group, iter, n_iter)
 
        if(trim(output_specific_energy_spectrum)=='all' .or. (trim(output_specific_energy_spectrum)=='last'.and.iter==n_iter)) then
 
-          ! The frequencies are a per-frequency quantity, not per-cell, so they
-          ! are written as a plain 1-D dataset rather than a grid array.
-          call mp_write_array(group, 'specific_energy_spectrum_frequencies', nu_bins)
+          ! The bin edges are a per-frequency quantity, not per-cell, so
+          ! they are written as a plain 1-D dataset rather than a grid array.
+          call mp_write_array(group, 'specific_energy_spectrum_bin_edges', nu_bin_edges)
 
           if(allocated(specific_energy_spectrum)) then
              select case(physics_io_type)
diff --git a/src/grid/grid_physics_3d.f90 b/src/grid/grid_physics_3d.f90
index 3a486186..851846ed 100644
--- a/src/grid/grid_physics_3d.f90
+++ b/src/grid/grid_physics_3d.f90
@@ -41,8 +41,11 @@ module grid_physics
   real(dp),allocatable, public :: specific_energy_spectrum(:,:,:) 
   real(dp),allocatable, public :: specific_energy_sum(:,:)
   real(dp),allocatable, public :: specific_energy_sum_spectrum(:,:,:)
-  real(dp),allocatable, public :: nu_bins(:)
-  real(dp),allocatable, public :: log_nu_bins(:)
+  ! Frequency bin edges for the specific energy spectrum: n_nu_bins+1
+  ! values that define the binning (photons outside the outer edges are not
+  ! counted), provided by the Python frontend.
+  real(dp),allocatable, public :: nu_bin_edges(:)
+  real(dp),allocatable, public :: log_nu_bin_edges(:)
 
   ! Fraction of the emissivity falling into each frequency bin, for each
   ! emissivity state and dust type - used to distribute energy deposited by
@@ -113,19 +116,15 @@ subroutine setup_grid_physics(group, use_mrw, use_pda, compute_specific_energy_s
     logical,intent(in) :: use_mrw, use_pda, compute_specific_energy_spectrum
     integer :: n_nu_bins
 
-    ! specific_energy_spectrum is binned onto a user-specified frequency grid if one was given,
-    ! otherwise onto the frequency grid of the first dust type. The spectrum
-    ! arrays are only allocated when the frequency-resolved specific energy is
-    ! actually requested (they can be large), so all whole-array operations on
-    ! them are guarded by allocated() or compute_specific_energy_spectrum
-    ! checks.
+    ! specific_energy_spectrum is binned onto the frequency bins defined by
+    ! the bin edges provided by the Python frontend. The spectrum arrays are
+    ! only allocated when the frequency-resolved specific energy is actually
+    ! requested (they can be large), so all whole-array operations on them
+    ! are guarded by allocated() or compute_specific_energy_spectrum checks.
     if (compute_specific_energy_spectrum) then
-       if (allocated(specific_energy_spectrum_frequencies)) then
-          n_nu_bins = size(specific_energy_spectrum_frequencies)
-       else if (n_dust > 0) then
-          n_nu_bins = d(1)%n_nu
-       else
-          n_nu_bins = 0
+       n_nu_bins = size(specific_energy_spectrum_bin_edges) - 1
+       if (any(specific_energy_spectrum_bin_edges(2:) <= specific_energy_spectrum_bin_edges(:n_nu_bins))) then
+          call error("setup_grid_physics", "specific_energy_spectrum_bin_edges should be strictly increasing")
        end if
     else
        n_nu_bins = 0
@@ -272,20 +271,15 @@ subroutine setup_grid_physics(group, use_mrw, use_pda, compute_specific_energy_s
        specific_energy_sum_spectrum = 0._dp
     end if
 
-    ! Cache the frequency grid (and its log) once so they do not have to be
-    ! rebuilt for every photon. Photons are binned to the nearest grid point in
-    ! log-frequency space (see grid_propagate).
+    ! Cache the frequency bins (and the log of the edges) once so they do
+    ! not have to be rebuilt for every photon. Photons are binned by
+    ! locating their frequency in the bin edges in log space (see
+    ! grid_propagate).
     if (compute_specific_energy_spectrum) then
-       allocate(nu_bins(n_nu_bins))
-       if (allocated(specific_energy_spectrum_frequencies)) then
-          nu_bins = specific_energy_spectrum_frequencies
-       else
-          do idx=1,n_nu_bins
-             nu_bins(idx) = d(1)%nu(idx)
-          end do
-       end if
-       allocate(log_nu_bins(n_nu_bins))
-       log_nu_bins = log10(nu_bins)
+       allocate(nu_bin_edges(n_nu_bins+1))
+       nu_bin_edges = specific_energy_spectrum_bin_edges
+       allocate(log_nu_bin_edges(n_nu_bins+1))
+       log_nu_bin_edges = log10(nu_bin_edges)
 
        call setup_j_nu_bin_fractions()
     end if
@@ -333,38 +327,23 @@ subroutine setup_j_nu_bin_fractions()
     ! Pre-compute, for each dust type and emissivity state, the fraction of
     ! the emissivity falling into each frequency bin, stored in
     ! j_nu_bin_frac. This is used to distribute the energy deposited by the
-    ! MRW over the specific energy spectrum without sampling. The bins are
-    ! defined by the same nearest-neighbour convention in log-frequency as
-    ! the photon binning in grid_propagate, so the edges are placed half-way
-    ! (in log space) between the bin frequencies, with the outer bins
-    ! extending to all lower/higher frequencies.
+    ! MRW over the specific energy spectrum without sampling.
 
     implicit none
 
-    integer :: inu, id, iv, n_nu_bins, n_jnu_max
-    real(dp),allocatable :: nu_bin_edges(:)
-
-    n_nu_bins = size(log_nu_bins)
-
-    allocate(nu_bin_edges(n_nu_bins+1))
-    nu_bin_edges(1) = 10._dp**(log_nu_bins(1) - 99._dp)
-    do inu=2,n_nu_bins
-       nu_bin_edges(inu) = 10._dp**(0.5_dp * (log_nu_bins(inu-1) + log_nu_bins(inu)))
-    end do
-    nu_bin_edges(n_nu_bins+1) = 10._dp**(log_nu_bins(n_nu_bins) + 99._dp)
+    integer :: id, iv, n_jnu_max
 
     n_jnu_max = 0
     do id=1,n_dust
        n_jnu_max = max(n_jnu_max, d(id)%n_jnu)
     end do
-    allocate(j_nu_bin_frac(n_nu_bins, n_jnu_max, n_dust))
+    allocate(j_nu_bin_frac(size(nu_bin_edges)-1, n_jnu_max, n_dust))
     j_nu_bin_frac = 0._dp
     do id=1,n_dust
        do iv=1,d(id)%n_jnu
           j_nu_bin_frac(:, iv, id) = get_j_nu_bin_fractions(d(id), iv, nu_bin_edges)
        end do
     end do
-    deallocate(nu_bin_edges)
 
   end subroutine setup_j_nu_bin_fractions
 
diff --git a/src/grid/grid_propagate_3d.f90 b/src/grid/grid_propagate_3d.f90
index 29b9f22f..162af0d5 100644
--- a/src/grid/grid_propagate_3d.f90
+++ b/src/grid/grid_propagate_3d.f90
@@ -5,7 +5,7 @@ module grid_propagate
   use type_grid_cell
   use dust_main, only : n_dust
   use grid_geometry, only : escaped, find_wall, in_correct_cell, next_cell, opposite_wall
-  use grid_physics, only : specific_energy_sum, specific_energy_sum_spectrum, log_nu_bins, density, n_photons, last_photon_id
+  use grid_physics, only : specific_energy_sum, specific_energy_sum_spectrum, log_nu_bin_edges, density, n_photons, last_photon_id
   use sources
   use counters
   use settings, only : frac_check => propagation_check_frequency, compute_specific_energy_spectrum
@@ -63,7 +63,9 @@ subroutine grid_integrate(p,tau_required,tau_achieved)
     ! (re-emission and scattering happen between calls, which is also when the
     ! cached opacities in p%current_kappa are refreshed), so the frequency bin
     ! only needs to be found once per call rather than at every cell crossing.
-    if (compute_specific_energy_spectrum) idx = minloc(abs(log_nu_bins - log10(p%nu)), DIM=1)
+    ! locate returns -1 for photons outside the outer bin edges, which are
+    ! not counted in the spectrum.
+    if (compute_specific_energy_spectrum) idx = locate(log_nu_bin_edges, log10(p%nu))
 
     radial = (p%r .dot. p%v) > 0.
 
@@ -147,7 +149,7 @@ subroutine grid_integrate(p,tau_required,tau_achieved)
              if(density(p%icell%ic, id) > 0._dp) then
                 specific_energy_sum(p%icell%ic, id) = &
                      & specific_energy_sum(p%icell%ic, id) + tmin * p%current_kappa(id) * p%energy
-                if (compute_specific_energy_spectrum) then
+                if (compute_specific_energy_spectrum .and. idx > 0) then
                    specific_energy_sum_spectrum(p%icell%ic, id, idx) = &
                         & specific_energy_sum_spectrum(p%icell%ic, id, idx) + tmin * p%current_kappa(id) * p%energy
                 end if
@@ -209,7 +211,7 @@ subroutine grid_integrate(p,tau_required,tau_achieved)
           ! energy but missing from the binned spectrum, which biases
           ! the spectrum low in any cell optically thick enough for
           ! photons to interact within a single crossing.
-          if (compute_specific_energy_spectrum) then
+          if (compute_specific_energy_spectrum .and. idx > 0) then
              do id=1,n_dust
                 if(density(p%icell%ic, id) > 0._dp) then
                    specific_energy_sum_spectrum(p%icell%ic, id, idx) = &
diff --git a/src/main/settings.f90 b/src/main/settings.f90
index 5fd3bf65..ee2fcf51 100644
--- a/src/main/settings.f90
+++ b/src/main/settings.f90
@@ -30,9 +30,10 @@ module settings
   real(dp) :: monochromatic_energy_threshold
   real(dp),allocatable :: frequencies(:)
 
-  ! Optional user-specified frequency grid for specific_energy_spectrum. If not
-  ! allocated, the frequency grid of the first dust type is used instead.
-  real(dp),allocatable :: specific_energy_spectrum_frequencies(:)
+  ! Frequency bin edges for specific_energy_spectrum, provided by the Python
+  ! frontend. The edges define the binning: photons outside the outer edges
+  ! are not counted.
+  real(dp),allocatable :: specific_energy_spectrum_bin_edges(:)
 
   integer :: physics_io_type
 
diff --git a/src/main/setup_rt.f90 b/src/main/setup_rt.f90
index 9bee928c..76225fdb 100644
--- a/src/main/setup_rt.f90
+++ b/src/main/setup_rt.f90
@@ -93,11 +93,14 @@ subroutine setup_initial(input_handle)
 
     compute_specific_energy_spectrum = trim(output_specific_energy_spectrum).ne.'none'
 
-    ! Optional user-specified frequency grid (otherwise the dust grid is used)
+    ! The frequency bin edges are written by the Python frontend when the
+    ! spectrum is requested
     if(compute_specific_energy_spectrum) then
-       if(mp_path_exists(input_handle, 'specific_energy_spectrum_frequencies')) then
-          call mp_table_read_column_auto(input_handle, 'specific_energy_spectrum_frequencies', 'nu', specific_energy_spectrum_frequencies)
+       if(.not.mp_path_exists(input_handle, 'specific_energy_spectrum_bin_edges')) then
+          call error("setup_initial","specific_energy_spectrum_bin_edges should be present in the input "// &
+               & "when output_specific_energy_spectrum is enabled")
        end if
+       call mp_table_read_column_auto(input_handle, 'specific_energy_spectrum_bin_edges', 'nu', specific_energy_spectrum_bin_edges)
     end if
 
     if(use_mrw) then

From ae50b0ae7abf70597968860d371b0052e207b0b1 Mon Sep 17 00:00:00 2001
From: Thomas Robitaille <thomas.robitaille@gmail.com>
Date: Sat, 1 Aug 2026 22:45:28 +0000
Subject: [PATCH 3/3] Initialize the frequency bin index to an out-of-range
 value so the spectrum deposit guards never read an undefined value when the
 spectrum is not being computed

---
 src/grid/grid_propagate_3d.f90 | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/src/grid/grid_propagate_3d.f90 b/src/grid/grid_propagate_3d.f90
index 162af0d5..ae9453e7 100644
--- a/src/grid/grid_propagate_3d.f90
+++ b/src/grid/grid_propagate_3d.f90
@@ -64,7 +64,10 @@ subroutine grid_integrate(p,tau_required,tau_achieved)
     ! cached opacities in p%current_kappa are refreshed), so the frequency bin
     ! only needs to be found once per call rather than at every cell crossing.
     ! locate returns -1 for photons outside the outer bin edges, which are
-    ! not counted in the spectrum.
+    ! not counted in the spectrum. idx is set to -1 up front so the deposit
+    ! guards below (idx > 0) are always reading a defined value, including
+    ! when the frequency-resolved spectrum is not being computed.
+    idx = -1
     if (compute_specific_energy_spectrum) idx = locate(log_nu_bin_edges, log10(p%nu))
 
     radial = (p%r .dot. p%v) > 0.