Source code for sbpy.spectroscopy.sources

# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""sbpy Spectroscopy Sources Module

Spectrophotometric classes that encasulate synphot.SpectralSource and
synphot.Observation in order to generate sbpy spectra and photometry.

Requires synphot.

"""

__all__ = [
    'BlackbodySource'
]

import numpy as np
from abc import ABC
import astropy.units as u
from astropy.utils.data import download_file, _is_url

try:
    import synphot
    from synphot import SpectralElement
except ImportError:
    synphot = None

    class SpectralElement:
        pass

from ..exceptions import SbpyException
from .. import bib

__doctest_requires__ = {
    'SpectralSource': ['synphot'],
    'BlackbodySource': ['synphot'],
}


class SinglePointSpectrumError(SbpyException):
    """Single point provided, but multiple values expected."""


class SynphotRequired(SbpyException):
    pass


class SpectralSource(ABC):
    """Abstract base class for SBPy spectral sources.

    Requires `~synphot`.


    Parameters
    ----------
    source : `~synphot.SourceSpectrum`
        The source spectrum.

    description : string, optional
        A brief description of the source spectrum.

    Attributes
    ----------
    wave        - Wavelengths of the source spectrum.
    fluxd       - Source spectrum.
    description - Brief description of the source spectrum.
    meta        - Meta data from ``source``, if any.

    """

    def __init__(self, source, description=None):
        if synphot is None:
            raise SynphotRequired(
                'synphot required for {}.'.format(self.__class__.__name__))

        self._source = source
        self._description = description

    @classmethod
    def from_array(cls, wave, fluxd, meta=None, **kwargs):
        """Create standard from arrays.


        Parameters
        ----------
        wave : `~astropy.units.Quantity`
            The spectral wavelengths.

        fluxd : `~astropy.units.Quantity`
            The spectral flux densities.

        meta : dict, optional
            Meta data.

        **kwargs
            Passed to object initialization.

        """
        if synphot is None:
            raise ImportError(
                'synphot required for {}.'.format(cls.__name__))

        source = synphot.SourceSpectrum(
            synphot.Empirical1D, points=wave, lookup_table=fluxd,
            meta=meta)

        return cls(source, **kwargs)

    @classmethod
    def from_file(cls, filename, wave_unit=None, flux_unit=None,
                  cache=True, **kwargs):
        """Load the source spectrum from a file.

        NaNs are dropped.


        Parameters
        ----------
        filename : string
            The name of the file.  See
            `~synphot.SourceSpectrum.from_file` for details.

        wave_unit, flux_unit : str or `~astropy.units.Unit`, optional
            Wavelength and flux units in the file.

        cache : bool, optional
            If ``True``, cache the contents of URLs.

        **kwargs
            Passed to object initialization.

        """

        if synphot is None:
            raise ImportError(
                'synphot required for {}.'.format(cls.__name__))

        if filename.lower().endswith(('.fits', '.fit', '.fz')):
            read_spec = synphot.specio.read_fits_spec
        else:
            read_spec = synphot.specio.read_ascii_spec

        # URL cache because synphot.SourceSpectrum.from_file does not
        if _is_url(filename):
            fn = download_file(filename, cache=True)
        else:
            fn = filename

        spec = read_spec(fn, wave_unit=wave_unit, flux_unit=flux_unit)
        i = np.isfinite(spec[1] * spec[2])
        source = synphot.SourceSpectrum(
            synphot.Empirical1D, points=spec[1][i], lookup_table=spec[2][i],
            meta={'header': spec[0]})

        return cls(source, **kwargs)

    @property
    def description(self):
        """Description of the source spectrum."""
        return self._description

    @property
    def wave(self):
        """Wavelengths of the source spectrum."""
        return self.source.waveset

    @property
    def fluxd(self):
        """The source spectrum."""
        return self.source(self._source.waveset, flux_unit='W / (m2 um)')

    @property
    def source(self):
        return self._source

    @property
    def meta(self):
        self._source.meta

    def __call__(self, wave_or_freq, unit=None):
        """Evaluate/interpolate the source spectrum.


        Parameters
        ----------
        wave_or_freq : `~astropy.units.Quantity`
            Requested wavelengths or frequencies of the resulting
            spectrum.

        unit : string, `~astropy.units.Unit`, optional
            Spectral units of the output (flux density).  If ``None``,
            the default depends on ``wave_or_freq``: W/(m2 μm) for
            wavelengths, Jy for frequencies.


        Returns
        -------
        fluxd : `~astropy.units.Quantity`
            The spectrum evaluated/interpolated to the requested
            wavelengths or frequencies.

        """

        from .. import units as sbu  # avoid circular dependency

        if unit is not None:
            unit = u.Unit(unit)
        elif wave_or_freq.unit.is_equivalent('m'):
            unit = u.Unit('W/(m2 um)')
        else:
            unit = u.Jy

        if unit.is_equivalent(sbu.VEGA):
            fluxd = self.source(wave_or_freq, 'W/(m2 um)').to(
                unit, sbu.spectral_density_vega(wave_or_freq))
        else:
            fluxd = self.source(wave_or_freq, unit)

        return fluxd

    def observe(self, wfb, unit=None, interpolate=False, **kwargs):
        """Observe source as through filters or spectrometer.

        Calls ``observe_bandpass``, ``observe_spectrum``, or
        ``self()``, as appropriate.


        Parameters
        ----------
        wfb : `~astropy.units.Quantity`, `~synphot.SpectralElement`
            Wavelengths, frequencies, or bandpasses.  May also be a
            list of ``SpectralElement``s.

        unit : string, `~astropy.units.Unit`, optional
            Spectral flux density units for the output.  If ``None``,
            the default depends on ``wfb``: W/(m2 μm) for wavelengths
            or bandpasses, Jy for frequencies.

        interpolate : bool, optional
            For wavelengths/frequencies, set to ``True`` for
            interpolation instead of rebinning.  Use this when the
            spectral resolution of the source is close to that of the
            requested wavelengths.

        **kwargs
            Additional keyword arguments for
            `~synphot.observation.Observation`, e.g., ``force``.


        Returns
        -------
        fluxd : `~astropy.units.Quantity`
            The spectrum rebinned.


        Raises
        ------
        SinglePointSpectrumError - If requested wavelengths or
            frequencies has only one value.


        Notes
        -----
        Method for spectra adapted from AstroBetter post by Jessica Lu:
        http://www.astrobetter.com/blog/2013/08/12/python-tip-re-sampling-spectra-with-pysynphot/

        """

        if isinstance(wfb, (list, tuple, SpectralElement)):
            lambda_eff, fluxd = self.observe_bandpass(
                wfb, unit=unit, **kwargs)
        elif isinstance(wfb, u.Quantity):
            if interpolate:
                fluxd = self(wfb, unit=unit)
            else:
                fluxd = self.observe_spectrum(wfb, unit=unit, **kwargs)
        else:
            raise TypeError('Unsupported type for `wfb` type: {}'
                            .format(type(wfb)))

        return fluxd

    def observe_bandpass(self, bp, unit=None, **kwargs):
        """Observe through a bandpass.


        Parameters
        ----------
        bp : `~synphot.SpectralElement`, list, or tuple
            Bandpass.

        unit : string, `~astropy.units.Unit`, optional
            Spectral flux density units for the output.  The default
            is W/(m2 μm).

        **kwargs
            Additional keyword arguments for
            `~synphot.observation.Observation`, e.g., ``force``.


        Returns
        -------
        lambda_eff : `~astropy.units.Quantity`
            Effective wavelength(s) of the observation(s).

        fluxd : `~astropy.units.Quantity`
            The spectrum rebinned.

        """

        from .. import units as sbu  # avoid circular dependency

        # promote single bandpasses to a list, but preserve number of
        # dimensions
        if isinstance(bp, (SpectralElement, str)):
            ndim = 0
            bp = [bp]
        else:
            ndim = np.ndim(bp)

        if unit is None:
            unit = u.Unit('W/(m2 um)')
        else:
            unit = u.Unit(unit)

        fluxd = np.ones(len(bp)) * unit
        for i in range(len(bp)):
            obs = synphot.Observation(self.source, bp[i], **kwargs)
            lambda_eff = obs.effective_wavelength()
            lambda_pivot = obs.bandpass.pivot()
            _fluxd = obs.effstim('W/(m2 um)')

            if unit.is_equivalent(sbu.VEGAmag):
                fluxd[i] = _fluxd.to(unit, sbu.spectral_density_vega(bp[i]))
            else:
                fluxd[i] = _fluxd.to(unit, u.spectral_density(lambda_pivot))

        if np.ndim(fluxd) != ndim:
            fluxd = fluxd.squeeze()

        return lambda_eff, fluxd

    def observe_spectrum(self, wave_or_freq, unit=None, **kwargs):
        """Observe source as through a spectrometer.

        .. Important:: This method works best when the requested
           spectral resolution is lower than the spectral resolution
           of the internal data.  If the requested
           wavelengths/frequencies are exactly the same as the
           internal spectrum, then the internal spectrum will be
           returned without binning.  This special case does not work
           for subsets of the wavelengths.


        Parameters
        ----------
        wave_or_freq : `~astropy.units.Quantity`
            Wavelengths or frequencies of the spectrum.  Spectral bins
            will be centered at these values.  The length must be
            larger than 1.

        unit : string, `~astropy.units.Unit`, optional
            Spectral flux density units for the output.  If ``None``,
            the default is W/(m2 μm) for wavelengths, Jy for
            frequencies.

        **kwargs
            Additional keyword arguments for
            `~synphot.observation.Observation`, e.g., ``force``.


        Returns
        -------
        fluxd : `~astropy.units.Quantity`
            The spectrum rebinned.


        Raises
        ------
        SinglePointSpectrumError - If requested wavelengths or
            frequencies has only one value.


        Notes
        -----
        Method for spectra adapted from AstroBetter post by Jessica Lu:
        http://www.astrobetter.com/blog/2013/08/12/python-tip-re-sampling-spectra-with-pysynphot/

        """

        from .. import units as sbu  # avoid circular dependency

        if np.size(wave_or_freq) == 1:
            raise SinglePointSpectrumError(
                'Multiple wavelengths or frequencies required for '
                'observe_spectrum.  Instead consider interpolation '
                'with {}().'
                .format(self.__class__.__name__))

        if unit is None:
            if wave_or_freq.unit.is_equivalent('m'):
                unit = u.Unit('W/(m2 um)')
            else:
                unit = u.Jy
        else:
            unit = u.Unit(unit)

        specele = synphot.SpectralElement(synphot.ConstFlux1D, amplitude=1)

        # Specele is defined over all wavelengths, but most spectral
        # standards are not.  force='taper' will affect retrieving
        # flux densities at the edges of the spectrum, but is
        # preferered to avoid wild extrapolation.
        kwargs['force'] = kwargs.get('force', 'taper')

        obs = synphot.Observation(
            self.source, specele, binset=wave_or_freq, **kwargs)

        if unit.is_equivalent(sbu.VEGAmag):
            fluxd = obs.sample_binned(flux_unit='W/(m2 um)').to(
                unit, sbu.spectral_density_vega(wave_or_freq))
        else:
            fluxd = obs.sample_binned(flux_unit=unit)

        return fluxd

    def color_index(self, wfb, unit):
        """Color index (magnitudes) and effective wavelengths.


        Parameters
        ----------
        wfb : `~astropy.units.Quantity` or tuple of `~synphot.SectralElement`
            Two wavelengths, frequencies, or bandpasses.

        unit : string or `~astropy.units.MagUnit`
            Units for the calculation, e.g., ``astropy.units.ABmag`` or
            ``sbpy.units.VEGAmag``.


        Returns
        -------
        eff_wave : `~astropy.units.Quantity`
            Effective wavelengths for each ``wfb``.

        ci : `~astropy.units.Quantity`
            Color index, ``m_0 - m_1``, where 0 and 1 are element
            indexes for ``wfb``.

        """

        eff_wave = []
        m = np.zeros(2) * u.Unit(unit)
        for i in range(2):
            if isinstance(wfb[i], u.Quantity):
                if wfb[i].unit.is_equivalent(u.Hz):
                    eff_wave.append(wfb[i].to(u.um, u.spectral()))
                else:
                    eff_wave.append(wfb[i])
                m[i] = self(eff_wave[i], unit=unit)
            elif isinstance(wfb[i], (list, tuple, SpectralElement)):
                w, m[i] = self.observe_bandpass(wfb[i], unit=unit)
                eff_wave.append(w)
            else:
                raise TypeError('Unsupported type for `wfb` type: {}'
                                .format(type(wfb[i])))

        ci = m[0] - m[1]

        return u.Quantity(eff_wave), ci


[docs]class BlackbodySource(SpectralSource): """Blackbody sphere. Spectral flux densities are calculated from ``pi * B(T)``, where ``B`` is the Planck function. Parameters ---------- T : `~astropy.units.Quantity`, required Temperature in Kelvin. """ def __init__(self, T=None): super().__init__(None, description='πB(T)') if T is None: raise TypeError('T is required.') self._T = u.Quantity(T, u.K) self._source = synphot.SourceSpectrum( synphot.BlackBody1D, temperature=self._T.value) * np.pi def __repr__(self): return '<BlackbodySource: T={}>'.format(self._T) @property def T(self): return self._T