Use bin edges internally for frequency-dependent specific energy

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::

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)

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):

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])

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])

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])

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])

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])

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'