# -*- coding: utf-8 -*-
# Time-stamp: <2016-02-14 21:29 ycopin@lyonovae03.in2p3.fr>
"""
.. _indices:
indices - Spectral indices
==========================
.. Warning:: there's no verification that the proper computation is used
for a given feature. It's the responsibility of the user to check
coherence.
.. TODO:: take into account covariance between bandpasses.
"""
__author__ = "Yannick Copin <y.copin@ipnl.in2p3.fr>"
import warnings
import numpy as N
from . import synth
[docs]class IndexComputationWarning(RuntimeWarning):
"""
Spectral index computation warning.
"""
pass
[docs]def flambda_integ(lrange, spec, p=0):
r"""
Integral of weighted flambda over wavelength range.
.. math::
s_X = \int_{\lambda_{\min}}^{\lambda_{\max}}
\left(\frac{\lambda}{hc}\right)^{p} f_{\lambda} d\lambda
"""
# Define continuum bandpasses
xfilt = synth.Filter.bandpass("X", lrange[0], lrange[1])
# Compute continuum fluxes
x, dx = xfilt.integrate(spec, p=p)
return x, dx
[docs]def flambda_ratio(feature, spec, redshift=0):
r"""
flambda ratio continuum index.
Compute flambda ratio for continuum index `feature` on spectrum `spec` at
`redshift`:
.. math::
r = \frac{\langle f_{\lambda}\rangle_{\textrm{red}}}
{\langle f_{\lambda}\rangle_{\textrm{blue}}}
where :math:`\langle f_{\lambda}\rangle_X` is the average flux density (in
erg/s/cm²/Å) in the wavelength interval *X*:
.. math::
\langle f_{\lambda}\rangle_X =
\frac{1}{\lambda_{\max} - \lambda_{\min}} \int_X f_{\lambda} d\lambda
**Reference:** `Spinrad+ 1997
<http://cdsads.u-strasbg.fr/abs/1997ApJ...484..581S>`_
"""
# Continuum bands
bcont = N.array(feature['blue'], dtype=float)
rcont = N.array(feature['red'], dtype=float)
# Redshift correction (computation is done restframe)
if redshift:
bcont *= (1 + redshift)
rcont *= (1 + redshift)
# Compute continuum fluxes (flambda: p=0)
b, db = flambda_integ(bcont, spec, p=0)
r, dr = flambda_integ(rcont, spec, p=0)
# Compute ratio
if N.isnan(b):
warnings.warn(
"Cannot compute blue continuum flux "
"[{b[0]:.2f}, {b[1]:.2f}]".format(b=bcont),
IndexComputationWarning)
return N.nan, N.nan
elif N.isnan(r):
warnings.warn(
"Cannot compute red continuum flux "
"[{r[0]:.2f}, {r[1]:.2f}]".format(r=rcont),
IndexComputationWarning)
return N.nan, N.nan
elif b <= 0 or r <= 0:
warnings.warn(
f"Cannot compute continuum ratio {r} / {b}",
IndexComputationWarning)
return N.nan, N.nan
i = (r / (rcont[1] - rcont[0])) / (b / (bcont[1] - bcont[0]))
di = i * N.hypot(dr / r, db / b)
return i, di
[docs]def flambda_normed(feature, spec, redshift=0):
r"""
flambda normed ratio continuum index.
Compute flambda normed ratio for continuum index `feature` on spectrum
`spec` at `redshift`:
.. math::
r = \frac{2 \int_C f_{\lambda} d\lambda}
{\int_B f_{\lambda} d\lambda +
\int_R f_{\lambda} d\lambda}
where *C*, *B* and *R* are respectively the central, blue and red
bandpasses.
**Reference:** `Daddi+ 2005
<http://cdsads.u-strasbg.fr/abs/2005ApJ...626..680D>`_
"""
# Continuum bands
bcont = N.array(feature['blue'], dtype=float)
rcont = N.array(feature['red'], dtype=float)
ccont = N.array(feature['line'], dtype=float)
# Redshift correction (computation is done on obsframe spectra)
if redshift:
bcont *= (1 + redshift)
rcont *= (1 + redshift)
ccont *= (1 + redshift)
# Compute continuum fluxes (flambda: p=0)
b, db = flambda_integ(bcont, spec, p=0)
r, dr = flambda_integ(rcont, spec, p=0)
c, dc = flambda_integ(ccont, spec, p=0)
# Compute ratio
if N.isnan(b):
warnings.warn(
"Cannot compute blue continuum flux "
"[{b[0]:.2f}, {b[1]:.2f}]".format(b=bcont),
IndexComputationWarning)
return N.nan, N.nan
elif N.isnan(r):
warnings.warn(
"Cannot compute red continuum flux "
"[{r[0]:.2f}, {r[1]:.2f}]".format(r=rcont),
IndexComputationWarning)
return N.nan, N.nan
elif N.isnan(c):
warnings.warn(
"Cannot compute line flux "
"[{c[0]:.2f}, {c[1]:.2f}]".format(c=rcont),
IndexComputationWarning)
return N.nan, N.nan
elif b <= 0 or r <= 0 or c <= 0:
warnings.warn(
f"Cannot compute continuum ratio {c} / ({b} + {r})")
return N.nan, N.nan
bpr = b + r
i = 2 * c / bpr
di = i * N.sqrt((dc / c) ** 2 + (db ** 2 + dr ** 2) / bpr ** 2)
return i, di
[docs]def D4000_ratio(feature, spec, redshift=0):
r"""
D4000-like continuum index.
Compute D4000-like continuum index `feature` on
spectrum `spec` at `redshift`:
.. math::
r = \frac{\langle f_{\nu}\rangle_{\textrm{red}}}
{\langle f_{\nu}\rangle_{\textrm{blue}}}
where :math:`\langle f_{\nu}\rangle_X` is defined for historical
reasons (see `Gorgas+ 1999
<http://cdsads.u-strasbg.fr/abs/1999A%26AS..139...29G>`_) as:
.. math::
\langle f_{\nu}\rangle_X =
\frac{1}{\lambda_{\max} - \lambda_{\min}}
\int_X \lambda^2 f_{\lambda} d\lambda
**Reference:** `Spinrad+ 1997
<http://cdsads.u-strasbg.fr/abs/1997ApJ...484..581S>`_
"""
# Continuum bands
bcont = N.array(feature['blue'], dtype=float)
rcont = N.array(feature['red'], dtype=float)
# Redshift correction (computation is done restframe)
if redshift:
bcont *= (1 + redshift)
rcont *= (1 + redshift)
# Compute continuum fluxes (lambda**2 * flambda: p=2)
b, db = flambda_integ(bcont, spec, p=2)
r, dr = flambda_integ(rcont, spec, p=2)
# Compute ratio
if N.isnan(b):
warnings.warn(
"Cannot compute blue continuum flux "
"[{b[0]:.2f}, {b[1]:.2f}]".format(b=bcont),
IndexComputationWarning)
return N.nan, N.nan
elif N.isnan(r):
warnings.warn(
"Cannot compute red continuum flux "
"[{r[0]:.2f}, {r[1]:.2f}]".format(r=rcont),
IndexComputationWarning)
return N.nan, N.nan
elif b <= 0 or r <= 0:
warnings.warn(
f"Cannot compute continuum ratio {r} / {b}",
IndexComputationWarning)
return N.nan, N.nan
i = (r / (rcont[1] - rcont[0])) / (b / (bcont[1] - bcont[0]))
di = i * N.hypot(dr / r, db / b)
return i, di