From 9f049057c00c67d2833545ccea26c25e837449e7 Mon Sep 17 00:00:00 2001
From: Martin Glatzle <mglatzle@mpa-garching.mpg.de>
Date: Fri, 24 Jul 2020 17:27:46 +0200
Subject: [PATCH 01/26] Interpolators depend on scipy.

---
 setup.cfg | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/setup.cfg b/setup.cfg
index e8acee07..d1fdcb19 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -17,7 +17,7 @@ github_project = heloises/hoki
 # install_requires should be formatted as a semicolon-separated list, e.g.:
 zip_safe = False
 packages =  find:
-install_requires = 
+install_requires =
 	astropy
 	numpy
 	pandas
@@ -26,6 +26,7 @@ install_requires =
 	wheel
 	pysynphot
 	numba
+	scipy
 setup_requires = setuptools_scm
 python_requires = >=3.5
 #use_2to3 = False

From aa8a8250a1695f0625c94766c1e10334229118b8 Mon Sep 17 00:00:00 2001
From: Martin Glatzle <mglatzle@mpa-garching.mpg.de>
Date: Fri, 24 Jul 2020 17:28:16 +0200
Subject: [PATCH 02/26] Implement base class for interpolation.

---
 hoki/interpolators.py            | 179 +++++++++++++++++++++++++++++++
 hoki/tests/test_interpolators.py |  26 +++++
 2 files changed, 205 insertions(+)
 create mode 100644 hoki/interpolators.py
 create mode 100644 hoki/tests/test_interpolators.py

diff --git a/hoki/interpolators.py b/hoki/interpolators.py
new file mode 100644
index 00000000..2d05ae05
--- /dev/null
+++ b/hoki/interpolators.py
@@ -0,0 +1,179 @@
+"""
+Author: Martin Glatzle
+
+Tools that allow to interpolate BPASS quantities on metallicity-age grids.
+"""
+
+import numpy as np
+from scipy import interpolate
+from hoki.constants import BPASS_NUM_METALLICITIES, BPASS_TIME_BINS
+import warnings
+
+
+class GridInterpolator():
+    """
+    Interpolate BPASS quantities on a metallicity-age grid.
+
+    Base class for the interpolation of BPASS data (possibly vector valued,
+    e.g. spectra) given on a grid of metallicities and ages. Using this class
+    directly is discouraged, it has no public methods or attributes. It is
+    advisable to implement convenience child classes for each specific
+    quantity.
+
+    Parameters
+    ----------
+    grid : `numpy.ndarray` (N_Z, N_a, D)
+        A 3D numpy array containing the `D`-dimensional quantity to be
+        interpolated as a function of metallicity and age.
+    metallicities : `numpy.ndarray` (N_Z,), optional
+        The metallicities at which `grid` is evaluated. Defaults to the full
+        array of BPASS output metallicities.
+    ages : `numpy.ndarray` (N_a,), optional
+        The ages at which `grid` is evaluated. Defaults to the full array of
+        BPASS output ages in log scale.
+    dtype : `type`, optional
+        The data type to be used by an instance of this class. Defaults to
+        `numpy.float64`. Can be used to reduce memory footprint.
+
+    Notes
+    -----
+    Uses `scipy.interpolate.LinearNDInterpolator`, which performs triangulation
+    of the input data and linear barycentric interpolation on each
+    triangle. Support for other interpolation methods should be easy to
+    implemented.
+    """
+    def __init__(self, grid,
+                 metallicities=BPASS_NUM_METALLICITIES,
+                 ages=BPASS_TIME_BINS,
+                 dtype=np.float64):
+        for arr, ndim in zip([grid, metallicities, ages], [3, 1, 1]):
+            if not np.ndim(arr) == ndim:
+                raise ValueError("Wrong dimensionality of input arrays.")
+        if not grid.shape[0] == len(metallicities):
+            raise ValueError(
+                "Shapes of `grid` and `metallicities` are incompatible."
+            )
+        if not grid.shape[1] == len(ages):
+            raise ValueError(
+                "Shapes of `grid` and `ages` are incompatible."
+            )
+        self._dtype = dtype
+        self._metallicities = metallicities.astype(self._dtype, copy=True)
+        self._zMin = np.amin(self._metallicities)
+        self._zMax = np.amax(self._metallicities)
+        self._ages = ages.astype(self._dtype, copy=True)
+        self._aMin = np.amin(self._ages)
+        self._aMax = np.amax(self._ages)
+        self._construct_interpolator(
+            grid.astype(dtype, copy=False),
+        )
+        return
+
+    def _construct_interpolator(self, grid):
+        """
+        Construct an interpolator on metallicity-age grid.
+
+        Parameters
+        ----------
+        grid : `numpy.ndarray` (N_Z, N_a, D)
+            A 3D numpy array containing the `D`-dimensional quantity to be
+            interpolated as a function of metallicity and age.
+        """
+        zz, aa = np.meshgrid(self._metallicities, self._ages, indexing='ij')
+        points = np.stack((zz, aa), -1).reshape((-1, 2))
+        self._interpolator = interpolate.LinearNDInterpolator(
+            points, grid.reshape((-1, grid.shape[2]))
+        )
+        return
+
+    def _interpolate(self, metallicities, ages, masses=1):
+        """
+        Perform interpolation on this instance's grid.
+
+        Parameters
+        ----------
+        metallicities : `numpy.ndarray` (N,)
+            Stellar metallicities at which `grid` to interpolate.
+        ages : `numpy.ndarray` (N,)
+            Stellar ages at which to interpolate.
+        masses : `numpy.ndarray` (N,) or `float`, optional
+            Stellar masses in units of 1e6 M_\\odot. Used to scale the
+            interpolation result. Defaults to unity.
+
+        Returns
+        -------
+         : `numpy.ndarray` (N, D)
+        """
+        # check dtypes
+        if metallicities.dtype != self._dtype:
+            warnings.warn(
+                "Input metallicities for interpolation of wrong dtype, "
+                "attempting copy and cast.",
+                UserWarning
+            )
+            metallicities = np.array(
+                metallicities, dtype=self._dtype
+            )
+        if ages.dtype != self._dtype:
+            warnings.warn(
+                "Input ages for interpolation of wrong dtype, "
+                "attempting copy and cast.",
+                UserWarning
+            )
+            ages = np.array(
+                ages, dtype=self._dtype
+            )
+
+        # clipping
+        if np.amax(metallicities) > self._zMax or \
+           np.amin(metallicities) < self._zMin:
+            warnings.warn(
+                "Input metallicities for interpolation outside of "
+                "available range "+str(self._zMin)+", "+str(self._zMax) +
+                " provided. They will be clipped.",
+                UserWarning
+            )
+            metallicities = np.clip(metallicities, self._zMin, self._zMax)
+        if np.amax(ages) > 10**(self._aMax) or \
+           np.amin(ages) < 10**(self._aMin):
+            warnings.warn(
+                "Input ages for interpolation outside of available "
+                "range "+str(10**self._aMin)+", "+str(10**self._aMax) +
+                " [yr] provided. They will be clipped.",
+                UserWarning
+            )
+            ages = np.clip(ages, self._aMin, self._aMax)
+        return self._interpolator(metallicities, ages) * \
+            self._check_masses(masses)
+
+    def _check_masses(self, masses):
+        """
+        Make sure `masses` has correct dimensionality.
+
+        Reshapes `masses` if it is a 1D array, does nothing if it is
+        scalar. Also warns about masses potentially being too small.
+
+        Parameters
+        ----------
+        masses : `numpy.ndarray` (N,) or `float`
+            Stellar masses in units of 1e6 M_\\odot.
+
+        Returns
+        -------
+        masses : `numpy.ndarray` (N, 1) or `float`
+            Input reshaped to be multiplicable by interpolation result.
+        """
+        if np.amin(masses) < 1e-3:
+            warnings.warn(
+                "Input masses below 1000 M_sun! For such small populations,"
+                " single stars can contribute a significant fraction of the"
+                " population mass and re-scaling BPASS values averaged over"
+                " more massive populations likely yields incorrect results.",
+                UserWarning
+            )
+        if np.ndim(masses) == 0:
+            return masses
+        elif np.ndim(masses) == 1:
+            return masses[:, None]
+        else:
+            raise ValueError("Wrong dimensionality of `masses`.")
diff --git a/hoki/tests/test_interpolators.py b/hoki/tests/test_interpolators.py
new file mode 100644
index 00000000..0d5e22d0
--- /dev/null
+++ b/hoki/tests/test_interpolators.py
@@ -0,0 +1,26 @@
+from unittest import TestCase
+from hoki import interpolators
+from hoki.constants import BPASS_NUM_METALLICITIES, BPASS_TIME_BINS
+import numpy as np
+
+
+class TestGridInterpolator(TestCase):
+
+    def test_init(self):
+        for i in [1, 2, 3, 4]:
+            grid = np.ones(
+                (len(BPASS_NUM_METALLICITIES), len(BPASS_TIME_BINS), i)
+            )
+            interpolators.GridInterpolator(grid)
+
+        with self.assertRaises(ValueError):
+            interpolators.GridInterpolator(grid, metallicities=1)
+        with self.assertRaises(ValueError):
+            interpolators.GridInterpolator(grid, ages=1)
+        with self.assertRaises(ValueError):
+            interpolators.GridInterpolator(grid,
+                                           metallicities=np.array([1, 2]))
+        with self.assertRaises(ValueError):
+            interpolators.GridInterpolator(grid,
+                                           ages=np.array([1, 2]))
+        return

From cae89c079cdf755ccfd839b74a6b24d2c2ce942d Mon Sep 17 00:00:00 2001
From: Martin Glatzle <mglatzle@mpa-garching.mpg.de>
Date: Wed, 29 Jul 2020 15:22:40 +0200
Subject: [PATCH 03/26] Docstring update.

---
 hoki/interpolators.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/hoki/interpolators.py b/hoki/interpolators.py
index 2d05ae05..d745e5d4 100644
--- a/hoki/interpolators.py
+++ b/hoki/interpolators.py
@@ -103,6 +103,7 @@ def _interpolate(self, metallicities, ages, masses=1):
         Returns
         -------
          : `numpy.ndarray` (N, D)
+            Interpolation result.
         """
         # check dtypes
         if metallicities.dtype != self._dtype:

From 28e932ccd3ee6641330611847354608306af5109 Mon Sep 17 00:00:00 2001
From: Martin Glatzle <mglatzle@mpa-garching.mpg.de>
Date: Wed, 19 Aug 2020 17:04:36 +0200
Subject: [PATCH 04/26] Use != where appropriate.

---
 hoki/interpolators.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/hoki/interpolators.py b/hoki/interpolators.py
index d745e5d4..d9e96dfe 100644
--- a/hoki/interpolators.py
+++ b/hoki/interpolators.py
@@ -47,13 +47,13 @@ def __init__(self, grid,
                  ages=BPASS_TIME_BINS,
                  dtype=np.float64):
         for arr, ndim in zip([grid, metallicities, ages], [3, 1, 1]):
-            if not np.ndim(arr) == ndim:
+            if np.ndim(arr) != ndim:
                 raise ValueError("Wrong dimensionality of input arrays.")
-        if not grid.shape[0] == len(metallicities):
+        if grid.shape[0] != len(metallicities):
             raise ValueError(
                 "Shapes of `grid` and `metallicities` are incompatible."
             )
-        if not grid.shape[1] == len(ages):
+        if grid.shape[1] != len(ages):
             raise ValueError(
                 "Shapes of `grid` and `ages` are incompatible."
             )

From b0d50d7204c58e470599c56b835f20b6ebcf2b27 Mon Sep 17 00:00:00 2001
From: Martin Glatzle <mglatzle@mpa-garching.mpg.de>
Date: Fri, 21 Aug 2020 11:07:16 +0200
Subject: [PATCH 05/26] Update docstrings.

---
 hoki/interpolators.py | 14 ++++++++------
 1 file changed, 8 insertions(+), 6 deletions(-)

diff --git a/hoki/interpolators.py b/hoki/interpolators.py
index d9e96dfe..b410d47e 100644
--- a/hoki/interpolators.py
+++ b/hoki/interpolators.py
@@ -1,7 +1,7 @@
 """
-Author: Martin Glatzle
-
 Tools that allow to interpolate BPASS quantities on metallicity-age grids.
+
+Author: Martin Glatzle
 """
 
 import numpy as np
@@ -39,8 +39,8 @@ class GridInterpolator():
     -----
     Uses `scipy.interpolate.LinearNDInterpolator`, which performs triangulation
     of the input data and linear barycentric interpolation on each
-    triangle. Support for other interpolation methods should be easy to
-    implemented.
+    triangle. Support for other interpolation methods should be fairly easy to
+    implement.
     """
     def __init__(self, grid,
                  metallicities=BPASS_NUM_METALLICITIES,
@@ -93,9 +93,11 @@ def _interpolate(self, metallicities, ages, masses=1):
         Parameters
         ----------
         metallicities : `numpy.ndarray` (N,)
-            Stellar metallicities at which `grid` to interpolate.
+            Stellar metallicities at which to interpolate. Same units as those
+            used in the construction of this instance.
         ages : `numpy.ndarray` (N,)
-            Stellar ages at which to interpolate.
+            Stellar ages at which to interpolate. Same units as those used in
+            the construction of this instance.
         masses : `numpy.ndarray` (N,) or `float`, optional
             Stellar masses in units of 1e6 M_\\odot. Used to scale the
             interpolation result. Defaults to unity.

From 771ef1a2d9b3778b42da88d7a2c3d0f93b6d169d Mon Sep 17 00:00:00 2001
From: Martin Glatzle <mglatzle@mpa-garching.mpg.de>
Date: Fri, 21 Aug 2020 11:08:22 +0200
Subject: [PATCH 06/26] Use f-strings and fix age units.

---
 hoki/interpolators.py | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/hoki/interpolators.py b/hoki/interpolators.py
index b410d47e..e030179d 100644
--- a/hoki/interpolators.py
+++ b/hoki/interpolators.py
@@ -131,18 +131,18 @@ def _interpolate(self, metallicities, ages, masses=1):
         if np.amax(metallicities) > self._zMax or \
            np.amin(metallicities) < self._zMin:
             warnings.warn(
-                "Input metallicities for interpolation outside of "
-                "available range "+str(self._zMin)+", "+str(self._zMax) +
-                " provided. They will be clipped.",
+                "Input metallicities for interpolation outside of available "
+                f"range {self._zMin} -- {self._zMax} provided. "
+                "They will be clipped.",
                 UserWarning
             )
             metallicities = np.clip(metallicities, self._zMin, self._zMax)
-        if np.amax(ages) > 10**(self._aMax) or \
-           np.amin(ages) < 10**(self._aMin):
+        if np.amax(ages) > self._aMax or \
+           np.amin(ages) < self._aMin:
             warnings.warn(
                 "Input ages for interpolation outside of available "
-                "range "+str(10**self._aMin)+", "+str(10**self._aMax) +
-                " [yr] provided. They will be clipped.",
+                f"range {self._aMin} -- {self._aMax} provided. "
+                "They will be clipped.",
                 UserWarning
             )
             ages = np.clip(ages, self._aMin, self._aMax)

From b1d7a12f1f5aaf35c6102298ef0e580981cde397 Mon Sep 17 00:00:00 2001
From: Martin Glatzle <mglatzle@mpa-garching.mpg.de>
Date: Fri, 21 Aug 2020 11:08:40 +0200
Subject: [PATCH 07/26] Make method static.

---
 hoki/interpolators.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/hoki/interpolators.py b/hoki/interpolators.py
index e030179d..07776a1a 100644
--- a/hoki/interpolators.py
+++ b/hoki/interpolators.py
@@ -149,7 +149,8 @@ def _interpolate(self, metallicities, ages, masses=1):
         return self._interpolator(metallicities, ages) * \
             self._check_masses(masses)
 
-    def _check_masses(self, masses):
+    @staticmethod
+    def _check_masses(masses):
         """
         Make sure `masses` has correct dimensionality.
 

From 3565e37278cb915fbdcfdf25120b600536941656 Mon Sep 17 00:00:00 2001
From: Martin Glatzle <mglatzle@mpa-garching.mpg.de>
Date: Fri, 21 Aug 2020 11:41:33 +0200
Subject: [PATCH 08/26] Split interpolation of mass scaled qtties off into
 child class.

---
 hoki/interpolators.py | 64 +++++++++++++++++++++++++++++++++++++------
 1 file changed, 56 insertions(+), 8 deletions(-)

diff --git a/hoki/interpolators.py b/hoki/interpolators.py
index 07776a1a..4535f89f 100644
--- a/hoki/interpolators.py
+++ b/hoki/interpolators.py
@@ -14,10 +14,10 @@ class GridInterpolator():
     """
     Interpolate BPASS quantities on a metallicity-age grid.
 
-    Base class for the interpolation of BPASS data (possibly vector valued,
-    e.g. spectra) given on a grid of metallicities and ages. Using this class
+    Base class for the interpolation of BPASS quantities (possibly vector
+    valued) given on a grid of metallicities and ages. Using this class
     directly is discouraged, it has no public methods or attributes. It is
-    advisable to implement convenience child classes for each specific
+    advisable to use/implement convenience child classes for each specific
     quantity.
 
     Parameters
@@ -86,7 +86,7 @@ def _construct_interpolator(self, grid):
         )
         return
 
-    def _interpolate(self, metallicities, ages, masses=1):
+    def _interpolate(self, metallicities, ages):
         """
         Perform interpolation on this instance's grid.
 
@@ -98,9 +98,6 @@ def _interpolate(self, metallicities, ages, masses=1):
         ages : `numpy.ndarray` (N,)
             Stellar ages at which to interpolate. Same units as those used in
             the construction of this instance.
-        masses : `numpy.ndarray` (N,) or `float`, optional
-            Stellar masses in units of 1e6 M_\\odot. Used to scale the
-            interpolation result. Defaults to unity.
 
         Returns
         -------
@@ -146,7 +143,58 @@ def _interpolate(self, metallicities, ages, masses=1):
                 UserWarning
             )
             ages = np.clip(ages, self._aMin, self._aMax)
-        return self._interpolator(metallicities, ages) * \
+        return self._interpolator(metallicities, ages)
+
+
+class GridInterpolatorMassScaled(GridInterpolator):
+    """
+    Interpolate BPASS quantities that scale with population mass on a
+    metallicity-age grid.
+
+    Base class for the interpolation of BPASS quantities (possibly vector
+    valued) which scale with stellar population mass given on a grid of
+    metallicities and ages. Using this class directly is discouraged, it has no
+    public methods or attributes. It is advisable to use/implement convenience
+    child classes for each specific quantity.
+
+    Parameters
+    ----------
+    grid : `numpy.ndarray` (N_Z, N_a, D)
+        A 3D numpy array containing the `D`-dimensional quantity to be
+        interpolated as a function of metallicity and age, normalized to a
+        population mass of 1e6 M_\\odot.
+    metallicities : `numpy.ndarray` (N_Z,), optional
+        The metallicities at which `grid` is evaluated. Defaults to the full
+        array of BPASS output metallicities.
+    ages : `numpy.ndarray` (N_a,), optional
+        The ages at which `grid` is evaluated. Defaults to the full array of
+        BPASS output ages in log scale.
+    dtype : `type`, optional
+        The data type to be used by an instance of this class. Defaults to
+        `numpy.float64`. Can be used to reduce memory footprint.
+    """
+    def _interpolate(self, metallicities, ages, masses=1):
+        """
+        Perform interpolation on this instance's grid.
+
+        Parameters
+        ----------
+        metallicities : `numpy.ndarray` (N,)
+            Stellar metallicities at which to interpolate. Same units as those
+            used in the construction of this instance.
+        ages : `numpy.ndarray` (N,)
+            Stellar ages at which to interpolate. Same units as those used in
+            the construction of this instance.
+        masses : `numpy.ndarray` (N,) or `float`, optional
+            Stellar population masses in units of 1e6 M_\\odot. Used to scale
+            the interpolation result. Defaults to unity.
+
+        Returns
+        -------
+         : `numpy.ndarray` (N, D)
+            Interpolation result.
+        """
+        return super()._interpolate(metallicities, ages) * \
             self._check_masses(masses)
 
     @staticmethod

From ea84e1f293c0ee491322fc5ec01f84ba4ca973d2 Mon Sep 17 00:00:00 2001
From: Martin Glatzle <mglatzle@mpa-garching.mpg.de>
Date: Fri, 21 Aug 2020 12:58:45 +0200
Subject: [PATCH 09/26] Enable scalar input.

---
 hoki/interpolators.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/hoki/interpolators.py b/hoki/interpolators.py
index 4535f89f..11881447 100644
--- a/hoki/interpolators.py
+++ b/hoki/interpolators.py
@@ -105,7 +105,7 @@ def _interpolate(self, metallicities, ages):
             Interpolation result.
         """
         # check dtypes
-        if metallicities.dtype != self._dtype:
+        if np.asarray(metallicities).dtype != self._dtype:
             warnings.warn(
                 "Input metallicities for interpolation of wrong dtype, "
                 "attempting copy and cast.",
@@ -114,7 +114,7 @@ def _interpolate(self, metallicities, ages):
             metallicities = np.array(
                 metallicities, dtype=self._dtype
             )
-        if ages.dtype != self._dtype:
+        if np.asarray(ages).dtype != self._dtype:
             warnings.warn(
                 "Input ages for interpolation of wrong dtype, "
                 "attempting copy and cast.",

From e35a89caec807389ee8a146ae974b2ee0f6888ce Mon Sep 17 00:00:00 2001
From: Martin Glatzle <mglatzle@mpa-garching.mpg.de>
Date: Fri, 21 Aug 2020 12:58:58 +0200
Subject: [PATCH 10/26] Expand test suite.

---
 hoki/tests/test_interpolators.py | 135 +++++++++++++++++++++++++++++--
 1 file changed, 128 insertions(+), 7 deletions(-)

diff --git a/hoki/tests/test_interpolators.py b/hoki/tests/test_interpolators.py
index 0d5e22d0..8248f806 100644
--- a/hoki/tests/test_interpolators.py
+++ b/hoki/tests/test_interpolators.py
@@ -6,21 +6,142 @@
 
 class TestGridInterpolator(TestCase):
 
+    @classmethod
+    def setUpClass(cls):
+        cls.tested_class = interpolators.GridInterpolator
+        return
+
     def test_init(self):
         for i in [1, 2, 3, 4]:
             grid = np.ones(
                 (len(BPASS_NUM_METALLICITIES), len(BPASS_TIME_BINS), i)
             )
-            interpolators.GridInterpolator(grid)
+            self.__class__.tested_class(grid)
 
         with self.assertRaises(ValueError):
-            interpolators.GridInterpolator(grid, metallicities=1)
+            self.__class__.tested_class(grid, metallicities=1)
         with self.assertRaises(ValueError):
-            interpolators.GridInterpolator(grid, ages=1)
+            self.__class__.tested_class(grid, ages=1)
         with self.assertRaises(ValueError):
-            interpolators.GridInterpolator(grid,
-                                           metallicities=np.array([1, 2]))
+            self.__class__.tested_class(grid,
+                                        metallicities=np.array([1, 2]))
         with self.assertRaises(ValueError):
-            interpolators.GridInterpolator(grid,
-                                           ages=np.array([1, 2]))
+            self.__class__.tested_class(grid,
+                                        ages=np.array([1, 2]))
+        return
+
+    def test__interpolate(self):
+        Z = np.array([1, 2, 3, 4])
+        ages = np.array([1e2, 1e3, 1e4])
+        grid = np.random.random((len(Z), len(ages), 1))
+        interpolator = self.__class__.tested_class(
+            grid, metallicities=Z, ages=ages
+        )
+
+        # single value
+        res = interpolator._interpolate(2., 1e3)
+        self.assertEqual(
+            (1,), res.shape
+        )
+        self.assertAlmostEqual(
+            grid[1, 1, 0], res
+        )
+        # in between points
+        res = interpolator._interpolate(1.5, 550.)
+        self.assertEqual(
+            (1,), res.shape
+        )
+        # Delaunay triangulation is not unique for regular grid, can get one of
+        # two possible answers:
+        opt1 = 0.5*(
+            grid[0, 0, 0] +
+            grid[1, 1, 0]
+        )
+        opt2 = 0.5*(
+            grid[1, 0, 0] +
+            grid[0, 1, 0]
+        )
+        self.assertTrue(
+            np.allclose(opt1, res) or np.allclose(opt2, res)
+        )
+
+        # array of values
+        res = interpolator._interpolate(
+            np.array([2., 3.]), np.array([1e2, 1e3])
+        )
+        self.assertEqual(
+            (2, 1), res.shape
+        )
+        self.assertAlmostEqual(
+            grid[1, 0, 0], res[0, 0]
+        )
+        self.assertAlmostEqual(
+            grid[2, 1, 0], res[1, 0]
+        )
+
+        # array valued qtty
+        grid = np.random.random((len(Z), len(ages), 3))
+        interpolator = self.__class__.tested_class(
+            grid, metallicities=Z, ages=ages
+        )
+        res = interpolator._interpolate(2., 1e3)
+        self.assertEqual(
+            (3,), res.shape
+        )
+        self.assertTrue(
+            np.allclose(grid[1, 1, :], res)
+        )
+        return
+
+
+class TestGridInterpolatorMassScaled(TestGridInterpolator):
+
+    @classmethod
+    def setUpClass(cls):
+        cls.tested_class = interpolators.GridInterpolatorMassScaled
+        return
+
+    def test__interpolate_with_masses(self):
+        Z = np.array([1, 2, 3, 4])
+        ages = np.array([1e2, 1e3, 1e4])
+        grid = np.random.random((len(Z), len(ages), 1))
+        interpolator = self.__class__.tested_class(
+            grid, metallicities=Z, ages=ages
+        )
+
+        # single value
+        res = interpolator._interpolate(2., 1e3, masses=2)
+        self.assertEqual(
+            (1,), res.shape
+        )
+        self.assertAlmostEqual(
+            2*grid[1, 1, 0], res
+        )
+
+        # array of values
+        res = interpolator._interpolate(
+            np.array([2., 3.]), np.array([1e2, 1e3]), masses=np.array([1, 2])
+        )
+        self.assertEqual(
+            (2, 1), res.shape
+        )
+        self.assertAlmostEqual(
+            grid[1, 0, 0], res[0, 0]
+        )
+        self.assertAlmostEqual(
+            2*grid[2, 1, 0], res[1, 0]
+        )
+
+        # array valued qtty
+        grid = np.random.random((len(Z), len(ages), 3))
+        interpolator = self.__class__.tested_class(
+            grid, metallicities=Z, ages=ages
+        )
+        res = interpolator._interpolate(2., 1e3, masses=2.5)
+        self.assertEqual(
+            (3,), res.shape
+        )
+        self.assertTrue(
+            np.allclose(2.5*grid[1, 1, :], res)
+        )
         return

From af8cc1a7b2fcbcd3e2dc544c415ebfb7c7b0a32b Mon Sep 17 00:00:00 2001
From: Martin Glatzle <mglatzle@mpa-garching.mpg.de>
Date: Tue, 25 Aug 2020 17:33:48 +0200
Subject: [PATCH 11/26] More verbose exceptions.

---
 hoki/interpolators.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/hoki/interpolators.py b/hoki/interpolators.py
index 11881447..4e642e59 100644
--- a/hoki/interpolators.py
+++ b/hoki/interpolators.py
@@ -51,11 +51,13 @@ def __init__(self, grid,
                 raise ValueError("Wrong dimensionality of input arrays.")
         if grid.shape[0] != len(metallicities):
             raise ValueError(
-                "Shapes of `grid` and `metallicities` are incompatible."
+                f"Shapes of `grid` {grid.shape} and "
+                f"`metallicities` {metallicities.shape} are incompatible."
             )
         if grid.shape[1] != len(ages):
             raise ValueError(
-                "Shapes of `grid` and `ages` are incompatible."
+                f"Shapes of `grid` {grid.shape} and "
+                f"`ages` {ages.shape} are incompatible."
             )
         self._dtype = dtype
         self._metallicities = metallicities.astype(self._dtype, copy=True)

From ee6150fd2e3549abc08d7e57a9e2df4e7179309a Mon Sep 17 00:00:00 2001
From: Martin Glatzle <mglatzle@mpa-garching.mpg.de>
Date: Tue, 25 Aug 2020 17:47:57 +0200
Subject: [PATCH 12/26] Initial implementation of SpectraInterpolator.

---
 hoki/constants.py                |  2 +
 hoki/interpolators.py            | 95 +++++++++++++++++++++++++++++++-
 hoki/tests/test_interpolators.py | 70 ++++++++++++++++++++++-
 3 files changed, 164 insertions(+), 3 deletions(-)

diff --git a/hoki/constants.py b/hoki/constants.py
index c640cc97..43cb9ceb 100644
--- a/hoki/constants.py
+++ b/hoki/constants.py
@@ -32,6 +32,8 @@
               "imf135_100", "imf135_300", "imf135all_100", "imf170_100",
               "imf170_300"]
 
+# wavelengths at which spectra are given [angstrom]
+BPASS_WAVELENGTHS = np.arange(1, 10001)
 
 dummy_dict = {'timestep': 0, 'age': 1, 'log(R1)': 2, 'log(T1)': 3, 'log(L1)': 4, 'M1': 5, 'He_core1': 6, 'CO_core1': 7,
               'ONe_core1': 8, 'X': 10, 'Y': 11, 'C': 14, 'N': 15, 'O': 16, 'Ne': 17, 'MH1': 18, 'MHe1': 19, 'MC1': 20,
diff --git a/hoki/interpolators.py b/hoki/interpolators.py
index 4e642e59..c4fb0915 100644
--- a/hoki/interpolators.py
+++ b/hoki/interpolators.py
@@ -6,7 +6,10 @@
 
 import numpy as np
 from scipy import interpolate
-from hoki.constants import BPASS_NUM_METALLICITIES, BPASS_TIME_BINS
+from hoki.constants import (
+    BPASS_NUM_METALLICITIES, BPASS_TIME_BINS, BPASS_WAVELENGTHS
+)
+from hoki import load
 import warnings
 
 
@@ -231,3 +234,93 @@ def _check_masses(masses):
             return masses[:, None]
         else:
             raise ValueError("Wrong dimensionality of `masses`.")
+
+
+class SpectraInterpolator(GridInterpolatorMassScaled):
+    """
+    Interpolate BPASS SSP spectra on a metallicity-age grid.
+
+    Interpolate a spectrum grid for single stellar populations (SSPs) with
+    fixed IMF provided by BPASS over its metallicity-age grid. The wavelength
+    range can be limited to something smaller than the BPASS default to reduce
+    memory footprint.
+
+    Parameters
+    ----------
+    data_path : `str`
+        The path to the folder containing the BPASS spectra. See
+        `load.all_spectra`.
+    imf : `str`
+        BPASS Identifier of the IMF to be used. See `load.all_spectra`.
+    binary : `bool`, optional
+        Use spectra including binaries or only single stars. Defaults to
+        `True`.
+    lam_min : `float`, optional
+        Limit the wavelength range on which this instance will perform
+        interpolation. Defaults to `None`, using full range available.
+    lam_max : `float`, optional
+        Limit the wavelength range on which this instance will perform
+        interpolation. Defaults to `None`, using full range available.
+    dtype : `type`, optional
+        The data type to be used by an instance of this class. Defaults to
+        `numpy.float64`. Can be used to reduce memory footprint.
+    """
+
+    def __init__(self, data_path, imf, binary=True,
+                 lam_min=None, lam_max=None,
+                 dtype=np.float64):
+        if lam_min is not None and lam_max is not None:
+            if lam_min >= lam_max:
+                raise ValueError("lam_min is larger than or equal to lam_max!")
+
+        lam = BPASS_WAVELENGTHS
+        if lam_min is not None:
+            idx_min = np.searchsorted(lam, lam_min, side='left')
+        else:
+            idx_min = None
+        if lam_max is not None:
+            idx_max = np.searchsorted(lam, lam_max, side='right')
+        else:
+            idx_max = None
+        self._wavelengths = lam[idx_min:idx_max].astype(
+            dtype, copy=True)
+
+        self._spectra = load.all_spectra(
+            data_path, imf, binary=binary)[:, :, idx_min:idx_max].astype(
+                dtype, copy=True
+            )
+
+        if len(self._wavelengths) != self._spectra.shape[2]:
+            raise ValueError(
+                "Incompatible dimesions for wavelengths "
+                f"{self._wavelengths.shape} and spectra {self._spectra.shape}."
+            )
+        super().__init__(self._spectra, dtype=dtype)
+
+        return
+
+    def interpolate(self, metallicities, ages, masses=1):
+        """
+        Perform interpolation on this instance's spectrum grid.
+
+        Parameters
+        ----------
+        metallicities : `numpy.ndarray` (N,)
+            Absolute initial stellar metallicities at which to
+            interpolate.
+        ages : `numpy.ndarray` (N,)
+            Stellar ages (in log scale) at which to interpolate.
+        masses : `numpy.ndarray` (N,) or `float`, optional
+            Stellar population masses in units of 1e6 M_\\odot. Used to scale
+            the interpolation result. Defaults to unity.
+
+        Returns
+        -------
+         : `numpy.ndarray` (N_lam,)
+            The wavelengths [angstrom] at which interpolated spectra are
+            provided.
+         : `numpy.ndarray` (N, N_lam)
+            Interpolated SEDs [L_\\odot/angstrom].
+        """
+        return self._wavelengths, \
+            self._interpolate(metallicities, ages, masses)
diff --git a/hoki/tests/test_interpolators.py b/hoki/tests/test_interpolators.py
index 8248f806..8111e322 100644
--- a/hoki/tests/test_interpolators.py
+++ b/hoki/tests/test_interpolators.py
@@ -1,6 +1,13 @@
-from unittest import TestCase
+"""
+Tests for the interpolators module.
+
+Author: Martin Glatzle
+"""
+from unittest import TestCase, mock
 from hoki import interpolators
-from hoki.constants import BPASS_NUM_METALLICITIES, BPASS_TIME_BINS
+from hoki.constants import (
+    BPASS_NUM_METALLICITIES, BPASS_TIME_BINS
+)
 import numpy as np
 
 
@@ -145,3 +152,62 @@ def test__interpolate_with_masses(self):
             np.allclose(2.5*grid[1, 1, :], res)
         )
         return
+
+
+class TestSpectraInterpolator(TestCase):
+
+    def setUp(self):
+        self.spectra = np.random.random((4, 4, 100))
+        self.lam = np.linspace(1, 10, num=100)
+        self.metallicities = np.linspace(1, 10, num=4)
+        self.ages = np.linspace(1, 10, num=4)
+
+        self.mock_lam = mock.patch.multiple(
+            'hoki.interpolators',
+            BPASS_WAVELENGTHS=self.lam,
+        )
+        self.mock_constructor_defaults = mock.patch(
+            'hoki.interpolators.GridInterpolator.__init__.__defaults__',
+            (self.metallicities, self.ages, np.float64)
+        )
+        self.mock_all_spectra = mock.patch(
+            'hoki.load.all_spectra',
+            return_value=self.spectra
+        )
+
+        return
+
+    def test_init(self):
+        with self.mock_lam, self.mock_constructor_defaults, \
+             self.mock_all_spectra:
+            # standard
+            interpolators.SpectraInterpolator(
+                '', '',
+            )
+            # limit lam range
+            interpolators.SpectraInterpolator(
+                '', '',
+                lam_min=3., lam_max=5.,
+            )
+        return
+
+    def test_interpolate(self):
+        with self.mock_lam, self.mock_constructor_defaults, \
+             self.mock_all_spectra:
+            interp = interpolators.SpectraInterpolator(
+                '', '',
+            )
+        res = interp.interpolate(1., 1.)
+        self.assertEqual(
+            len(res[0]), len(res[1])
+        )
+        interp.interpolate(1., 1., 1.)
+        res = interp.interpolate(
+            np.array([1., 2.]),
+            np.array([1., 2.]),
+            np.array([1., 2.])
+        )
+        self.assertEqual(
+            len(res[0]), res[1].shape[1]
+        )
+        return

From dd65f4f3b6824afcb68d77bab9a81e3bb7e3e3ff Mon Sep 17 00:00:00 2001
From: Martin Glatzle <mglatzle@mpa-garching.mpg.de>
Date: Fri, 28 Aug 2020 11:24:52 +0200
Subject: [PATCH 13/26] Fix patch.

---
 hoki/tests/test_interpolators.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/hoki/tests/test_interpolators.py b/hoki/tests/test_interpolators.py
index 8111e322..bc49b355 100644
--- a/hoki/tests/test_interpolators.py
+++ b/hoki/tests/test_interpolators.py
@@ -162,9 +162,9 @@ def setUp(self):
         self.metallicities = np.linspace(1, 10, num=4)
         self.ages = np.linspace(1, 10, num=4)
 
-        self.mock_lam = mock.patch.multiple(
-            'hoki.interpolators',
-            BPASS_WAVELENGTHS=self.lam,
+        self.mock_lam = mock.patch(
+            'hoki.interpolators.BPASS_WAVELENGTHS',
+            self.lam,
         )
         self.mock_constructor_defaults = mock.patch(
             'hoki.interpolators.GridInterpolator.__init__.__defaults__',

From d83c643edeba21e246aeaa1553bfa15afd23ce6d Mon Sep 17 00:00:00 2001
From: Martin Glatzle <mglatzle@mpa-garching.mpg.de>
Date: Fri, 28 Aug 2020 18:17:40 +0200
Subject: [PATCH 14/26] Tests for spectra interpolator.

---
 hoki/tests/test_interpolators.py | 41 ++++++++++++++++++++++++++++----
 1 file changed, 36 insertions(+), 5 deletions(-)

diff --git a/hoki/tests/test_interpolators.py b/hoki/tests/test_interpolators.py
index bc49b355..96394a95 100644
--- a/hoki/tests/test_interpolators.py
+++ b/hoki/tests/test_interpolators.py
@@ -6,7 +6,7 @@
 from unittest import TestCase, mock
 from hoki import interpolators
 from hoki.constants import (
-    BPASS_NUM_METALLICITIES, BPASS_TIME_BINS
+    BPASS_NUM_METALLICITIES, BPASS_TIME_BINS, BPASS_WAVELENGTHS
 )
 import numpy as np
 
@@ -159,6 +159,7 @@ class TestSpectraInterpolator(TestCase):
     def setUp(self):
         self.spectra = np.random.random((4, 4, 100))
         self.lam = np.linspace(1, 10, num=100)
+        self.lam2 = np.linspace(1, 10, num=10)
         self.metallicities = np.linspace(1, 10, num=4)
         self.ages = np.linspace(1, 10, num=4)
 
@@ -166,6 +167,10 @@ def setUp(self):
             'hoki.interpolators.BPASS_WAVELENGTHS',
             self.lam,
         )
+        self.mock_lam2 = mock.patch(
+            'hoki.interpolators.BPASS_WAVELENGTHS',
+            self.lam2,
+        )
         self.mock_constructor_defaults = mock.patch(
             'hoki.interpolators.GridInterpolator.__init__.__defaults__',
             (self.metallicities, self.ages, np.float64)
@@ -174,7 +179,6 @@ def setUp(self):
             'hoki.load.all_spectra',
             return_value=self.spectra
         )
-
         return
 
     def test_init(self):
@@ -191,23 +195,50 @@ def test_init(self):
             )
         return
 
+    def test_exceptions(self):
+        with self.mock_lam, self.mock_constructor_defaults, \
+             self.mock_all_spectra:
+            with self.assertRaises(ValueError):
+                interpolators.SpectraInterpolator(
+                    '', '', lam_min=3, lam_max=2,
+                )
+        with self.mock_lam2, self.mock_constructor_defaults, \
+             self.mock_all_spectra:
+            with self.assertRaises(ValueError):
+                interpolators.SpectraInterpolator(
+                    '', '',
+                )
+        return
+
     def test_interpolate(self):
         with self.mock_lam, self.mock_constructor_defaults, \
              self.mock_all_spectra:
             interp = interpolators.SpectraInterpolator(
                 '', '',
             )
+
+        # simple case
         res = interp.interpolate(1., 1.)
         self.assertEqual(
             len(res[0]), len(res[1])
         )
-        interp.interpolate(1., 1., 1.)
+
+        # with mass value
+        res = interp.interpolate(1., 1., 1.)
+
+        # multiple values
         res = interp.interpolate(
-            np.array([1., 2.]),
-            np.array([1., 2.]),
+            np.array([3., 3.]),
+            np.array([2., 2.]),
             np.array([1., 2.])
         )
         self.assertEqual(
             len(res[0]), res[1].shape[1]
         )
+        self.assertEqual(
+            2, res[1].shape[0]
+        )
+        self.assertTrue(
+            np.allclose(2*res[1][0, :], res[1][1, :])
+        )
         return

From a5850c732ef6523f675b6cbfbb72a3f31daba0a9 Mon Sep 17 00:00:00 2001
From: Martin Glatzle <mglatzle@mpa-garching.mpg.de>
Date: Fri, 28 Aug 2020 18:18:10 +0200
Subject: [PATCH 15/26] BUG: fix wavelenghts.

---
 hoki/constants.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/hoki/constants.py b/hoki/constants.py
index 43cb9ceb..bb2da1f8 100644
--- a/hoki/constants.py
+++ b/hoki/constants.py
@@ -33,7 +33,7 @@
               "imf170_300"]
 
 # wavelengths at which spectra are given [angstrom]
-BPASS_WAVELENGTHS = np.arange(1, 10001)
+BPASS_WAVELENGTHS = np.arange(1, 100001)
 
 dummy_dict = {'timestep': 0, 'age': 1, 'log(R1)': 2, 'log(T1)': 3, 'log(L1)': 4, 'M1': 5, 'He_core1': 6, 'CO_core1': 7,
               'ONe_core1': 8, 'X': 10, 'Y': 11, 'C': 14, 'N': 15, 'O': 16, 'Ne': 17, 'MH1': 18, 'MHe1': 19, 'MC1': 20,

From 4076f4426f4c665c91a0d97b7b295caf4f474294 Mon Sep 17 00:00:00 2001
From: Martin Glatzle <mglatzle@mpa-garching.mpg.de>
Date: Tue, 6 Oct 2020 16:46:16 +0200
Subject: [PATCH 16/26] Fix for all_spectra rename.

---
 hoki/interpolators.py            | 6 +++---
 hoki/tests/test_interpolators.py | 2 +-
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/hoki/interpolators.py b/hoki/interpolators.py
index c4fb0915..1a3e274b 100644
--- a/hoki/interpolators.py
+++ b/hoki/interpolators.py
@@ -249,9 +249,9 @@ class SpectraInterpolator(GridInterpolatorMassScaled):
     ----------
     data_path : `str`
         The path to the folder containing the BPASS spectra. See
-        `load.all_spectra`.
+        `load.spectra_all_z`.
     imf : `str`
-        BPASS Identifier of the IMF to be used. See `load.all_spectra`.
+        BPASS Identifier of the IMF to be used. See `load.spectra_all_z`.
     binary : `bool`, optional
         Use spectra including binaries or only single stars. Defaults to
         `True`.
@@ -285,7 +285,7 @@ def __init__(self, data_path, imf, binary=True,
         self._wavelengths = lam[idx_min:idx_max].astype(
             dtype, copy=True)
 
-        self._spectra = load.all_spectra(
+        self._spectra = load.spectra_all_z(
             data_path, imf, binary=binary)[:, :, idx_min:idx_max].astype(
                 dtype, copy=True
             )
diff --git a/hoki/tests/test_interpolators.py b/hoki/tests/test_interpolators.py
index 96394a95..67fcd333 100644
--- a/hoki/tests/test_interpolators.py
+++ b/hoki/tests/test_interpolators.py
@@ -176,7 +176,7 @@ def setUp(self):
             (self.metallicities, self.ages, np.float64)
         )
         self.mock_all_spectra = mock.patch(
-            'hoki.load.all_spectra',
+            'hoki.load.spectra_all_z',
             return_value=self.spectra
         )
         return

From 3e735c47ca87737846890b12cd02ea3a33eda2a5 Mon Sep 17 00:00:00 2001
From: Martin Glatzle <mglatzle@mpa-garching.mpg.de>
Date: Thu, 28 Jan 2021 15:28:02 +0100
Subject: [PATCH 17/26] Update docstrings.

---
 hoki/interpolators.py | 11 +++++------
 1 file changed, 5 insertions(+), 6 deletions(-)

diff --git a/hoki/interpolators.py b/hoki/interpolators.py
index 1a3e274b..a26a2ae9 100644
--- a/hoki/interpolators.py
+++ b/hoki/interpolators.py
@@ -33,7 +33,7 @@ class GridInterpolator():
         array of BPASS output metallicities.
     ages : `numpy.ndarray` (N_a,), optional
         The ages at which `grid` is evaluated. Defaults to the full array of
-        BPASS output ages in log scale.
+        BPASS output ages [yr] in log scale.
     dtype : `type`, optional
         The data type to be used by an instance of this class. Defaults to
         `numpy.float64`. Can be used to reduce memory footprint.
@@ -173,7 +173,7 @@ class GridInterpolatorMassScaled(GridInterpolator):
         array of BPASS output metallicities.
     ages : `numpy.ndarray` (N_a,), optional
         The ages at which `grid` is evaluated. Defaults to the full array of
-        BPASS output ages in log scale.
+        BPASS output ages [yr] in log scale.
     dtype : `type`, optional
         The data type to be used by an instance of this class. Defaults to
         `numpy.float64`. Can be used to reduce memory footprint.
@@ -248,10 +248,9 @@ class SpectraInterpolator(GridInterpolatorMassScaled):
     Parameters
     ----------
     data_path : `str`
-        The path to the folder containing the BPASS spectra. See
-        `load.spectra_all_z`.
+        The path to the folder containing the BPASS spectra files.
     imf : `str`
-        BPASS Identifier of the IMF to be used. See `load.spectra_all_z`.
+        BPASS Identifier of the IMF to be used, e.g. `"imf_chab100"`.
     binary : `bool`, optional
         Use spectra including binaries or only single stars. Defaults to
         `True`.
@@ -288,7 +287,7 @@ def __init__(self, data_path, imf, binary=True,
         self._spectra = load.spectra_all_z(
             data_path, imf, binary=binary)[:, :, idx_min:idx_max].astype(
                 dtype, copy=True
-            )
+        )
 
         if len(self._wavelengths) != self._spectra.shape[2]:
             raise ValueError(

From 21470cea68a7d000d7123c0b79970d03bdf265b1 Mon Sep 17 00:00:00 2001
From: Martin Glatzle <mglatzle@mpa-garching.mpg.de>
Date: Thu, 28 Jan 2021 15:33:13 +0100
Subject: [PATCH 18/26] Fix unused import.

---
 hoki/tests/test_interpolators.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/hoki/tests/test_interpolators.py b/hoki/tests/test_interpolators.py
index 67fcd333..40a2c1f4 100644
--- a/hoki/tests/test_interpolators.py
+++ b/hoki/tests/test_interpolators.py
@@ -6,7 +6,7 @@
 from unittest import TestCase, mock
 from hoki import interpolators
 from hoki.constants import (
-    BPASS_NUM_METALLICITIES, BPASS_TIME_BINS, BPASS_WAVELENGTHS
+    BPASS_NUM_METALLICITIES, BPASS_TIME_BINS
 )
 import numpy as np
 

From a142265738c861dd78302bca911e955ccf555543 Mon Sep 17 00:00:00 2001
From: Martin Glatzle <mglatzle@mpa-garching.mpg.de>
Date: Fri, 29 Jan 2021 11:54:50 +0100
Subject: [PATCH 19/26] First version of emissivity interpolator.

---
 hoki/interpolators.py            | 65 ++++++++++++++++++++++++++++++++
 hoki/tests/test_interpolators.py | 51 +++++++++++++++++++++++++
 2 files changed, 116 insertions(+)

diff --git a/hoki/interpolators.py b/hoki/interpolators.py
index a26a2ae9..33153c3a 100644
--- a/hoki/interpolators.py
+++ b/hoki/interpolators.py
@@ -323,3 +323,68 @@ def interpolate(self, metallicities, ages, masses=1):
         """
         return self._wavelengths, \
             self._interpolate(metallicities, ages, masses)
+
+
+class EmissivitiesInterpolator(GridInterpolatorMassScaled):
+    """
+    Interpolate BPASS SSP emissivities on a metallicity-age grid.
+
+    Interpolate a grid of emissivities for single stellar populations (SSPs)
+    with fixed IMF provided by BPASS over its metallicity-age grid.
+
+    Parameters
+    ----------
+    data_path : `str`
+        The path to the folder containing the BPASS emissivity files.
+    imf : `str`
+        BPASS Identifier of the IMF to be used, e.g. `"imf_chab100"`.
+    binary : `bool`, optional
+        Use spectra including binaries or only single stars. Defaults to
+        `True`.
+    dtype : `type`, optional
+        The data type to be used by an instance of this class. Defaults to
+        `numpy.float64`. Can be used to reduce memory footprint.
+    """
+
+    def __init__(self, data_path, imf, binary=True,
+                 dtype=np.float64):
+
+        self._emissivities = load.emissivities_all_z(
+            data_path, imf, binary=binary
+        ).astype(dtype, copy=True)
+
+        super().__init__(self._emissivities, dtype=dtype)
+
+        return
+
+    def interpolate(self, metallicities, ages, masses=1):
+        """
+        Perform interpolation on this instance's emissivities grid.
+
+        Parameters
+        ----------
+        metallicities : `numpy.ndarray` (N,)
+            Absolute initial stellar metallicities at which to
+            interpolate.
+        ages : `numpy.ndarray` (N,)
+            Stellar ages (in log scale) at which to interpolate.
+        masses : `numpy.ndarray` (N,) or `float`, optional
+            Stellar population masses in units of 1e6 M_\\odot. Used to scale
+            the interpolation result. Defaults to unity.
+
+        Returns
+        -------
+         : `numpy.ndarray` (N, 4)
+            Interpolated emissivities:
+
+            Nion in 1/s:
+                ionizing photon production rate
+            L_Halpha in ergs/s:
+                Balmer H line luminosity, assuming =log(Nion/s)-11.87
+            L_FUV in ergs/s/A:
+                luminosity in the FUV band (mean flux from 1556 to 1576A)
+            L_NUV in ergs/s/A:
+                luminosity in the NUV band (mean flux from 2257 to 2277A)
+
+        """
+        return self._interpolate(metallicities, ages, masses)
diff --git a/hoki/tests/test_interpolators.py b/hoki/tests/test_interpolators.py
index 40a2c1f4..338ea53e 100644
--- a/hoki/tests/test_interpolators.py
+++ b/hoki/tests/test_interpolators.py
@@ -242,3 +242,54 @@ def test_interpolate(self):
             np.allclose(2*res[1][0, :], res[1][1, :])
         )
         return
+
+
+class TestEmissivitiesInterpolator(TestCase):
+
+    def setUp(self):
+        # simulate 4 emssivities evaluated at 5 metallicity values and 10 age
+        # values
+        self.emissivities = np.random.random((5, 10, 4))
+        self.metallicities = np.linspace(1, 10, num=5)
+        self.ages = np.linspace(1, 10, num=10)
+
+        self.mock_constructor_defaults = mock.patch(
+            'hoki.interpolators.GridInterpolator.__init__.__defaults__',
+            (self.metallicities, self.ages, np.float64)
+        )
+        self.mock_all_emissivities = mock.patch(
+            'hoki.load.emissivities_all_z',
+            return_value=self.emissivities
+        )
+        return
+
+    def test_init(self):
+        with self.mock_constructor_defaults, self.mock_all_emissivities:
+            # standard
+            interpolators.EmissivitiesInterpolator('', '')
+        return
+
+    def test_interpolate(self):
+        with self.mock_constructor_defaults, self.mock_all_emissivities:
+            interp = interpolators.EmissivitiesInterpolator(
+                '', '',
+            )
+
+        # simple case
+        res = interp.interpolate(1., 1.)
+        self.assertEqual(4, len(res))
+
+        # with mass value
+        res = interp.interpolate(1., 1., 1.)
+        self.assertEqual(4, len(res))
+
+        # multiple values
+        res = interp.interpolate(
+            np.array([3., 3.]),
+            np.array([2., 2.]),
+            np.array([1., 2.])
+        )
+        self.assertTrue(
+            np.allclose(2*res[0, :], res[1, :])
+        )
+        return

From 6237f279c522bed94e0e3bd9d52a0b7aa0c9ed56 Mon Sep 17 00:00:00 2001
From: Martin Glatzle <mglatzle@mpa-garching.mpg.de>
Date: Fri, 29 Jan 2021 12:09:43 +0100
Subject: [PATCH 20/26] Make interpolators callable.

---
 hoki/interpolators.py            |  4 ++--
 hoki/tests/test_interpolators.py | 12 ++++++------
 2 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/hoki/interpolators.py b/hoki/interpolators.py
index 33153c3a..864f2062 100644
--- a/hoki/interpolators.py
+++ b/hoki/interpolators.py
@@ -298,7 +298,7 @@ def __init__(self, data_path, imf, binary=True,
 
         return
 
-    def interpolate(self, metallicities, ages, masses=1):
+    def __call__(self, metallicities, ages, masses=1):
         """
         Perform interpolation on this instance's spectrum grid.
 
@@ -357,7 +357,7 @@ def __init__(self, data_path, imf, binary=True,
 
         return
 
-    def interpolate(self, metallicities, ages, masses=1):
+    def __call__(self, metallicities, ages, masses=1):
         """
         Perform interpolation on this instance's emissivities grid.
 
diff --git a/hoki/tests/test_interpolators.py b/hoki/tests/test_interpolators.py
index 338ea53e..3d1c687a 100644
--- a/hoki/tests/test_interpolators.py
+++ b/hoki/tests/test_interpolators.py
@@ -218,16 +218,16 @@ def test_interpolate(self):
             )
 
         # simple case
-        res = interp.interpolate(1., 1.)
+        res = interp(1., 1.)
         self.assertEqual(
             len(res[0]), len(res[1])
         )
 
         # with mass value
-        res = interp.interpolate(1., 1., 1.)
+        res = interp(1., 1., 1.)
 
         # multiple values
-        res = interp.interpolate(
+        res = interp(
             np.array([3., 3.]),
             np.array([2., 2.]),
             np.array([1., 2.])
@@ -276,15 +276,15 @@ def test_interpolate(self):
             )
 
         # simple case
-        res = interp.interpolate(1., 1.)
+        res = interp(1., 1.)
         self.assertEqual(4, len(res))
 
         # with mass value
-        res = interp.interpolate(1., 1., 1.)
+        res = interp(1., 1., 1.)
         self.assertEqual(4, len(res))
 
         # multiple values
-        res = interp.interpolate(
+        res = interp(
             np.array([3., 3.]),
             np.array([2., 2.]),
             np.array([1., 2.])

From c6be523f9ae9d8d7b014c93d52e920d911dae931 Mon Sep 17 00:00:00 2001
From: Martin Glatzle <mglatzle@mpa-garching.mpg.de>
Date: Mon, 1 Feb 2021 14:00:11 +0100
Subject: [PATCH 21/26] Small docstring fix.

---
 hoki/interpolators.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/hoki/interpolators.py b/hoki/interpolators.py
index 864f2062..50e18833 100644
--- a/hoki/interpolators.py
+++ b/hoki/interpolators.py
@@ -308,7 +308,7 @@ def __call__(self, metallicities, ages, masses=1):
             Absolute initial stellar metallicities at which to
             interpolate.
         ages : `numpy.ndarray` (N,)
-            Stellar ages (in log scale) at which to interpolate.
+            Stellar ages [yr] in log scale at which to interpolate.
         masses : `numpy.ndarray` (N,) or `float`, optional
             Stellar population masses in units of 1e6 M_\\odot. Used to scale
             the interpolation result. Defaults to unity.
@@ -367,7 +367,7 @@ def __call__(self, metallicities, ages, masses=1):
             Absolute initial stellar metallicities at which to
             interpolate.
         ages : `numpy.ndarray` (N,)
-            Stellar ages (in log scale) at which to interpolate.
+            Stellar ages [yr] in log scale at which to interpolate.
         masses : `numpy.ndarray` (N,) or `float`, optional
             Stellar population masses in units of 1e6 M_\\odot. Used to scale
             the interpolation result. Defaults to unity.

From 8c9fe633f195ce5481acb52eb89e17eb086ebaac Mon Sep 17 00:00:00 2001
From: Martin Glatzle <mglatzle@mpa-garching.mpg.de>
Date: Thu, 4 Feb 2021 17:41:07 +0100
Subject: [PATCH 22/26] Fix docstring.

---
 hoki/interpolators.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/hoki/interpolators.py b/hoki/interpolators.py
index 50e18833..912932ec 100644
--- a/hoki/interpolators.py
+++ b/hoki/interpolators.py
@@ -375,7 +375,7 @@ def __call__(self, metallicities, ages, masses=1):
         Returns
         -------
          : `numpy.ndarray` (N, 4)
-            Interpolated emissivities:
+            Log of interpolated emissivities:
 
             Nion in 1/s:
                 ionizing photon production rate

From e735e7cf2dd0ebdac20bcf3811d179b85a679008 Mon Sep 17 00:00:00 2001
From: Martin Glatzle <mglatzle@mpa-garching.mpg.de>
Date: Thu, 4 Feb 2021 17:41:07 +0100
Subject: [PATCH 23/26] Fix docstring.

---
 hoki/interpolators.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/hoki/interpolators.py b/hoki/interpolators.py
index 50e18833..57119763 100644
--- a/hoki/interpolators.py
+++ b/hoki/interpolators.py
@@ -143,7 +143,7 @@ def _interpolate(self, metallicities, ages):
            np.amin(ages) < self._aMin:
             warnings.warn(
                 "Input ages for interpolation outside of available "
-                f"range {self._aMin} -- {self._aMax} provided. "
+                f"log range {self._aMin} -- {self._aMax} [yr] provided. "
                 "They will be clipped.",
                 UserWarning
             )
@@ -375,7 +375,7 @@ def __call__(self, metallicities, ages, masses=1):
         Returns
         -------
          : `numpy.ndarray` (N, 4)
-            Interpolated emissivities:
+            Log of interpolated emissivities:
 
             Nion in 1/s:
                 ionizing photon production rate

From 12d23e3f56099fdf60f4d50c1326026ddc40695f Mon Sep 17 00:00:00 2001
From: Martin Glatzle <mglatzle@mpa-garching.mpg.de>
Date: Mon, 8 Mar 2021 11:04:06 +0100
Subject: [PATCH 24/26] Comments/docstrings.

---
 hoki/interpolators.py            | 1 -
 hoki/tests/test_interpolators.py | 1 -
 2 files changed, 2 deletions(-)

diff --git a/hoki/interpolators.py b/hoki/interpolators.py
index 57119763..22c688f7 100644
--- a/hoki/interpolators.py
+++ b/hoki/interpolators.py
@@ -385,6 +385,5 @@ def __call__(self, metallicities, ages, masses=1):
                 luminosity in the FUV band (mean flux from 1556 to 1576A)
             L_NUV in ergs/s/A:
                 luminosity in the NUV band (mean flux from 2257 to 2277A)
-
         """
         return self._interpolate(metallicities, ages, masses)
diff --git a/hoki/tests/test_interpolators.py b/hoki/tests/test_interpolators.py
index 3d1c687a..bdb58d8d 100644
--- a/hoki/tests/test_interpolators.py
+++ b/hoki/tests/test_interpolators.py
@@ -265,7 +265,6 @@ def setUp(self):
 
     def test_init(self):
         with self.mock_constructor_defaults, self.mock_all_emissivities:
-            # standard
             interpolators.EmissivitiesInterpolator('', '')
         return
 

From c4bdbcf35ec02e1dffdc59cecd2f96c880568811 Mon Sep 17 00:00:00 2001
From: Martin Glatzle <mglatzle@mpa-garching.mpg.de>
Date: Mon, 8 Mar 2021 11:04:44 +0100
Subject: [PATCH 25/26] Convert emssivities to linear so that mass scaling
 works.

---
 hoki/interpolators.py            | 6 +++---
 hoki/tests/test_interpolators.py | 4 +++-
 2 files changed, 6 insertions(+), 4 deletions(-)

diff --git a/hoki/interpolators.py b/hoki/interpolators.py
index 22c688f7..42d44817 100644
--- a/hoki/interpolators.py
+++ b/hoki/interpolators.py
@@ -349,8 +349,8 @@ class EmissivitiesInterpolator(GridInterpolatorMassScaled):
     def __init__(self, data_path, imf, binary=True,
                  dtype=np.float64):
 
-        self._emissivities = load.emissivities_all_z(
-            data_path, imf, binary=binary
+        self._emissivities = 10**(
+            load.emissivities_all_z(data_path, imf, binary=binary)
         ).astype(dtype, copy=True)
 
         super().__init__(self._emissivities, dtype=dtype)
@@ -375,7 +375,7 @@ def __call__(self, metallicities, ages, masses=1):
         Returns
         -------
          : `numpy.ndarray` (N, 4)
-            Log of interpolated emissivities:
+            Interpolated emissivities:
 
             Nion in 1/s:
                 ionizing photon production rate
diff --git a/hoki/tests/test_interpolators.py b/hoki/tests/test_interpolators.py
index bdb58d8d..0b4e1f33 100644
--- a/hoki/tests/test_interpolators.py
+++ b/hoki/tests/test_interpolators.py
@@ -277,10 +277,12 @@ def test_interpolate(self):
         # simple case
         res = interp(1., 1.)
         self.assertEqual(4, len(res))
+        self.assertTrue(np.allclose(10**self.emissivities[0, 0, :], res))
 
         # with mass value
-        res = interp(1., 1., 1.)
+        res = interp(1., 1., 2.)
         self.assertEqual(4, len(res))
+        self.assertTrue(np.allclose(2*10**self.emissivities[0, 0, :], res))
 
         # multiple values
         res = interp(

From 7b0a1400d5660949b06589045a3965880ef60184 Mon Sep 17 00:00:00 2001
From: Martin Glatzle <findesgh@yandex.ru>
Date: Sat, 19 Feb 2022 16:26:04 +0100
Subject: [PATCH 26/26] Make docstrings more precise.

---
 hoki/interpolators.py | 15 ++++++++++-----
 1 file changed, 10 insertions(+), 5 deletions(-)

diff --git a/hoki/interpolators.py b/hoki/interpolators.py
index 42d44817..1d36f96f 100644
--- a/hoki/interpolators.py
+++ b/hoki/interpolators.py
@@ -95,6 +95,8 @@ def _interpolate(self, metallicities, ages):
         """
         Perform interpolation on this instance's grid.
 
+        Input parameters will be clipped to available range.
+
         Parameters
         ----------
         metallicities : `numpy.ndarray` (N,)
@@ -243,7 +245,8 @@ class SpectraInterpolator(GridInterpolatorMassScaled):
     Interpolate a spectrum grid for single stellar populations (SSPs) with
     fixed IMF provided by BPASS over its metallicity-age grid. The wavelength
     range can be limited to something smaller than the BPASS default to reduce
-    memory footprint.
+    memory footprint. Provided limits beyond the available range will be
+    clipped.
 
     Parameters
     ----------
@@ -255,11 +258,13 @@ class SpectraInterpolator(GridInterpolatorMassScaled):
         Use spectra including binaries or only single stars. Defaults to
         `True`.
     lam_min : `float`, optional
-        Limit the wavelength range on which this instance will perform
-        interpolation. Defaults to `None`, using full range available.
+        Lower limit of the wavelength range on which this instance will perform
+        interpolation. Defaults to `None`, using the minimal wavelength
+        available.
     lam_max : `float`, optional
-        Limit the wavelength range on which this instance will perform
-        interpolation. Defaults to `None`, using full range available.
+        Upper limit of the wavelength range on which this instance will perform
+        interpolation. Defaults to `None`, using the maximal wavelength
+        available.
     dtype : `type`, optional
         The data type to be used by an instance of this class. Defaults to
         `numpy.float64`. Can be used to reduce memory footprint.