Repository: skypyproject/skypy

Issue:
move spectrum and magnitude functions to photometry module
We have amassed a number of functions in `skypy.galaxy.spectrum` and `skypy.galaxy.luminosity` for working with spectra and magnitudes, but these functions are more generally useful and equally apply to the simulation of stars. I propose that we create a new `skypy.photometry` module where such functions can live.

Discussion:

Patch:
diff --git a/docs/galaxy.rst b/docs/galaxy.rst
index b480e9b..6d5ab66 100644
--- a/docs/galaxy.rst
+++ b/docs/galaxy.rst
@@ -68,8 +68,6 @@ The following models are found in the `skypy.galaxy.spectrum` package.
    :nosignatures:
 
    dirichlet_coefficients
-   mag_ab
-   SpectrumTemplates
    KCorrectTemplates
    kcorrect
 
diff --git a/docs/utils/index.rst b/docs/utils/index.rst
index 727d077..d1747b8 100644
--- a/docs/utils/index.rst
+++ b/docs/utils/index.rst
@@ -120,6 +120,20 @@ should be placed below decorators which modify any of the independent
 arguments.
 
 
+Photometry (`skypy.utils.photometry`)
+=====================================
+
+This module contains methods that model spectral energy distributions and
+calculate photometric properties.
+
+.. currentmodule:: skypy.utils.photometry
+.. autosummary::
+   :nosignatures:
+
+   mag_ab
+   SpectrumTemplates
+
+
 Random sampling (`skypy.utils.random`)
 ======================================
 
@@ -136,3 +150,4 @@ Reference/API
 =============
 
 .. automodapi:: skypy.utils
+.. automodapi:: skypy.utils.photometry
diff --git a/skypy/galaxy/spectrum.py b/skypy/galaxy/spectrum.py
index 1fff2b7..99bece3 100644
--- a/skypy/galaxy/spectrum.py
+++ b/skypy/galaxy/spectrum.py
@@ -2,28 +2,19 @@ r"""Galaxy spectrum module.
 
 """
 
-import numpy as np
 from astropy import units
 from astropy.io import fits
+import numpy as np
 from pkg_resources import resource_filename
-from abc import ABCMeta, abstractmethod
+from skypy.utils.photometry import SpectrumTemplates
 
 
 __all__ = [
     'dirichlet_coefficients',
-    'mag_ab',
-    'SpectrumTemplates',
     'KCorrectTemplates',
     'kcorrect',
 ]
 
-try:
-    import speclite.filters
-except ImportError:
-    HAS_SPECLITE = False
-else:
-    HAS_SPECLITE = True
-
 
 def dirichlet_coefficients(redshift, alpha0, alpha1, z1=1., weight=None):
     r"""Dirichlet-distributed SED coefficients.
@@ -118,212 +109,6 @@ def dirichlet_coefficients(redshift, alpha0, alpha1, z1=1., weight=None):
     return coeff
 
 
-def mag_ab(wavelength, spectrum, filters, *, redshift=None, coefficients=None,
-           distmod=None, interpolate=1000):
-    r'''Compute AB magnitudes from spectra and filters.
-
-    This function takes *emission* spectra and observation filters and computes
-    bandpass AB magnitudes [1]_.
-
-    The filter specification in the `filters` argument is passed unchanged to
-    `speclite.filters.load_filters`. See there for the syntax, and the list of
-    supported values.
-
-    The emission spectra can optionally be redshifted. If the `redshift`
-    parameter is given, the output array will have corresponding axes. If the
-    `interpolate` parameter is not `False`, at most that number of redshifted
-    spectra are computed, and the remainder is interpolated from the results.
-
-    The spectra can optionally be combined. If the `coefficients` parameter is
-    given, its shape must match `spectra`, and the corresponding axes are
-    contracted using a product sum. If the spectra are redshifted, the
-    coefficients array can contain axes for each redshift.
-
-    By default, absolute magnitudes are returned. To compute apparent magnitudes
-    instead, provide the `distmod` argument with the distance modulus for each
-    redshift. The distance modulus is applied after redshifts and coefficients
-    and should match the shape of the `redshift` array.
-
-    Parameters
-    ----------
-    wavelength : (nw,) `~astropy.units.Quantity` or array_like
-        Wavelength array of the spectrum.
-    spectrum : ([ns,] nw,) `~astropy.units.Quantity` or array_like
-        Emission spectrum. Can be multidimensional for computing with multiple
-        spectra of the same wavelengths. The last axis is the wavelength axis.
-    filters : str or list of str
-        Filter specification, loaded filters are array_like of shape (nf,).
-    redshift : (nz,) array_like, optional
-        Optional array of redshifts. Can be multidimensional.
-    coefficients : ([nz,] [ns,]) array_like
-        Optional coefficients for combining spectra. Axes must be compatible
-        with all redshift and spectrum dimensions.
-    distmod : (nz,) array_like, optional
-        Optional distance modulus for each redshift. Shape must be compatible
-        with redshift dimensions.
-    interpolate : int or `False`, optional
-        Maximum number of redshifts to compute explicitly. Default is `1000`.
-
-    Returns
-    -------
-    mag_ab : ([nz,] [ns,] nf,) array_like
-        The AB magnitude of each redshift (if given), each spectrum (if not
-        combined), and each filter.
-
-    Warnings
-    --------
-    The :mod:`speclite` package must be installed to use this function.
-
-    References
-    ----------
-    .. [1] M. R. Blanton et al., 2003, AJ, 125, 2348
-
-    '''
-
-    # load the filters
-    if np.ndim(filters) == 0:
-        filters = (filters,)
-    filters = speclite.filters.load_filters(*filters)
-    if np.shape(filters) == (1,):
-        filters = filters[0]
-
-    # number of dimensions for each input
-    nd_s = len(np.shape(spectrum)[:-1])  # last axis is spectral axis
-    nd_f = len(np.shape(filters))
-    nd_z = len(np.shape(redshift))
-
-    # check if interpolation is necessary
-    if interpolate and np.size(redshift) <= interpolate:
-        interpolate = False
-
-    # if interpolating, split the redshift range into `interpolate` bits
-    if interpolate:
-        redshift_ = np.linspace(np.min(redshift), np.max(redshift), interpolate)
-    else:
-        redshift_ = redshift if redshift is not None else 0
-
-    # working array shape
-    m_shape = np.shape(redshift_) + np.shape(spectrum)[:-1] + np.shape(filters)
-
-    # compute AB maggies for every redshift, spectrum, and filter
-    m = np.empty(m_shape)
-    for i, z in np.ndenumerate(redshift_):
-        for j, f in np.ndenumerate(filters):
-            # create a shifted filter for redshift
-            fs = f.create_shifted(z)
-            m[i+(...,)+j] = fs.get_ab_maggies(spectrum, wavelength)
-
-    # if interpolating, compute the full set of redshifts
-    if interpolate:
-        # diy interpolation keeps memory use to a minimum
-        dm = np.diff(m, axis=0, append=m[-1:])
-        u, n = np.modf(np.interp(redshift, redshift_, np.arange(redshift_.size)))
-        n = n.astype(int)
-        u = u.reshape(u.shape + (1,)*(nd_s+nd_f))
-        m = np.ascontiguousarray(m[n])
-        m += u*dm[n]
-        del(dm, n, u)
-
-    # combine spectra if asked to
-    if coefficients is not None:
-        # contraction over spectrum axes (`nd_z` to `nd_z+nd_s`)
-        if nd_z == 0:
-            m = np.matmul(coefficients, m)
-        else:
-            c = np.reshape(coefficients, np.shape(coefficients) + (1,)*nd_f)
-            m = np.sum(m*c, axis=tuple(range(nd_z, nd_z+nd_s)))
-        # no spectrum axes left
-        nd_s = 0
-
-    # convert maggies to magnitudes
-    np.log10(m, out=m)
-    m *= -2.5
-
-    # apply the redshift K-correction if necessary
-    if redshift is not None:
-        kcorr = -2.5*np.log10(1 + redshift)
-        m += np.reshape(kcorr, kcorr.shape + (1,)*(nd_s+nd_f))
-
-    # add distance modulus if given
-    if distmod is not None:
-        m += np.reshape(distmod, np.shape(distmod) + (1,)*(nd_s+nd_f))
-
-    # done
-    return m
-
-
-class SpectrumTemplates(metaclass=ABCMeta):
-    '''Base class for composite galaxy spectra from a set of basis templates'''
-
-    @abstractmethod
-    def __init__(self):
-        raise NotImplementedError
-
-    def absolute_magnitudes(self, coefficients, filters, stellar_mass=None):
-        '''Galaxy AB absolute magnitudes from template spectra.
-
-        This function calculates photometric AB absolute magnitudes for
-        galaxies whose spectra are modelled as a linear combination of a set of
-        template spectra.
-
-        Parameters
-        ----------
-        coefficients : (ng, nt) array_like
-            Array of spectrum coefficients.
-        filters : str or list of str
-            Bandpass filter specification for `~speclite.filters.load_filters`.
-        stellar_mass : (ng,) array_like, optional
-            Optional array of stellar masses for each galaxy in template units.
-
-        Returns
-        -------
-        magnitudes : (ng, nf) array_like
-            The absolute AB magnitude of each object in each filter, where
-            ``nf`` is the number of loaded filters.
-        '''
-        mass_modulus = -2.5*np.log10(stellar_mass) if stellar_mass is not None else 0
-        M = mag_ab(self.wavelength, self.templates, filters, coefficients=coefficients)
-        return (M.T + mass_modulus).T
-
-    def apparent_magnitudes(self, coefficients, redshift, filters, cosmology, *,
-                            stellar_mass=None, resolution=1000):
-        '''Galaxy AB apparent magnitudes from template spectra.
-
-        This function calculates photometric AB apparent magnitudes for
-        galaxies whose spectra are modelled as a linear combination of a set of
-        template spectra.
-
-        Parameters
-        ----------
-        coefficients : (ng, nt) array_like
-            Array of spectrum coefficients.
-        redshifts : (ng,) array_like
-            Array of redshifts for each galaxy used to calculte the distance
-            modulus and k-correction.
-        filters : str or list of str
-            Bandpass filter specification for `~speclite.filters.load_filters`.
-        cosmology : Cosmology
-            Astropy Cosmology object to calculate distance modulus.
-        stellar_mass : (ng,) array_like, optional
-            Optional array of stellar masses for each galaxy in template units.
-        resolution : integer, optional
-            Redshift resolution for intepolating magnitudes. Default is 1000. If
-            the number of objects is less than resolution their magnitudes are
-            calculated directly without interpolation.
-
-        Returns
-        -------
-        magnitudes : (ng, nf) array_like
-            The apparent AB magnitude of each object in each filter, where
-            ``nf`` is the number of loaded filters.
-        '''
-        distmod = cosmology.distmod(redshift).value
-        mass_modulus = -2.5*np.log10(stellar_mass) if stellar_mass is not None else 0
-        m = mag_ab(self.wavelength, self.templates, filters, redshift=redshift,
-                   coefficients=coefficients, distmod=distmod, interpolate=resolution)
-        return (m.T + mass_modulus).T
-
-
 class KCorrectTemplates(SpectrumTemplates):
     '''Galaxy spectra from kcorrect templates.
 
diff --git a/skypy/utils/__init__.py b/skypy/utils/__init__.py
index fba5334..62c729d 100644
--- a/skypy/utils/__init__.py
+++ b/skypy/utils/__init__.py
@@ -4,6 +4,7 @@ This module contains utility functions.
 
 __all__ = []
 
+from . import photometry
 from . import random
 from . import special
 
diff --git a/skypy/utils/photometry.py b/skypy/utils/photometry.py
new file mode 100644
index 0000000..2f3e7dc
--- /dev/null
+++ b/skypy/utils/photometry.py
@@ -0,0 +1,226 @@
+"""Photometry module.
+
+"""
+
+
+from abc import ABCMeta, abstractmethod
+import numpy as np
+
+
+__all__ = [
+    'mag_ab',
+    'SpectrumTemplates',
+]
+
+try:
+    import speclite.filters
+except ImportError:
+    HAS_SPECLITE = False
+else:
+    HAS_SPECLITE = True
+
+
+def mag_ab(wavelength, spectrum, filters, *, redshift=None, coefficients=None,
+           distmod=None, interpolate=1000):
+    r'''Compute AB magnitudes from spectra and filters.
+
+    This function takes *emission* spectra and observation filters and computes
+    bandpass AB magnitudes [1]_.
+
+    The filter specification in the `filters` argument is passed unchanged to
+    `speclite.filters.load_filters`. See there for the syntax, and the list of
+    supported values.
+
+    The emission spectra can optionally be redshifted. If the `redshift`
+    parameter is given, the output array will have corresponding axes. If the
+    `interpolate` parameter is not `False`, at most that number of redshifted
+    spectra are computed, and the remainder is interpolated from the results.
+
+    The spectra can optionally be combined. If the `coefficients` parameter is
+    given, its shape must match `spectra`, and the corresponding axes are
+    contracted using a product sum. If the spectra are redshifted, the
+    coefficients array can contain axes for each redshift.
+
+    By default, absolute magnitudes are returned. To compute apparent magnitudes
+    instead, provide the `distmod` argument with the distance modulus for each
+    redshift. The distance modulus is applied after redshifts and coefficients
+    and should match the shape of the `redshift` array.
+
+    Parameters
+    ----------
+    wavelength : (nw,) `~astropy.units.Quantity` or array_like
+        Wavelength array of the spectrum.
+    spectrum : ([ns,] nw,) `~astropy.units.Quantity` or array_like
+        Emission spectrum. Can be multidimensional for computing with multiple
+        spectra of the same wavelengths. The last axis is the wavelength axis.
+    filters : str or list of str
+        Filter specification, loaded filters are array_like of shape (nf,).
+    redshift : (nz,) array_like, optional
+        Optional array of redshifts. Can be multidimensional.
+    coefficients : ([nz,] [ns,]) array_like
+        Optional coefficients for combining spectra. Axes must be compatible
+        with all redshift and spectrum dimensions.
+    distmod : (nz,) array_like, optional
+        Optional distance modulus for each redshift. Shape must be compatible
+        with redshift dimensions.
+    interpolate : int or `False`, optional
+        Maximum number of redshifts to compute explicitly. Default is `1000`.
+
+    Returns
+    -------
+    mag_ab : ([nz,] [ns,] nf,) array_like
+        The AB magnitude of each redshift (if given), each spectrum (if not
+        combined), and each filter.
+
+    Warnings
+    --------
+    The :mod:`speclite` package must be installed to use this function.
+
+    References
+    ----------
+    .. [1] M. R. Blanton et al., 2003, AJ, 125, 2348
+
+    '''
+
+    # load the filters
+    if np.ndim(filters) == 0:
+        filters = (filters,)
+    filters = speclite.filters.load_filters(*filters)
+    if np.shape(filters) == (1,):
+        filters = filters[0]
+
+    # number of dimensions for each input
+    nd_s = len(np.shape(spectrum)[:-1])  # last axis is spectral axis
+    nd_f = len(np.shape(filters))
+    nd_z = len(np.shape(redshift))
+
+    # check if interpolation is necessary
+    if interpolate and np.size(redshift) <= interpolate:
+        interpolate = False
+
+    # if interpolating, split the redshift range into `interpolate` bits
+    if interpolate:
+        redshift_ = np.linspace(np.min(redshift), np.max(redshift), interpolate)
+    else:
+        redshift_ = redshift if redshift is not None else 0
+
+    # working array shape
+    m_shape = np.shape(redshift_) + np.shape(spectrum)[:-1] + np.shape(filters)
+
+    # compute AB maggies for every redshift, spectrum, and filter
+    m = np.empty(m_shape)
+    for i, z in np.ndenumerate(redshift_):
+        for j, f in np.ndenumerate(filters):
+            # create a shifted filter for redshift
+            fs = f.create_shifted(z)
+            m[i+(...,)+j] = fs.get_ab_maggies(spectrum, wavelength)
+
+    # if interpolating, compute the full set of redshifts
+    if interpolate:
+        # diy interpolation keeps memory use to a minimum
+        dm = np.diff(m, axis=0, append=m[-1:])
+        u, n = np.modf(np.interp(redshift, redshift_, np.arange(redshift_.size)))
+        n = n.astype(int)
+        u = u.reshape(u.shape + (1,)*(nd_s+nd_f))
+        m = np.ascontiguousarray(m[n])
+        m += u*dm[n]
+        del(dm, n, u)
+
+    # combine spectra if asked to
+    if coefficients is not None:
+        # contraction over spectrum axes (`nd_z` to `nd_z+nd_s`)
+        if nd_z == 0:
+            m = np.matmul(coefficients, m)
+        else:
+            c = np.reshape(coefficients, np.shape(coefficients) + (1,)*nd_f)
+            m = np.sum(m*c, axis=tuple(range(nd_z, nd_z+nd_s)))
+        # no spectrum axes left
+        nd_s = 0
+
+    # convert maggies to magnitudes
+    np.log10(m, out=m)
+    m *= -2.5
+
+    # apply the redshift K-correction if necessary
+    if redshift is not None:
+        kcorr = -2.5*np.log10(1 + redshift)
+        m += np.reshape(kcorr, kcorr.shape + (1,)*(nd_s+nd_f))
+
+    # add distance modulus if given
+    if distmod is not None:
+        m += np.reshape(distmod, np.shape(distmod) + (1,)*(nd_s+nd_f))
+
+    # done
+    return m
+
+
+class SpectrumTemplates(metaclass=ABCMeta):
+    '''Base class for composite galaxy spectra from a set of basis templates'''
+
+    @abstractmethod
+    def __init__(self):
+        raise NotImplementedError
+
+    def absolute_magnitudes(self, coefficients, filters, stellar_mass=None):
+        '''Galaxy AB absolute magnitudes from template spectra.
+
+        This function calculates photometric AB absolute magnitudes for
+        galaxies whose spectra are modelled as a linear combination of a set of
+        template spectra.
+
+        Parameters
+        ----------
+        coefficients : (ng, nt) array_like
+            Array of spectrum coefficients.
+        filters : str or list of str
+            Bandpass filter specification for `~speclite.filters.load_filters`.
+        stellar_mass : (ng,) array_like, optional
+            Optional array of stellar masses for each galaxy in template units.
+
+        Returns
+        -------
+        magnitudes : (ng, nf) array_like
+            The absolute AB magnitude of each object in each filter, where
+            ``nf`` is the number of loaded filters.
+        '''
+        mass_modulus = -2.5*np.log10(stellar_mass) if stellar_mass is not None else 0
+        M = mag_ab(self.wavelength, self.templates, filters, coefficients=coefficients)
+        return (M.T + mass_modulus).T
+
+    def apparent_magnitudes(self, coefficients, redshift, filters, cosmology, *,
+                            stellar_mass=None, resolution=1000):
+        '''Galaxy AB apparent magnitudes from template spectra.
+
+        This function calculates photometric AB apparent magnitudes for
+        galaxies whose spectra are modelled as a linear combination of a set of
+        template spectra.
+
+        Parameters
+        ----------
+        coefficients : (ng, nt) array_like
+            Array of spectrum coefficients.
+        redshifts : (ng,) array_like
+            Array of redshifts for each galaxy used to calculte the distance
+            modulus and k-correction.
+        filters : str or list of str
+            Bandpass filter specification for `~speclite.filters.load_filters`.
+        cosmology : Cosmology
+            Astropy Cosmology object to calculate distance modulus.
+        stellar_mass : (ng,) array_like, optional
+            Optional array of stellar masses for each galaxy in template units.
+        resolution : integer, optional
+            Redshift resolution for intepolating magnitudes. Default is 1000. If
+            the number of objects is less than resolution their magnitudes are
+            calculated directly without interpolation.
+
+        Returns
+        -------
+        magnitudes : (ng, nf) array_like
+            The apparent AB magnitude of each object in each filter, where
+            ``nf`` is the number of loaded filters.
+        '''
+        distmod = cosmology.distmod(redshift).value
+        mass_modulus = -2.5*np.log10(stellar_mass) if stellar_mass is not None else 0
+        m = mag_ab(self.wavelength, self.templates, filters, redshift=redshift,
+                   coefficients=coefficients, distmod=distmod, interpolate=resolution)
+        return (m.T + mass_modulus).T


Test patch:
diff --git a/skypy/galaxy/tests/test_spectrum.py b/skypy/galaxy/tests/test_spectrum.py
index 0862351..45338de 100644
--- a/skypy/galaxy/tests/test_spectrum.py
+++ b/skypy/galaxy/tests/test_spectrum.py
@@ -1,7 +1,7 @@
 import numpy as np
 import scipy.stats
 import pytest
-from skypy.galaxy.spectrum import HAS_SPECLITE
+from skypy.utils.photometry import HAS_SPECLITE
 
 
 @pytest.mark.flaky
@@ -64,144 +64,6 @@ def test_sampling_coefficients():
         dirichlet_coefficients(redshift, [2.5, 2.5], [2.5, 2.5], weight=[1, 2, 3])
 
 
-@pytest.mark.skipif(not HAS_SPECLITE, reason='test requires speclite')
-def test_mag_ab_standard_source():
-
-    from astropy import units
-    from speclite.filters import FilterResponse
-    from skypy.galaxy.spectrum import mag_ab
-
-    # create a filter
-    filt_lam = np.logspace(0, 4, 1000)*units.AA
-    filt_tx = np.exp(-((filt_lam - 1000*units.AA)/(100*units.AA))**2)
-    filt_tx[[0, -1]] = 0
-    FilterResponse(wavelength=filt_lam, response=filt_tx,
-                   meta=dict(group_name='test', band_name='filt'))
-
-    # test that the AB standard source has zero magnitude
-    lam = filt_lam  # same grid to prevent interpolation issues
-    flam = 0.10885464149979998*units.Unit('erg s-1 cm-2 AA')/lam**2
-
-    m = mag_ab(lam, flam, 'test-filt')
-
-    assert np.isclose(m, 0)
-
-
-@pytest.mark.skipif(not HAS_SPECLITE, reason='test requires speclite')
-def test_mag_ab_redshift_dependence():
-
-    from astropy import units
-    from speclite.filters import FilterResponse
-    from skypy.galaxy.spectrum import mag_ab
-
-    # make a wide tophat bandpass
-    filt_lam = [1.0e-10, 1.1e-10, 1.0e0, 0.9e10, 1.0e10]
-    filt_tx = [0., 1., 1., 1., 0.]
-    FilterResponse(wavelength=filt_lam, response=filt_tx,
-                   meta=dict(group_name='test', band_name='filt'))
-
-    # create a narrow gaussian source
-    lam = np.logspace(-11, 11, 1000)*units.AA
-    flam = np.exp(-((lam - 100*units.AA)/(10*units.AA))**2)*units.Unit('erg s-1 cm-2 AA-1')
-
-    # array of redshifts
-    z = np.linspace(0, 1, 11)
-
-    # compute the AB magnitude at different redshifts
-    m = mag_ab(lam, flam, 'test-filt', redshift=z)
-
-    # compare with expected redshift dependence
-    np.testing.assert_allclose(m, m[0] - 2.5*np.log10(1 + z))
-
-
-@pytest.mark.skipif(not HAS_SPECLITE, reason='test requires speclite')
-def test_mag_ab_multi():
-
-    from astropy import units
-    from skypy.galaxy.spectrum import mag_ab
-    from speclite.filters import FilterResponse
-
-    # 5 redshifts
-    z = np.linspace(0, 1, 5)
-
-    # 2 Gaussian bandpasses
-    filt_lam = np.logspace(0, 4, 1000) * units.AA
-    filt_mean = np.array([[1000], [2000]]) * units.AA
-    filt_width = np.array([[100], [10]]) * units.AA
-    filt_tx = np.exp(-((filt_lam-filt_mean)/filt_width)**2)
-    filt_tx[:, [0, -1]] = 0
-    FilterResponse(wavelength=filt_lam, response=filt_tx[0],
-                   meta=dict(group_name='test', band_name='filt0'))
-    FilterResponse(wavelength=filt_lam, response=filt_tx[1],
-                   meta=dict(group_name='test', band_name='filt1'))
-
-    # 3 Flat Spectra
-    # to prevent issues with interpolation, collect all redshifted filt_lam
-    lam = []
-    for z_ in z:
-        lam = np.union1d(lam, filt_lam.value/(1+z_))
-    lam = lam*filt_lam.unit
-    A = np.array([[2], [3], [4]])
-    flam = A * 0.10885464149979998*units.Unit('erg s-1 cm-2 AA')/lam**2
-
-    # Compare calculated magnitudes with truth
-    magnitudes = mag_ab(lam, flam, ['test-filt0', 'test-filt1'], redshift=z)
-    truth = -2.5 * np.log10(A * (1+z)).T[:, :, np.newaxis]
-    assert magnitudes.shape == (5, 3, 2)
-    np.testing.assert_allclose(*np.broadcast_arrays(magnitudes, truth), rtol=1e-4)
-
-
-@pytest.mark.skipif(not HAS_SPECLITE, reason='test requires speclite')
-def test_template_spectra():
-
-    from astropy import units
-    from skypy.galaxy.spectrum import mag_ab, SpectrumTemplates
-    from astropy.cosmology import Planck15
-    from speclite.filters import FilterResponse
-
-    class TestTemplates(SpectrumTemplates):
-        '''Three flat templates'''
-
-        def __init__(self):
-            self.wavelength = np.logspace(-1, 4, 1000)*units.AA
-            A = np.array([[2], [3], [4]]) * 0.10885464149979998
-            self.templates = A * units.Unit('erg s-1 cm-2 AA') / self.wavelength**2
-
-    test_templates = TestTemplates()
-    lam, flam = test_templates.wavelength, test_templates.templates
-
-    # Gaussian bandpass
-    filt_lam = np.logspace(0, 4, 1000)*units.AA
-    filt_tx = np.exp(-((filt_lam - 1000*units.AA)/(100*units.AA))**2)
-    filt_tx[[0, -1]] = 0
-    FilterResponse(wavelength=filt_lam, response=filt_tx,
-                   meta=dict(group_name='test', band_name='filt'))
-
-    # Each test galaxy is exactly one of the templates
-    coefficients = np.eye(3)
-
-    # Test absolute magnitudes
-    mt = test_templates.absolute_magnitudes(coefficients, 'test-filt')
-    m = mag_ab(lam, flam, 'test-filt')
-    np.testing.assert_allclose(mt, m)
-
-    # Test apparent magnitudes
-    redshift = np.array([1, 2, 3])
-    dm = Planck15.distmod(redshift).value
-    mt = test_templates.apparent_magnitudes(coefficients, redshift, 'test-filt', Planck15)
-    np.testing.assert_allclose(mt, m - 2.5*np.log10(1+redshift) + dm)
-
-    # Redshift interpolation test; linear interpolation sufficient over a small
-    # redshift range at low relative tolerance
-    z = np.linspace(0.1, 0.2, 3)
-    m_true = test_templates.apparent_magnitudes(coefficients, z, 'test-filt',
-                                                Planck15, resolution=4)
-    m_interp = test_templates.apparent_magnitudes(coefficients, z, 'test-filt',
-                                                  Planck15, resolution=2)
-    np.testing.assert_allclose(m_true, m_interp, rtol=1e-5)
-    assert not np.all(m_true == m_interp)
-
-
 @pytest.mark.skipif(not HAS_SPECLITE, reason='test requires speclite')
 def test_kcorrect_magnitudes():
 
diff --git a/skypy/utils/tests/test_decorators.py b/skypy/utils/tests/test_decorators.py
index bb96ba0..2f2832b 100644
--- a/skypy/utils/tests/test_decorators.py
+++ b/skypy/utils/tests/test_decorators.py
@@ -1,6 +1,6 @@
 import numpy as np
 import pytest
-from skypy.galaxy.spectrum import HAS_SPECLITE
+from skypy.utils.photometry import HAS_SPECLITE
 
 
 def test_broadcast_arguments():
diff --git a/skypy/utils/tests/test_photometry.py b/skypy/utils/tests/test_photometry.py
new file mode 100644
index 0000000..6f1d1ec
--- /dev/null
+++ b/skypy/utils/tests/test_photometry.py
@@ -0,0 +1,141 @@
+import numpy as np
+import pytest
+from skypy.utils.photometry import HAS_SPECLITE
+
+
+@pytest.mark.skipif(not HAS_SPECLITE, reason='test requires speclite')
+def test_mag_ab_standard_source():
+
+    from astropy import units
+    from speclite.filters import FilterResponse
+    from skypy.utils.photometry import mag_ab
+
+    # create a filter
+    filt_lam = np.logspace(0, 4, 1000)*units.AA
+    filt_tx = np.exp(-((filt_lam - 1000*units.AA)/(100*units.AA))**2)
+    filt_tx[[0, -1]] = 0
+    FilterResponse(wavelength=filt_lam, response=filt_tx,
+                   meta=dict(group_name='test', band_name='filt'))
+
+    # test that the AB standard source has zero magnitude
+    lam = filt_lam  # same grid to prevent interpolation issues
+    flam = 0.10885464149979998*units.Unit('erg s-1 cm-2 AA')/lam**2
+
+    m = mag_ab(lam, flam, 'test-filt')
+
+    assert np.isclose(m, 0)
+
+
+@pytest.mark.skipif(not HAS_SPECLITE, reason='test requires speclite')
+def test_mag_ab_redshift_dependence():
+
+    from astropy import units
+    from speclite.filters import FilterResponse
+    from skypy.utils.photometry import mag_ab
+
+    # make a wide tophat bandpass
+    filt_lam = [1.0e-10, 1.1e-10, 1.0e0, 0.9e10, 1.0e10]
+    filt_tx = [0., 1., 1., 1., 0.]
+    FilterResponse(wavelength=filt_lam, response=filt_tx,
+                   meta=dict(group_name='test', band_name='filt'))
+
+    # create a narrow gaussian source
+    lam = np.logspace(-11, 11, 1000)*units.AA
+    flam = np.exp(-((lam - 100*units.AA)/(10*units.AA))**2)*units.Unit('erg s-1 cm-2 AA-1')
+
+    # array of redshifts
+    z = np.linspace(0, 1, 11)
+
+    # compute the AB magnitude at different redshifts
+    m = mag_ab(lam, flam, 'test-filt', redshift=z)
+
+    # compare with expected redshift dependence
+    np.testing.assert_allclose(m, m[0] - 2.5*np.log10(1 + z))
+
+
+@pytest.mark.skipif(not HAS_SPECLITE, reason='test requires speclite')
+def test_mag_ab_multi():
+
+    from astropy import units
+    from skypy.utils.photometry import mag_ab
+    from speclite.filters import FilterResponse
+
+    # 5 redshifts
+    z = np.linspace(0, 1, 5)
+
+    # 2 Gaussian bandpasses
+    filt_lam = np.logspace(0, 4, 1000) * units.AA
+    filt_mean = np.array([[1000], [2000]]) * units.AA
+    filt_width = np.array([[100], [10]]) * units.AA
+    filt_tx = np.exp(-((filt_lam-filt_mean)/filt_width)**2)
+    filt_tx[:, [0, -1]] = 0
+    FilterResponse(wavelength=filt_lam, response=filt_tx[0],
+                   meta=dict(group_name='test', band_name='filt0'))
+    FilterResponse(wavelength=filt_lam, response=filt_tx[1],
+                   meta=dict(group_name='test', band_name='filt1'))
+
+    # 3 Flat Spectra
+    # to prevent issues with interpolation, collect all redshifted filt_lam
+    lam = []
+    for z_ in z:
+        lam = np.union1d(lam, filt_lam.value/(1+z_))
+    lam = lam*filt_lam.unit
+    A = np.array([[2], [3], [4]])
+    flam = A * 0.10885464149979998*units.Unit('erg s-1 cm-2 AA')/lam**2
+
+    # Compare calculated magnitudes with truth
+    magnitudes = mag_ab(lam, flam, ['test-filt0', 'test-filt1'], redshift=z)
+    truth = -2.5 * np.log10(A * (1+z)).T[:, :, np.newaxis]
+    assert magnitudes.shape == (5, 3, 2)
+    np.testing.assert_allclose(*np.broadcast_arrays(magnitudes, truth), rtol=1e-4)
+
+
+@pytest.mark.skipif(not HAS_SPECLITE, reason='test requires speclite')
+def test_template_spectra():
+
+    from astropy import units
+    from skypy.utils.photometry import mag_ab, SpectrumTemplates
+    from astropy.cosmology import Planck15
+    from speclite.filters import FilterResponse
+
+    class TestTemplates(SpectrumTemplates):
+        '''Three flat templates'''
+
+        def __init__(self):
+            self.wavelength = np.logspace(-1, 4, 1000)*units.AA
+            A = np.array([[2], [3], [4]]) * 0.10885464149979998
+            self.templates = A * units.Unit('erg s-1 cm-2 AA') / self.wavelength**2
+
+    test_templates = TestTemplates()
+    lam, flam = test_templates.wavelength, test_templates.templates
+
+    # Gaussian bandpass
+    filt_lam = np.logspace(0, 4, 1000)*units.AA
+    filt_tx = np.exp(-((filt_lam - 1000*units.AA)/(100*units.AA))**2)
+    filt_tx[[0, -1]] = 0
+    FilterResponse(wavelength=filt_lam, response=filt_tx,
+                   meta=dict(group_name='test', band_name='filt'))
+
+    # Each test galaxy is exactly one of the templates
+    coefficients = np.eye(3)
+
+    # Test absolute magnitudes
+    mt = test_templates.absolute_magnitudes(coefficients, 'test-filt')
+    m = mag_ab(lam, flam, 'test-filt')
+    np.testing.assert_allclose(mt, m)
+
+    # Test apparent magnitudes
+    redshift = np.array([1, 2, 3])
+    dm = Planck15.distmod(redshift).value
+    mt = test_templates.apparent_magnitudes(coefficients, redshift, 'test-filt', Planck15)
+    np.testing.assert_allclose(mt, m - 2.5*np.log10(1+redshift) + dm)
+
+    # Redshift interpolation test; linear interpolation sufficient over a small
+    # redshift range at low relative tolerance
+    z = np.linspace(0.1, 0.2, 3)
+    m_true = test_templates.apparent_magnitudes(coefficients, z, 'test-filt',
+                                                Planck15, resolution=4)
+    m_interp = test_templates.apparent_magnitudes(coefficients, z, 'test-filt',
+                                                  Planck15, resolution=2)
+    np.testing.assert_allclose(m_true, m_interp, rtol=1e-5)
+    assert not np.all(m_true == m_interp)