# -*- coding: utf-8 -*-
# Time-stamp: <2021-03-16 16:50:18 ycopin>
"""
.. _lines:
lines - Emission/absorption line adjustment
===========================================
.. autosummary::
GaussianNormed1D
ComplexHa
DoubletOIII
"""
__author__ = "Yannick Copin <y.copin@ipnl.in2p3.fr>"
import numpy as N
import astropy.modeling as AM
SQRT2PI = 2.5066282746310002
[docs]class GaussianNormed1D(AM.Fittable1DModel):
r"""
*Flux-normalized* 1D-Gaussian.
.. math::
G_N(\lambda; f, \mu, \sigma) = \frac{f}{\sigma\sqrt{2\pi}}
\exp\left(-\frac{(\lambda - \mu)^2}{2\sigma^2}\right)
"""
# Due to Python issue #25731 (http://bugs.python.org/issue25731), one
# cannot mock a class with parameters:
# TypeError: object() takes no parameters
# flux = AM.Parameter(default=1)
# mean = AM.Parameter(default=0)
# stddev = AM.Parameter(default=1)
flux = AM.Parameter() #: Gaussian integrated flux
mean = AM.Parameter() #: Gaussian mean
stddev = AM.Parameter() #: Gaussian stddev
def bounding_box(self, factor=5.5):
x0 = self.mean.value
dx = factor * self.stddev
return (x0 - dx, x0 + dx)
[docs] @staticmethod
def evaluate(x, flux, mean, stddev):
return flux / (stddev * SQRT2PI) * \
N.exp(-0.5 * ((x - mean) / stddev)**2)
[docs] @staticmethod
def fit_deriv(x, flux, mean, stddev):
xnorm = (x - mean) / stddev
exp = N.exp(-0.5 * xnorm**2)
d_flux = exp / (stddev * SQRT2PI) # dG/dflux
d_mean = flux * d_flux / stddev # Temporary placeholder
d_stddev = d_mean * (xnorm**2 - 1) # dG/dstddev
d_mean *= xnorm # dG/dmean (final)
return [d_flux, d_mean, d_stddev]
[docs]class ComplexHa(AM.Fittable1DModel):
r"""
[NII] + Hα complex fittable model.
The [NII] + Hα complex is modeled as a tied sum of three
:class:`astropy.modeling.models.Gaussian1D`:
.. math::
P(\lambda) &= I_{H\alpha} \times G(\lambda; \lambda_{H\alpha}, \sigma) \\
& \quad + 3/4\, I_{NII} \times G(\lambda; \lambda_{NIIa}, \sigma) \\
& \quad + 1/4\, I_{NII} \times G(\lambda; \lambda_{NIIb}, \sigma) \\
G(\lambda; \mu, \sigma) &=
\exp\left(-\frac{(\lambda - \mu)^2}{2\sigma^2}\right)
where the 4 adjustable parameters are:
* :math:`\lambda_{H\alpha}` = `wHa`: Hα wavelength [Å]
* :math:`I_{H\alpha}` = `iHa`: Hα amplitude
* :math:`I_{NII}` = `iNII`: [NII] amplitude
* :math:`\sigma` = `sigma`: line width [Å]
and where [NII] doublet wavelengths :math:`\lambda_{NIIa}` and
:math:`\lambda_{NIIa}` are tied to the Hα wavelength
:math:`\lambda_{H\alpha}`.
.. Warning:: the amplitudes `iHa` and `iNII` are *not* the line fluxes,
since the standard :class:`astropy.modeling.models.Gaussian1D` is not
flux-normalized. Still, the amplitude ratio holds since the stddev is
common to all lines.
"""
Halpha = 6564.614 #: Hα reference wavelength [Å]
NIIa = 6585.27 #: NIIa reference wavelength [Å]
NIIb = 6549.86 #: NIIb reference wavelength [Å]
rNIIa = NIIa / Halpha # Wavelength wrt Halpha ref. wavelength
rNIIb = NIIb / Halpha # Wavelength wrt Halpha ref. wavelength
fNIIa = 0.75 #: Intensity ratio wrt NII amplitude
fNIIb = 0.25 #: Intensity ratio wrt NII amplitude
# List of parameters
# wHa = AM.Parameter(default=6564.614) #: Hα wavelength [Å]
# iHa = AM.Parameter(default=1) #: Hα amplitude
# iNII = AM.Parameter(default=0.33) #: [NII] amplitude
# sigma = AM.Parameter(default=1) #: Line width [Å]
wHa = AM.Parameter() #: Hα wavelength [Å]
iHa = AM.Parameter() #: Hα amplitude
iNII = AM.Parameter() #: [NII] amplitude
sigma = AM.Parameter() #: Line width [Å]
[docs] @staticmethod
def Ha_line(x, wHa, iHa, iNII, sigma):
"""
Hα line modeled as a single :class:`astropy.modeling.models.Gaussian1D`.
"""
return AM.models.Gaussian1D(amplitude=iHa, mean=wHa, stddev=sigma)(x)
[docs] @staticmethod
def NII_lines(x, wHa, iHa, iNII, sigma):
"""
[NII] doublet modeled as a tied sum of two
:class:`astropy.modeling.models.Gaussian1D`.
"""
return (AM.models.Gaussian1D(amplitude=iNII * ComplexHa.fNIIa,
mean=wHa * ComplexHa.rNIIa,
stddev=sigma)(x) +
AM.models.Gaussian1D(amplitude=iNII * ComplexHa.fNIIb,
mean=wHa * ComplexHa.rNIIb,
stddev=sigma)(x))
[docs] @staticmethod
def evaluate(x, wHa, iHa, iNII, sigma):
"""
Functional form of [NII]+Hα complex.
"""
return (ComplexHa.Ha_line(x, wHa, iHa, iNII, sigma) +
ComplexHa.NII_lines(x, wHa, iHa, iNII, sigma))
[docs] @staticmethod
def fit_deriv(x, wHa, iHa, iNII, sigma):
"""
Derivatives of [NII]+Hα complex.
"""
# Halpha derivatives
d_iHa, d_wHa, d_sHa = AM.models.Gaussian1D.fit_deriv(
x, amplitude=iHa, mean=wHa, stddev=sigma)
# NIIa derivatives
d_iNIIa, d_wNIIa, d_sNIIa = AM.models.Gaussian1D.fit_deriv(
x,
amplitude=iNII * ComplexHa.fNIIa,
mean=wHa * ComplexHa.rNIIa,
stddev=sigma)
# NIIb derivatives
d_iNIIb, d_wNIIb, d_sNIIb = AM.models.Gaussian1D.fit_deriv(
x,
amplitude=iNII * ComplexHa.fNIIb,
mean=wHa * ComplexHa.rNIIb,
stddev=sigma)
# Total derivatives
return [d_wHa + d_wNIIa * ComplexHa.rNIIa + d_wNIIb * ComplexHa.rNIIb,
d_iHa,
d_iNIIa * ComplexHa.fNIIa + d_iNIIb * ComplexHa.fNIIb,
d_sHa + d_sNIIa + d_sNIIb]
[docs]class DoubletOIII(AM.Fittable1DModel):
r"""
[OIII] doublet fittable model.
The [OIII] doublet is modeled as a tied sum of two
:class:`astropy.modeling.models.Gaussian1D`:
.. math::
P(\lambda) &= 3/4\, I_{OIII} \times G(\lambda; \lambda_{OIII}, \sigma) \\
& \quad + 1/4\, I_{OIII} \times G(\lambda; \lambda_{OIIIb}, \sigma) \\
G(\lambda; \mu, \sigma) &=
\exp\left(-\frac{(\lambda - \mu)^2}{2\sigma^2}\right)
where the 3 adjustable parameters are:
* :math:`\lambda_{OIII}` = `wOIII`: OIIIa wavelength [Å]
* :math:`I_{OIII}` = `iOIII`: [OIII] doublet amplitude
* :math:`\sigma` = `sigma`: line width [Å]
and where :math:`\lambda_{OIIIb}` wavelength is tied to the OIIIa
wavelength :math:`\lambda_{OIII}`.
.. Warning:: the amplitude `iOIII` is *not* the integrated line flux,
since the standard :class:`astropy.modeling.models.Gaussian1D` is not
flux-normalized. Still, the amplitude ratio holds since the stddev is
common to all lines.
"""
OIII = 5008.24 #: [OIIIa] reference wavelength [Å]
OIIIb = 4960.295 #: [OIIIb] reference wavelength [Å]
rOIIIb = OIIIb / OIII # Wavelength wrt OIIIa
fOIII = 0.75 #: Intensity ratio wrt OIII amplitude
fOIIIb = 0.25 #: Intensity ratio wrt OIII amplitude
# List of parameters
# wOIII = AM.Parameter(default=5008.24) #: OIIIa wavelength [Å]
# iOIII = AM.Parameter(default=1) #: [OIII] amplitude
# sigma = AM.Parameter(default=1) #: Line width [Å]
wOIII = AM.Parameter() #: OIIIa wavelength [Å]
iOIII = AM.Parameter() #: [OIII] amplitude
sigma = AM.Parameter() #: Line width [Å]
[docs] @staticmethod
def evaluate(x, wOIII, iOIII, sigma):
"""
Functional form of [OIII] doublet.
"""
return (AM.models.Gaussian1D(amplitude=iOIII * DoubletOIII.fOIII,
mean=wOIII,
stddev=sigma)(x) +
AM.models.Gaussian1D(amplitude=iOIII * DoubletOIII.fOIIIb,
mean=wOIII * DoubletOIII.rOIIIb,
stddev=sigma)(x))
[docs] @staticmethod
def fit_deriv(x, wOIII, iOIII, sigma):
"""
Derivatives of [OIII] doublet.
"""
# OIIIa derivatives
d_iOIII, d_wOIII, d_sOIII = AM.models.Gaussian1D.fit_deriv(
x,
amplitude=iOIII * DoubletOIII.fOIII,
mean=wOIII,
stddev=sigma)
# OIIIb derivatives
d_iOIIIb, d_wOIIIb, d_sOIIIb = AM.models.Gaussian1D.fit_deriv(
x,
amplitude=iOIII * DoubletOIII.fOIIIb,
mean=wOIII * DoubletOIII.rOIIIb,
stddev=sigma)
# Total derivatives
return [d_wOIII + d_wOIIIb * DoubletOIII.rOIIIb,
d_iOIII * DoubletOIII.fOIII + d_iOIIIb * DoubletOIII.fOIIIb,
d_sOIII + d_sOIIIb]
[docs]class FitterWithPrior(AM.fitting.LevMarLSQFitter):
"""
LevMarLSQFitter with additional bayesian priors on parameters.
"""
def __init__(self, model, priors=None):
super().__init__()
assert isinstance(model, AM.FittableModel)
self.model = model
#: List of prior functions (or None)
self.priors = [ None for p in self.model.parameters ]
if priors:
self.set_priors(priors) # {param_name: prior_fn}
[docs] def objective_function(self, fps, *args):
"""
Standard objective function, with additional bayesian prior terms.
"""
objf = super().objective_function(fps, *args) # Standard objective function
hterms = self.get_hterms() # Hyper-terms
objf += sum(hterms) # Add hyper-terms to objf
return objf
[docs] def set_priors(self, priors):
"""
Set priors {param_name: prior_fn}.
"""
for param_name in priors: # Basic check on parameter names
try:
i = self.model.param_names.index(param_name)
except ValueError:
raise KeyError(
f"Unknown model parameter '{param_name}'.")
self.priors[i] = priors[param_name] # Update prior list
[docs] def get_hterms(self, param_names=[]):
"""
Compute prior hypter-terms.
"""
if not param_names:
param_names = self.model.param_names
else:
assert all( param_name in self.model.param_names
for param_name in param_names ), \
"Some parameter names are unknown."
ipars = [ self.model.param_names.index(param_name)
for param_name in param_names ]
# Restrict to parameters with priors
ipars = [ ipar for ipar in ipars if self.priors[ipar] is not None ]
hterms = [ self.priors[ipar](self.model.parameters[ipar])
for ipar in ipars ]
return hterms # [ hyper-terms ]