"""Definitions for the `Line` class."""
import numexpr as ne
import numpy as np
from astropy import constants as c
from astropy import units as u
from mosfit.modules.seds.sed import SED
from mosfit.constants import BOL_BAND_INDEX, SQRT_2_PI
# Important: Only define one ``Module`` class per file.
[docs]
class Line(SED):
"""Line spectral energy distribution, modifies existing SED."""
C_CONST = c.c.cgs.value
[docs]
def process(self, **kwargs):
"""Process module."""
kwargs = self.prepare_input(self.key('luminosities'), **kwargs)
prt = self._printer
self._rest_t_explosion = kwargs[self.key('resttexplosion')]
self._times = kwargs[self.key('rest_times')]
self._seds = kwargs.get(self.key('seds'))
self._bands = kwargs['all_bands']
self._band_indices = kwargs['all_band_indices']
self._sample_wavelengths = kwargs['sample_wavelengths']
self._frequencies = kwargs['all_frequencies']
self._luminosities = kwargs[self.key('luminosities')]
self._line_wavelength = kwargs[self.key('line_wavelength')]
self._line_width = kwargs[self.key('line_width')]
self._line_time = kwargs[self.key('line_time')]
self._line_duration = kwargs[self.key('line_duration')]
self._line_amplitude = kwargs[self.key('line_amplitude')]
lw = self._line_wavelength # noqa: F841
ls = self._line_width
cc = self.C_CONST
# Some temp vars for speed.
zp1 = 1.0 + kwargs[self.key('redshift')]
czp1A = cc / (zp1 * u.Angstrom.cgs.scale)
amps = self._line_amplitude * np.array([
np.exp(-0.5 * (
(x - self._rest_t_explosion - self._line_time) /
self._line_duration) ** 2) for x in self._times])
if self._seds is None:
raise ValueError(prt.message('line_sed'))
seds = self.as_rectangular_seds(self._seds)
amps = np.asarray(amps, dtype=float)
seds = seds * (1.0 - amps)[:, None]
amps *= self._luminosities / (ls * SQRT_2_PI)
amps_dict = {}
evaled = False
for li, lum in enumerate(self._luminosities):
bi = self._band_indices[li]
if lum == 0.0:
continue
if bi == BOL_BAND_INDEX:
continue
bind = czp1A / self._frequencies[li] if bi < 0 else bi
if bind not in amps_dict:
# Leave `rest_wavs` in Angstroms.
if bi >= 0:
rest_wavs = self._sample_wavelengths[bi] / zp1
else:
rest_wavs = np.array([bind]) # noqa: F841
if not evaled:
amps_dict[bind] = ne.evaluate(
'exp(-0.5 * ((rest_wavs - lw) / ls) ** 2)')
evaled = True
else:
amps_dict[bind] = ne.re_evaluate()
profile = np.asarray(amps_dict[bind], dtype=float).ravel()
if profile.size == 1:
seds[li, 0] += amps[li] * profile[0]
else:
n = min(profile.size, seds.shape[1])
seds[li, :n] += amps[li] * profile[:n]
# seds[li][np.isnan(seds[li])] = 0.0
# Units of `seds` is ergs / s / Angstrom.
return {'sample_wavelengths': self._sample_wavelengths,
self.key('seds'): seds}