Source code for mosfit.modules.seds.blackbody

"""Definitions for the `Blackbody` class."""
from math import pi

import numexpr as ne
import numpy as np
from astropy import constants as c
from astropy import units as u
from mosfit.constants import FOUR_PI
from mosfit.modules.seds.sed import SED


# Important: Only define one ``Module`` class per file.


[docs]class Blackbody(SED): """Blackbody spectral energy dist. for given temperature and radius.""" C_CONST = c.c.cgs.value FLUX_CONST = FOUR_PI * ( 2.0 * c.h * c.c ** 2 * pi).cgs.value * u.Angstrom.cgs.scale X_CONST = (c.h * c.c / c.k_B).cgs.value STEF_CONST = (4.0 * pi * c.sigma_sb).cgs.value
[docs] def process(self, **kwargs): """Process module.""" lum_key = self.key('luminosities') kwargs = self.prepare_input(lum_key, **kwargs) self._luminosities = kwargs[lum_key] self._bands = kwargs['all_bands'] self._band_indices = kwargs['all_band_indices'] self._frequencies = kwargs['all_frequencies'] self._radius_phot = kwargs[self.key('radiusphot')] self._temperature_phot = kwargs[self.key('temperaturephot')] xc = self.X_CONST # noqa: F841 fc = self.FLUX_CONST # noqa: F841 cc = self.C_CONST temperature_phot = self._temperature_phot # Some temp vars for speed. zp1 = 1.0 + kwargs[self.key('redshift')] Azp1 = u.Angstrom.cgs.scale / zp1 czp1 = cc / zp1 seds = [] rest_wavs_dict = {} evaled = False for li, lum in enumerate(self._luminosities): bi = self._band_indices[li] if lum == 0.0: seds.append(np.zeros(len( self._sample_wavelengths[bi]) if bi >= 0 else 1)) continue if bi >= 0: rest_wavs = rest_wavs_dict.setdefault( bi, self._sample_wavelengths[bi] * Azp1) else: rest_wavs = np.array( # noqa: F841 [czp1 / self._frequencies[li]]) radius_phot = self._radius_phot[li] # noqa: F841 temperature_phot = self._temperature_phot[li] # noqa: F841 if not evaled: seds.append(ne.evaluate( 'fc * radius_phot**2 / rest_wavs**5 / ' 'expm1(xc / rest_wavs / temperature_phot)')) evaled = True else: try: seds.append(ne.re_evaluate()) except: seds.append(fc * radius_phot**2. / rest_wavs**5. / np.expm1(xc / rest_wavs / temperature_phot)) seds[-1][np.isnan(seds[-1])] = 0.0 seds = self.add_to_existing_seds(seds, **kwargs) # Units of `seds` is ergs / s / Angstrom. return {'sample_wavelengths': self._sample_wavelengths, 'seds': seds}