# functions.py: Various spectral functions and base function class
#
# Authors: William Cleveland (USRA),
# Adam Goldstein (USRA) and
# Daniel Kocevski (NASA)
#
# Portions of the code are Copyright 2020 William Cleveland and
# Adam Goldstein, Universities Space Research Association
# All rights reserved.
#
# Written for the Fermi Gamma-ray Burst Monitor (Fermi-GBM)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <https://www.gnu.org/licenses/>.
#
import numpy as np
from scipy.integrate import trapz, quad
from scipy.special import erfc
[docs]class Function:
"""Abstract class for a 1-D function (e.g. photon model).
Sub-classes should be designed with the following requirements:
1. Attribute :attr:`nparams` containing the total number of parameters of the
function (even parameters that are fixed to a particular value during fitting)
2. Method :meth:`eval` that accepts two arguments: an array of parameter values
of length :attr:`nparams`; and the array of abscissa values.
Sub-classes can optionally add the following for complete and enhanced
functionality:
1. Attribute :attr:`param_list` that is a list of parameter names in the order
they are accepted in the :meth:`eval` method. Each entry in :attr:`param_list` is
a tuple of the form (<param_name>, <units>, <description>). If no units
or description are needed, those fields should be a null string ''. If
omitted, these will be filled to default values ('param0', '', ''), etc.
2. Attribute :attr:`default_values` that is a list of default parameter values.
This is primarily used for parameter estimation. Default is 1.0 for each
parameter
3. Attribute :attr:`delta_abs` that is a list of max absolute fitting step for
each parameter. This is used for parameter estimation. Default is 0.1
for each parameter.
4. Attribute :attr:`delta_rel` that is a list of max relative fitting step for
each parameter. This is used for paramter estimation. Default is 0.01
for each parameter.
5. Attribute :attr:`min_values` that is a list of minimum values each parameter
can take. This is used primarily for parameter estimation. Default is
-np.inf for each parameter.
6) Attribute :attr:`max_values` that is a list of maximum values each parameter
can take. This is used primarily for parameter estimation. Default is
np.inf for each parameter.
7) Attribute :attr:`free` that is a boolean list indicating if each parameter
is allowed to be free for fitting. False indicates that the parameter
will be fixed at its default value. This is used for parameter estimation.
Default is True for each parameter.
Attributes:
default_values (list): The default values for each parameter.
If not defined, each defaults to 1.0
delta_abs (list): The maximum absolute fitting step size for each
parameter. If not defined, each default to 0.1
delta_rel (list): The maximum relative fitting step size for each
parameter. If not defined, each default to 0.01
free (list): For each parameter, mark True if left free for fitting,
otherwise mark False to fix at the default value.
If not defined, each default to True.
max_values (list): The maximum possible value for each parameter.
If not defined, each default to np.inf
min_values (list): The minimum possible value for each parameter.
If not defined, each default to -np.inf
name (str): The name of the function
nparams (int): Number of parameters in the function.
num_components (int): Number of function components. For a normal
function this is one, for a :class:`SuperFunction`,
this can be > 1.
param_list (list of 3-tuples):
A list of tuples corresponding to the metadata for the
parameters: (<param_name>, <units>, <description>)
"""
__slots__ = ['nparams', 'name', 'param_list', 'default_values',
'delta_abs',
'delta_rel', 'min_values', 'max_values', 'free']
def __init__(self):
self._num_components = 1
if not hasattr(self, 'nparams'):
raise AttributeError('Function MUST have nparams defined')
if not hasattr(self, 'param_list'):
self.param_list = [('param' + str(i), '', '') for i in
range(self.nparams)]
if not hasattr(self, 'name'):
self.name = self.__class__.__name__
if not hasattr(self, 'default_values'):
self.default_values = [1.0] * self.nparams
if not hasattr(self, 'delta_abs'):
self.delta_abs = [0.1] * self.nparams
if not hasattr(self, 'delta_rel'):
self.delta_rel = [0.01] * self.nparams
if not hasattr(self, 'min_values'):
self.min_values = [-np.inf] * self.nparams
if not hasattr(self, 'max_values'):
self.max_values = [np.inf] * self.nparams
if not hasattr(self, 'free'):
self.free = [True] * self.nparams
@property
def num_components(self):
return self._num_components
[docs] def fit_eval(self, params, x, **kwargs):
"""Evaluate the function for a fit. This differs from the :meth:`eval`
because the ``params`` vector should only contain the free fit parameters.
The parameters that are set to fixed will automatically be passed to
:meth:`eval` at their fixed value. This enables a fit to be performed
withfixed parameters and the associated jacobian and covariance matrices
will have the appropriate shape.
Args:
params (iterable): The parameter vector of the free parameters
x (np.array): The values at which to evaluate the function
Returns:
np.array: The function evaluated at the ``x`` points
"""
full_params = np.zeros(self.nparams)
free = np.array(self.free)
full_params[free] = params
full_params[~free] = np.array(self.default_values)[~free]
return self.eval(full_params, x, **kwargs)
[docs] def eval(self, params, x):
"""Evaluate the function
Args:
params (iterable): The parameter vector of the free parameters
x (np.array): The values at which to evaluate the function
Returns:
np.array: The function evaluated at the ``x`` points
"""
raise NotImplementedError('eval() has not been defined')
[docs] def integrate(self, params, erange, num_points=1000, log=True,
energy=False):
"""Integrate the function over energy to produce a flux.
Args:
params (iterable): The parameter vector, which can be of length
:attr:`nparams` or equal to the length of the
non-fixed (free) parameters.
erange (float, float): The energy range over which to integrate
num_points (int, optional): The number of points used for
integration. Default is 1000.
log (bool, optional): If True, the integration grid is made in
log-space, False it is linear. Default is True.
energy (bool, optional): If True, perform integration of E*F(E) to
produce an energy flux, otherwise
integration is over F(E) to produce a
photon flux. Default is False.
Returns:
float: The photon or energy flux
"""
# If given all params, use eval
if len(params) == self.nparams:
func = self.eval
# otherwise, use fit_eval (auto apply fixed params)
elif len(params) == np.sum(self.free):
func = self.fit_eval
else:
raise ValueError('Incorrect number of parameters')
# evenly-spaced energies in log space
if log:
energies = np.logspace(*np.log10(erange), num_points)
# evenly-spaced energies in linear space
else:
energies = np.linspace(*erange, num_points)
spectrum = func(params, energies)
# photon flux
if not energy:
flux = trapz(spectrum, energies)
# energy flux
else:
flux = trapz(energies * spectrum, energies) * 1.602e-9
return flux
[docs] def parameter_bounds(self, apply_state=True):
"""The valid parameter bounds
Args:
apply_state (bool, optional):
If True, will only return the bounds for free parameters.
Default is True.
Returns:
[(float, float), ...]: The list of bounds [(min, max), ...]
"""
bounds = list(zip(*(self.min_values, self.max_values)))
if apply_state:
bounds = [bounds[i] for i in range(self.nparams) if self.free[i]]
return bounds
[docs] @classmethod
def add(cls, func1, func2):
"""Add two Functions. Can also use the ``+`` operator.
Args:
func1 (:class:`Function`): A function
func2 (:class:`Function`): A function to add to ``func1``.
Returns:
:class:`SuperFunction`: A function containing the sum of the two functions
"""
obj = cls._super_function(func1, func2, np.sum)
obj.names = [func1.name, func2.name]
obj.name = '+'.join(obj.names)
return obj
[docs] @classmethod
def multiply(cls, func1, func2):
"""Multiply two Functions. Can also use the ``*`` operator.
Args:
func1 (:class:`Function`): A function
func2 (:class:`Function`): A function to multiply to ``func1``.
Returns:
:class:`SuperFunction`: A function containing the product of the \
two functions
"""
obj = cls._super_function(func1, func2, np.prod)
obj.names = [func1.name, func2.name]
obj.name = '*'.join(obj.names)
return obj
@classmethod
def _super_function(cls, func1, func2, operator):
"""Creates a SuperFunction from two functions (one may itself be
a SuperFunction) and an operator.
Args:
func1 (:class:`Function` or :class:`SuperFunction`): A function
func2 (:class:`Function` or :class:`SuperFunction`): Another function.
operator (<function>): The operator to be performed on the two functions.
Returns:
:class:`SuperFunction`: A function containing the product of the \
two functions
"""
obj = SuperFunction()
obj._num_components = func1._num_components + func2._num_components
obj.nparams = func1.nparams + func2.nparams
# func1 could be a SuperFunction or just a simple Function
if isinstance(func1, SuperFunction):
obj._operator = func1._operator
obj._sub_names = func1._sub_names
obj._sub_param_list = func1._sub_param_list
obj._sub_nparams = func1._sub_nparams
obj._sub_functions = func1._sub_functions
obj._operator.append(operator)
obj._sub_names.append(func2.name)
obj._sub_param_list.append(func2.param_list)
obj._sub_nparams.append(func2.nparams)
obj._sub_functions.append(func2.eval)
elif isinstance(func2, SuperFunction):
obj._operator = [operator]
obj._sub_names = [func1.name]
obj._sub_param_list = [func1.param_list]
obj._sub_nparams = [func1.nparams]
obj._sub_functions = [func1.eval]
obj._operator.extend(func2._operator)
obj._sub_names.extend(func2._sub_names)
obj._sub_param_list.extend(func2._sub_param_list)
obj._sub_nparams.extend(func2._sub_nparams)
obj._sub_functions.extend(func2._sub_functions)
else:
obj._operator = [operator]
obj._sub_names = [func1.name]
obj._sub_param_list = [func1.param_list]
obj._sub_nparams = [func1.nparams]
obj._sub_functions = [func1.eval]
obj._sub_names.append(func2.name)
obj._sub_param_list.append(func2.param_list)
obj._sub_nparams.append(func2.nparams)
obj._sub_functions.append(func2.eval)
# modify param lists to have component names
obj.param_list = obj._modify_param_list()
obj.default_values = func1.default_values + func2.default_values
obj.delta_abs = func1.delta_abs + func2.delta_abs
obj.delta_rel = func1.delta_rel + func2.delta_rel
obj.min_values = func1.min_values + func2.min_values
obj.max_values = func1.max_values + func2.max_values
obj.free = func1.free + func2.free
return obj
def __add__(self, other):
return Function.add(self, other)
def __mul__(self, other):
return Function.multiply(self, other)
[docs]class SuperFunction(Function):
"""Abstract class for a sum or product of functions.
Attributes:
default_values (list): The default values for each parameter.
If not defined, each defaults to 1.0
delta_abs (list): The maximum absolute fitting step size for each
parameter. If not defined, each default to 0.1
delta_rel (list): The maximum relative fitting step size for each
parameter. If not defined, each default to 0.01
free (list): For each parameter, mark True if left free for fitting,
otherwise mark False to fix at the default value.
If not defined, each default to True.
max_values (list): The maximum possible value for each parameter.
If not defined, each default to np.inf
min_values (list): The minimum possible value for each parameter.
If not defined, each default to -np.inf
name (str): The name of the function
nparams (int): Number of parameters in the function.
num_components (int): Number of function components. For a normal
function this is one, for a :class:`SuperFunction`,
this can be > 1.
param_list (list of 3-tuples):
A list of tuples corresponding to the metadata for the
parameters: (<param_name>, <units>, <description>)
"""
nparams = 0
def __init__(self):
super().__init__()
self._sub_names = []
self._sub_nparams = []
self._sub_param_list = []
self._sub_functions = []
self._operator = []
[docs] def eval(self, params, x, components=False):
"""Evaluate the function
Args:
params (iterable): The parameter vector of the free parameters
x (np.array): The values at which to evaluate the function
components (bool, optional): If True, evaluate for each component.
Default is False.
Returns:
np.array or list of np.array:
The function evaluated at ``x`` points or the function
evaluated at ``x`` points for each component.
"""
ind = np.concatenate(([0], np.cumsum(self._sub_nparams)))
evals = [self._sub_functions[i](params[ind[i]:ind[i + 1]], x) \
for i in range(len(self._sub_functions))]
if components:
return evals
else:
the_eval = evals[0]
for i in range(self.num_components - 1):
the_eval = self._operator[i]((the_eval, evals[i + 1]), axis=0)
return the_eval
def _modify_param_list(self):
"""Modify the parameter list so the parmeter names contain the parent
function name.
Returns:
[(str, str, str), ...]: The modified parameter list
"""
param_list = []
for name, plist in zip(self._sub_names, self._sub_param_list):
new_list = [(name + ': ' + p[0], p[1], p[2]) for p in plist]
param_list.extend(new_list)
return param_list
# ---- functions
[docs]class PowerLaw(Function):
r"""A power law function:
:math:`F(E) = A \bigl(\frac{E}{E_{piv}}\bigr)^{\lambda}`,
where
* :math:`A` is the amplitude in ph/s/cm^2/keV
* :math:`\lambda` is the power-law index
* :math:`E_{piv}` is the pivot energy in keV. Nominally fixed to 100 keV.
"""
# The total number of parameters
nparams = 3
# The list of parameters: (parameter name, units, description)
param_list = [('A', 'ph/s/cm^2/keV', 'Amplitude'),
('index', '', 'Photon index'),
('Epiv', '', 'Pivot energy')]
# -----------------------------------------------------------------------
# The following attributes are used for fitting and parameter estimation
# The default initialization values for each parameter
default_values = [0.1, -2.0, 100.0]
# The maximum absolute fitting step for each parameter
delta_abs = [0.1, 0.1, 0.1]
# The maximum relative fitting step for each parameter
delta_rel = [0.01, 0.01, 0.01]
# The minimum value that each parameter can take
min_values = [1e-10, -20.0, 0.01]
# The maximum value that each parameter can take
max_values = [np.inf, np.inf, np.inf]
# If True, then the parameter is a free parameter, otherwise it is fixed
free = [True, True, False]
def eval(self, params, x):
return params[0] * (x / params[2]) ** params[1]
[docs]class Comptonized(Function):
r"""An exponentially cut-off power law, parametrized by the
peak of the :math:`\nu F_{\nu}` spectrum:
:math:`F(E) = A e^{-(2+\lambda)E/E_{peak}} \bigl(\frac{E}{E_{piv}}\bigr)^\lambda`
where
* :math:`A` is the amplitude in ph/s/cm^2/keV
* :math:`E_{peak}` is the :math:`\nu F_{\nu}` peak in keV
* :math:`\lambda` is the power-law index
* :math:`E_{piv}` is the pivot energy in keV. Nominally fixed to 100 keV.
"""
nparams = 4
param_list = [('A', 'ph/s/cm^2/keV', 'Amplitude'),
('Epeak', 'keV', 'SED Peak'),
('index', '', 'Photon index'),
('Epiv', 'keV', 'Pivot energy')]
default_values = [0.01, 300.0, -0.5, 100.0]
delta_abs = [0.1, 0.1, 0.1, 0.1]
delta_rel = [0.01, 0.01, 0.01, 0.01]
min_values = [1e-10, 0.01, -1.9, 0.01]
max_values = [np.inf, np.inf, 20.0, np.inf]
free = [True, True, True, False]
def eval(self, params, x):
A, Epeak, index, Epiv = params
if index == -2.0:
return np.zeros_like(x)
logfxn = -x * (2.0 + index) / Epeak + index * np.log(x / Epiv)
return A * np.exp(logfxn)
[docs]class Band(Function):
r"""The Band GRB function (a type of smoothly broken power law), parametrized
by the peak of the :math:`\nu F_{\nu}` spectrum:
:math:`F(E) = A
\begin{cases}
\bigl(\frac{E}{E_{piv}}\bigr)^\alpha e^{-(2+\alpha)E/E_{peak}},
\text{ if } {E < \frac{(\alpha-\beta)E_{peak}}{2+\alpha}}\\
\bigl(\frac{(\alpha-\beta)E_{peak}}{(2+\alpha)E_{piv}}\bigr)^{\alpha-\beta}
\exp{\beta-\alpha}\bigl(\frac{E}{E_{piv}}\bigr)^\beta, \text{ otherwise }
\end{cases}`
where
* :math:`A` is the amplitude in ph/s/cm^2/keV
* :math:`E_{peak}` is the :math:`\nu F_{\nu}` peak in keV
* :math:`\alpha` is the low-energy power-law index
* :math:`\beta` is the high-energy power-law index
* :math:`E_{piv}` is the pivot energy in keV. Nominally fixed to 100 keV.
"""
nparams = 5
param_list = [('A', 'ph/s/cm^2/keV', 'Amplitude'),
('Epeak', 'keV', 'SED Peak'),
('alpha', '', 'Low-Energy Photon index'),
('beta', '', 'High-Energy Photon index'),
('Epiv', 'keV', 'Pivot energy')]
default_values = [0.01, 500.0, -0.5, -2.5, 100.0]
delta_abs = [0.1, 0.1, 0.1, 0.1, 0.1]
delta_rel = [0.01, 0.01, 0.01, 0.01, 0.01]
min_values = [1e-10, 0.01, -1.9, -10.0, 0.01]
max_values = [np.inf, np.inf, 20.0, -2.0001, np.inf]
free = [True, True, True, True, False]
def eval(self, params, x):
A, Epeak, alpha, beta, Epiv = params
e0 = Epeak / (2.0 + alpha)
ebreak = (alpha - beta) * e0
idx = (x < ebreak)
logfxn = np.zeros(len(x), dtype=float)
logfxn[idx] = np.log(A) + alpha * np.log(x[idx] / Epiv) - x[idx] / e0
dindex = alpha - beta
logfxn[~idx] = np.log(A) + dindex * np.log(dindex * e0 / Epiv) - \
dindex + beta * np.log(x[~idx] / Epiv)
return np.exp(logfxn)
[docs]class BandOld(Function):
r"""The Band GRB function (a type of smoothly broken power law), using the
original parameterization:
:math:`F(E) = A
\begin{cases}
\bigl(\frac{E}{E_{piv}}\bigr)^\alpha e^{-\frac{E}{E_0}},
\text{ if } {E < (\alpha-\beta) E_0}\\
\bigl(\frac{(\alpha-\beta) E_0}{E_{piv}}\bigr)^{(\alpha-\beta)}
e^{(\beta-\alpha)} \bigl(\frac{E}{E_{piv}}\bigr)^\beta, \text{ otherwise }
\end{cases}`
where
* :math:`A` is the amplitude in ph/s/cm^2/keV
* :math:`E_0` is the break energy in keV
* :math:`\alpha` is the low-energy power-law index
* :math:`\beta` is the high-energy power-law index
* :math:`E_{piv}` is the pivot energy in keV. Nominally fixed to 100 keV.
References:
`Band, D. et al. 1993, ApJ, 413, 281
<http://adsabs.harvard.edu/full/1993ApJ...413..281B>`_
"""
nparams = 5
param_list = [('A', 'ph/s/cm^2/keV', 'Amplitude'),
('E0', 'keV', 'Break Energy'),
('alpha', '', 'Low-Energy Photon index'),
('beta', '', 'High-Energy Photon index'),
('Epiv', 'keV', 'Pivot energy')]
default_values = [0.01, 200.0, -0.5, -2.5, 100.0]
delta_abs = [0.1, 0.1, 0.1, 0.1, 0.1]
delta_rel = [0.01, 0.01, 0.01, 0.01, 0.01]
min_values = [1e-10, 0.01, -1.9, -10.0, 0.01]
max_values = [np.inf, np.inf, 20.0, -2.0001, np.inf]
free = [True, True, True, True, False]
def eval(self, params, x):
A, e0, alpha, beta, Epiv = params
ebreak = (alpha - beta) * e0
idx = (x < ebreak)
logfxn = np.zeros(len(x), dtype=float)
logfxn[idx] = np.log(A) + alpha * np.log(x[idx] / Epiv) - x[idx] / e0
dindex = alpha - beta
logfxn[~idx] = np.log(A) + dindex * np.log(dindex * e0 / Epiv) - \
dindex + beta * np.log(x[~idx] / Epiv)
return np.exp(logfxn)
[docs]class BrokenPowerLaw(Function):
r"""A sharply-broken power law:
:math:`F(E) = A
\begin{cases}
\bigl(\frac{E}{E_{piv}}\bigr)^\alpha, & \text{ if } {E \leq E_b}\\
\bigl(\frac{E_b}{E_{piv}}\bigr)^\alpha \bigl(\frac{E}{E_b}\bigr)^\beta, &
\text{ otherwise }
\end{cases}`
where
* :math:`A` is the amplitude in ph/s/cm^2/keV
* :math:`E_{b}` is the break energy in keV
* :math:`\alpha` is the low-energy power-law index
* :math:`\beta` is the high-energy power-law index
* :math:`E_{piv}` is the pivot energy in keV. Nominally fixed to 100 keV.
"""
nparams = 5
param_list = [('A', 'ph/s/cm^2/keV', 'Amplitude'),
('Eb', 'keV', 'Break Energy'),
('alpha', '', 'Low-Energy Photon index'),
('beta', '', 'High-Energy Photon index'),
('Epiv', 'keV', 'Pivot energy')]
default_values = [0.01, 700.0, -1.0, -2.0, 100.0]
delta_abs = [0.1, 0.1, 0.1, 0.1, 0.1]
delta_rel = [0.01, 0.01, 0.01, 0.01, 0.01]
min_values = [1e-10, 0.01, -1.9, -10.0, 0.01]
max_values = [np.inf, np.inf, 20.0, 20.0, np.inf]
free = [True, True, True, True, False]
def eval(self, params, x):
A, ebreak, alpha, beta, Epiv = params
mask = (x <= ebreak)
fxn = np.zeros(len(x), dtype=float)
fxn[mask] = (x[mask] / Epiv) ** alpha
fxn[~mask] = (ebreak / Epiv) ** alpha * (x[~mask] / ebreak) ** beta
return A * fxn
[docs]class DoubleBrokenPowerLaw(Function):
r"""A sharply-broken power law with two breaks:
:math:`F(E) = A
\begin{cases}
\bigl(\frac{E}{E_{piv}}\bigr)^\alpha,
\text{ if } {E \leq E_{b1}}\\
\bigl(\frac{E_{b1}}{E_{piv}}\bigr)^\alpha \bigl(\frac{E}{E_b}\bigr)^\beta,
\text{ if } {E > E_{b1} \text{ and } E \leq E_{b2}}\\
\bigl(\frac{E_{b1}}{E_{piv}}\bigr)^\alpha \bigl(\frac{E_{b1}}{E_{b2}}\bigr)^\beta
\bigl(\frac{E}{E_{b2}}\bigr)^\gamma,
\text{ otherwise}
\end{cases}`
where
* :math:`A` is the amplitude in ph/s/cm^2/keV
* :math:`E_{b1}` is the first break energy in keV
* :math:`E_{b2}` is the second break energy in keV
* :math:`\alpha` is the low-energy power-law index
* :math:`\beta` is the mid-energy power-law index
* :math:`\gamma` is the high-energy power-law index
* :math:`E_{piv}` is the pivot energy in keV. Nominally fixed to 100 keV.
"""
nparams = 7
param_list = [('A', 'ph/s/cm^2/keV', 'Amplitude'),
('Eb1', 'keV', 'Lower Break Energy'),
('Eb2', 'keV', 'Upper Break Energy'),
('alpha', '', 'Low-Energy Photon index'),
('beta', '', 'Mid-Energy Photon index'),
('gamma', '', 'High-Energy Photon index'),
('Epiv', 'keV', 'Pivot energy')]
default_values = [0.01, 200.0, 500, -0.5, -1.0, -2.0, 100.0]
delta_abs = [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]
delta_rel = [0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01]
min_values = [1e-10, 0.01, 0.01, -1.9, -10.0, -20.0, 0.01]
max_values = [np.inf, np.inf, np.inf, 20.0, 20.0, 20.0, np.inf]
free = [True, True, True, True, True, True, False]
def eval(self, params, x):
A, ebreak1, ebreak2, alpha, beta, gamma, Epiv = params
fxn = np.zeros(len(x), dtype=float)
mask1 = (x <= ebreak1)
fxn[mask1] = (x[mask1] / Epiv) ** alpha
mask2 = (x > ebreak1) & (x <= ebreak2)
fxn[mask2] = (ebreak1 / Epiv) ** alpha * (x[mask2] / ebreak1) ** beta
mask3 = (x > ebreak2)
fxn[mask3] = (ebreak1 / Epiv) ** alpha * (ebreak1 / ebreak2) ** beta * \
(x[mask3] / ebreak2) ** gamma
return A * fxn
[docs]class SmoothlyBrokenPowerLaw(Function):
r"""A smoothly-broken power law with a break scale:
:math:`F(E) = A \bigl(\frac{E}{E_{piv}}\bigr)^\beta 10^{(\beta-\beta_{piv})}, \\
\text{where }\\
\beta = m \Delta \ln\bigl(\frac{e^\alpha + e^{-\alpha}}{2}\bigr); \text{ }
\beta_{piv} = m \Delta \ln\bigl(\frac{e^{\alpha_{piv}} + e^{-\alpha_{piv}}}{2}\bigr); \\
\alpha = \frac{\log(E/E_b)}{\Delta}; \text{ }
\alpha_{piv} = \frac{\log(E_{piv}/E_b)}{\Delta}; \\
m = \frac{\lambda_2 - \lambda_1}{2}; \text{ }
b = \frac{\lambda_1 + \lambda_2}{2};`
and where
* :math:`A` is the amplitude in ph/s/cm^2/keV
* :math:`E_{b}` is the break energy in keV
* :math:`\Delta` is the break scale in decades of energy
* :math:`\lambda_1` is the low-energy power-law index
* :math:`\lambda_2` is the high-energy power-law index
* :math:`E_{piv}` is the pivot energy in keV. Nominally fixed to 100 keV.
"""
nparams = 6
param_list = [('A', 'ph/s/cm^2/keV', 'Amplitude'),
('Eb', 'keV', 'Break Energy'),
('Delta', 'dex', 'Break Scale'),
('lambda1', '', 'Low-Energy Photon index'),
('lambda2', '', 'High-Energy Photon index'),
('Epiv', 'keV', 'Pivot energy')]
default_values = [0.01, 700.0, 0.3, -1.0, -2.0, 100.0]
delta_abs = [0.1, 0.1, 0.1, 0.1, 0.1, 0.1]
delta_rel = [0.01, 0.01, 0.01, 0.01, 0.01, 0.01]
min_values = [1e-10, 0.01, 0.01, -1.9, -10.0, 0.01]
max_values = [np.inf, np.inf, 10.0, 20.0, 20.0, np.inf]
free = [True, True, True, True, True, False]
def eval(self, params, x):
A, ebreak, delta, lam1, lam2, Epiv = params
m = (lam2 - lam1) / 2.0
b = (lam1 + lam2) / 2.0
apiv = np.log10(Epiv / ebreak) / delta
a = np.log10(x / ebreak) / delta
bpiv = m * delta * np.log((np.exp(apiv) + np.exp(-apiv)) / 2.0)
beta = m * delta * np.log((np.exp(a) + np.exp(-a)) / 2.0)
fxn = A * (x / Epiv) ** b * 10.0 ** (beta - bpiv)
return fxn
[docs]class LogNormal(Function):
r"""A Log-Normal function:
:math:`F(E) = \frac{A}{\sqrt{2\pi} \sigma E}
e^{-\bigl(\frac{\ln E - \mu}{2\sigma}\bigr)^2}`
where
* :math:`A` is the amplitude in ph/s/cm^2/keV
* :math:`\mu` is the natural log of energy
* :math:`\sigma` is the logarithmic standard deviation
"""
nparams = 3
param_list = [('A', 'ph/s/cm^2/keV', 'Amplitude'),
('mu', 'Ln(keV)', 'Ln(Energy)'),
('sigma', 'Ln(keV)', 'Log Standard Deviation')]
default_values = [0.01, 5.0, 1.0]
delta_abs = [0.1, 0.1, 0.1]
delta_rel = [0.01, 0.01, 0.01]
min_values = [1e-10, -1.0, 0.0]
max_values = [np.inf, 10.0, np.inf]
free = [True, True, True]
def eval(self, params, x):
A, mu, sig = params
fxn = A / (np.sqrt(2.0 * np.pi) * sig * x) * np.exp(
-0.5 * ((np.log(x) - mu) / sig) ** 2)
return fxn
[docs]class GaussianLog(Function):
r"""A Gaussian with log_10(E):
:math:`F(E) = \frac{A}{\sqrt{2\pi} \sigma}
e^{-\bigl(\frac{\log E-\log E_{cen} }{2\sigma}\bigr)^2}, \\
\text{ where } \\
\sigma = \frac{W}{2.35482}`
and where
* :math:`A` is the amplitude in ph/s/cm^2/keV
* :math:`E_{cen}` is the centroid energy in keV
* :math:`W` is the :math:`\log_{10}` FWHM at :math:`E_{cen}` in decades
of energy
"""
nparams = 3
param_list = [('A', 'ph/s/cm^2/keV', 'Amplitude'),
('Ecen', 'keV', 'Centroid'),
('W', 'dex', 'Log(FWHM)')]
default_values = [0.01, 100.0, 1.0]
delta_abs = [0.1, 0.1, 0.1]
delta_rel = [0.01, 0.01, 0.01]
min_values = [1e-10, 0.1, 0.1]
max_values = [np.inf, np.inf, np.inf]
free = [True, True, True]
def eval(self, params, x):
A, Ecen, w = params
sig = w / 2.35482
fxn = A / (np.sqrt(2.0 * np.pi) * sig) * \
np.exp(-0.5 * ((np.log10(x) - np.log10(Ecen)) / sig) ** 2)
return fxn
[docs]class GaussianLogVaryingFWHM(Function):
r"""A Gaussian with log_10(E) and linearly-varying FWHM:
:math:`F(E) = \frac{A}{\sqrt{2\pi} \sigma}
e^{-\bigl(\frac{\log E-\log E_{cen} }{2\sigma}\bigr)^2}, \\
\text{ where } \\
\sigma = \frac{W + s (\log E - \log E_{cen})}{2.35482}`
and where
* :math:`A` is the amplitude in ph/s/cm^2/keV
* :math:`E_{cen}` is the centroid energy in keV
* :math:`W` is the :math:`\log_{10}` FWHM at :math:`E_{cen}` in decades
of energy
* :math:`s` is the slope of :math:`W` with respect to :math:`\log_{10}(E)`
in decades per :math:`\log_{10}(E)`
"""
nparams = 4
param_list = [('A', 'ph/s/cm^2/keV', 'Amplitude'),
('Ecen', 'keV', 'Centroid'),
('W', 'dex', 'Log(FWHM)'),
('s', 'dex/log(keV)', 'Slope of W')]
default_values = [0.01, 100.0, 1.0, 0.0]
delta_abs = [0.1, 0.1, 0.1, 0.1]
delta_rel = [0.01, 0.01, 0.01, 0.1]
min_values = [1e-10, 0.1, 0.1, -1.0]
max_values = [np.inf, np.inf, np.inf, np.inf]
free = [True, True, True, True]
def eval(self, params, x):
A, Ecen, w, s = params
sig = (w + s * (np.log10(x) - np.log10(Ecen))) / 2.35482
fxn = A / (np.sqrt(2.0 * np.pi) * sig) * \
np.exp(-0.5 * ((np.log10(x) - np.log10(Ecen)) / sig) ** 2)
return fxn
[docs]class SunyaevTitarchuk(Function):
r"""Sunyaev-Titarchuk Comptonization spectra:
:math:`F(E) = A E^2 e^{-E} \int_0^\infty t^{n-1}
e^{-t \bigl(1+\frac{t}{E}\bigr)^{n+3}} \text{d}t, \\
\text{ where } \\
n = -\frac{3}{2} + \sqrt{\gamma + \frac{9}{4}}; \\
\gamma = \frac{511 \pi^2}{G kT \bigl(\tau + \frac{2}{3} \bigr)^2}`
and where
* :math:`A` is the amplitude in ph/s/cm^2/keV
* :math:`kT` is the electron energy in keV
* :math:`\tau` is the optical depth
* :math:`G` is the geometry factor, fixed at 3 for a spherical cloud and
at 12 for a disk of electrons
References:
`Sunyaev, R. A. and Titarchuk, L. G. 1980, A&A, 86, 121
<http://adsabs.harvard.edu/full/1980A%26A....86..121S>`_
"""
nparams = 4
param_list = [('A', 'ph/s/cm^2/keV', 'Amplitude'),
('kT', 'keV', 'Electron energy'),
('tau', '', 'Optical depth'),
('G', '', 'Geometry factor')]
default_values = [0.01, 30.0, 10.0, 3.0]
delta_abs = [0.1, 0.1, 0.1, 0.1]
delta_rel = [0.01, 0.01, 0.01, 0.01]
min_values = [1e-10, 0.1, 0.0, 3]
max_values = [np.inf, np.inf, np.inf, 12]
free = [True, True, True, False]
def eval(self, params, x):
A, kT, tau, G = params
gamma = 511.0 * np.pi ** 2 / (G * kT * (tau + 2. / 3.) ** 2)
n = -1.5 + np.sqrt(gamma + 2.25)
integrand = lambda t, n, x: t ** (n - 1.0) * np.exp(
-t * (1.0 + t / x) ** (n + 3.))
integral = [quad(integrand, 0, np.inf, args=(n, x1))[0] for x1 in x]
fxn = A * x ** 2 * np.exp(-x) * np.array(integral)
return fxn
[docs]class OTTB(Function):
r"""Optically-Thin Thermal Bremsstrahlung:
:math:`F(E) = A e^{-\frac{E}{kT}} e^{\frac{E_{piv}}{kT}}
\biggl(\frac{E_{piv}}{E} \biggr)`
where
* :math:`A` is the amplitude in ph/s/cm^2/keV
* :math:`kT` is the electron energy in keV
* :math:`E_{piv}` is the pivot energy in keV
"""
nparams = 3
param_list = [('A', 'ph/s/cm^2/keV', 'Amplitude'),
('kT', 'keV', 'Electron energy'),
('Epiv', 'keV', 'Pivot energy')]
default_values = [0.01, 30.0, 100.0]
delta_abs = [0.1, 0.1, 0.1]
delta_rel = [0.01, 0.01, 0.01]
min_values = [1e-10, 0.1, 0.0]
max_values = [np.inf, np.inf, np.inf]
free = [True, True, False]
def eval(self, params, x):
A, kT, Epiv = params
fxn = A * np.exp(-x / kT) * np.exp(Epiv / kT) * (Epiv / x)
return fxn
[docs]class BlackBody(Function):
r"""Black Body spectrum:
:math:`F(E) = A \frac{E^2}{e^{E/kT}-1}`
where
* :math:`A` is the amplitude in ph/s/cm^2/keV
* :math:`kT` is the temperature in keV
"""
nparams = 2
param_list = [('A', 'ph/s/cm^2/keV', 'Amplitude'),
('kT', 'keV', 'Temperature')]
default_values = [0.01, 30.0]
delta_abs = [0.1, 0.1]
delta_rel = [0.01, 0.01]
min_values = [1e-10, 0.1]
max_values = [np.inf, np.inf]
free = [True, True]
def eval(self, params, x):
A, kT = params
fxn = A * x ** 2 / (np.exp(x / kT) - 1.0)
return fxn
[docs]class YangSoongPulsar(Function):
r"""Yang Soong's pulsar spectral form
:math:`F(E) = A C (1 - E_W G), \\
\text{ where} \\
C = \begin{cases}
\bigl(\frac{E}{E_{piv}}\bigr)^{-\alpha}, & \text{ if } {E \leq E_b} \\
\frac{M}{E} e^{-E/E_{piv}}, & \text{ otherwise }
\end{cases}; \\
M = e^{E_b/E_F} E_b \bigl(\frac{E_b}{E_{piv}}\bigr)^{-\alpha}; \\
G = \frac{0.94}{W} e^{-2.76 (\frac{E-E_{cen}}{W})^2};`
and where
* :math:`A` is the amplitude in ph/s/cm^2/keV
* :math:`\alpha` is the power-law index
* :math:`E_b` is the break energy in keV
* :math:`E_F` is the e-folding energy in keV
* :math:`E_W` is the equivalent line width in keV
* :math:`E_{cen}` is the line centroid in keV
* :math:`W` is the FWHM of the line, in keV
* :math:`E_{piv}` is the pivot energy in keV
References:
`Soong, Y., et al. 1990, ApJ, 348, 641
<http://adsabs.harvard.edu/full/1990ApJ...348..641S>`_
"""
nparams = 8
param_list = [('A', 'ph/s/cm^2/keV', 'Amplitude'),
('alpha', '', 'Index'),
('Eb', 'keV', 'Break Energy'),
('EF', 'keV', 'E-folding Energy'),
('EW', 'keV', 'Line Width'),
('Ecen', 'keV', 'Line Centroid'),
('W', 'keV', 'Line FWHM'),
('Epiv', 'keV', 'Pivot Energy')]
default_values = [0.01, -1.0, 200.0, 300.0, 50.0, 200.0, 50.0, 100.0]
delta_abs = [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]
delta_rel = [0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01]
min_values = [1e-10, -10.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
max_values = [np.inf, 10.0, np.inf, np.inf, np.inf, np.inf, np.inf, np.inf]
free = [True, True, True, True, True, True, True, False]
def eval(self, params, x):
A, alpha, Eb, EF, EW, Ecen, W, Epiv = params
M = np.exp(Eb / EF) * Eb * (Eb / Epiv) ** (-alpha)
G = 0.94 / W * np.exp(-2.76 * ((x - Ecen) / W) ** 2)
C = np.zeros(len(x), dtype=float)
mask = (x <= Eb)
C[mask] = A * (x[mask] / Epiv) ** (-alpha)
C[~mask] = A * M * np.exp(-x[~mask] / Epiv) / x[~mask]
fxn = C * (1.0 - EW * G)
return fxn
[docs]class TanakaPulsar(Function):
r"""Tanaka's pulsar model
:math:`F(E) = A \frac{C}{C_{piv}} e^{L_{piv}-L}, \\
\text{ where} \\
C = A E^{-\alpha} e^{-E/kT}; \text{ }
C_{piv} E_{piv}^{-\alpha} e^{-E_{piv}/kT}; \\
L = \sum_{i=1}^{N_L} \frac{\tau_i (i W)^2 [E/(i E_L)]^2}{(E-i E_L)^2 + (i W)^2}; \\
L_{piv} = \sum_{i=1}^{N_L} \frac{\tau_i (i W)^2 [E_{piv}/(i E_L)]^2}
{(E_{piv}-i E_L)^2 + (i W)^2}; \\`
and where
* :math:`A` is the amplitude in ph/s/cm^2/keV
* :math:`\alpha` is the power-law index
* :math:`kT` is the temperature in keV
* :math:`\tau_1` is the optical depth of the first line
* :math:`\tau_2` is the optical depth of the second line
* :math:`N_L` is the number of lines; fixed at either 0, 1, or 2
* :math:`E_L` is the line centroid in keV of the first line
* :math:`W` is the line FWHM in keV
* :math:`E_{piv}` is the pivot energy in keV
References:
`Grove, J. E., et al. 1995, ApJ, 438, L25
<https://apps.dtic.mil/docs/citations/ADA461878>`_
"""
nparams = 9
param_list = [('A', 'ph/s/cm^2/keV', 'Amplitude'),
('alpha', '', 'Index'),
('kT', 'keV', 'Temperature'),
('tau1', '', 'Optical Depth of 1st line'),
('tau2', '', 'Optical Depth of 2nd line'),
('NL', '', 'Number of lines'),
('EL', 'keV', '1st Line Centroid'),
('W', 'keV', '1st Line FWHM'),
('Epiv', 'keV', 'Pivot Energy')]
default_values = [0.01, -1.0, 200.0, 1.0, 1.0, 1, 100.0, 50.0, 100.0]
delta_abs = [0.1, 0.1, 0.1, 0.1, 0.1, 0.0, 0.1, 0.1, 0.1]
delta_rel = [0.01, 0.01, 0.01, 0.01, 0.01, 0.00, 0.01, 0.01, 0.01]
min_values = [1e-10, -10.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
max_values = [np.inf, 10.0, np.inf, np.inf, np.inf, 2, np.inf, np.inf,
np.inf]
free = [True, True, True, True, True, False, True, True, False]
def eval(self, params, x):
A, alpha, kT, tau1, tau2, NL, EL, W, Epiv = params
C = A * x ** (-alpha) * np.exp(-x / kT)
Cpiv = Epiv ** (-alpha) * np.exp(-Epiv / kT)
if NL > 0:
taus = np.array((tau1, tau2))
idx = np.arange(2) + 1
L = (taus * (idx[np.newaxis, :] * W) ** 2 * \
(x[:, np.newaxis] / (idx[np.newaxis, :] * EL)) ** 2) / \
((x[:, np.newaxis] - idx[np.newaxis, :] * EL) ** 2 + (
idx[np.newaxis, :] * W) ** 2)
L = L[:, :NL].sum(axis=1)
Lpiv = (taus * (idx * W) ** 2 * (Epiv / (idx * EL)) ** 2) / (
(Epiv - idx * EL) ** 2 + (idx * W) ** 2)
Lpiv = Lpiv[:NL].sum()
else:
L = Lpiv = 0.0
fxn = A * C / Cpiv * np.exp(Lpiv - L)
return fxn
[docs]class OTTS(Function):
r"""Optically-Thin Thermal Synchrotron:
:math:`F(E) = A e^{(E/E_C)^{1/3}}`
where
* :math:`A` is the amplitude in ph/s/cm^2/keV
* :math:`E_C` is the energy scale in keV
References:
`Liang, E. P., et al., 1983, ApJ, 271, 776
<http://adsabs.harvard.edu/full/1983ApJ...271..766L>`_
"""
nparams = 2
param_list = [('A', 'ph/s/cm^2/keV', 'Amplitude'),
('EC', 'keV', 'Energy Scale')]
default_values = [0.01, 100.0]
delta_abs = [0.1, 0.1]
delta_rel = [0.01, 0.01]
min_values = [1e-10, 0]
max_values = [np.inf, np.inf]
free = [True, True]
def eval(self, params, x):
A, Ec = params
fxn = A * np.exp((x / Ec) ** (1.0 / 3.0))
return fxn
[docs]class GaussLine(Function):
r"""A Gaussian Line:
:math:`F(E) = A \frac{\text{erfc}(a_l)-\text{erfc}(a_r)}{2(E_r-E_l)}, \\
\text{ where } \\
a_l = \frac{E_l - E_{cen}}{\sqrt{2}\sigma}; \text{ }
a_r = \frac{E_r - E_{cen}}{\sqrt{2}\sigma}; \\
\sigma = W/2.35482`
and where
* :math:`A` is the amplitude in ph/s/cm^2/keV
* :math:`E_{cen}` is the centroid energy in keV
* :math:`W` is the FHWM in keV
"""
nparams = 3
param_list = [('A', 'ph/s/cm^2/keV', 'Amplitude'),
('Ecen', 'keV', 'Centroid Energy'),
('W', 'keV', 'FWHM')]
default_values = [0.01, 100.0, 8.0]
delta_abs = [0.1, 0.1, 0.1]
delta_rel = [0.01, 0.01, 0.01]
min_values = [1e-10, 0.0, 0.0]
max_values = [np.inf, np.inf, np.inf]
free = [True, True, True]
def eval(self, params, x):
A, Ecen, W = params
# calculate edges; need to add an ending edge
edges = np.asarray(x)
dE = np.log10(x[-1]) - np.log10(x[-2])
edges = np.append(edges, [x[-1] * 10.0 ** dE])
el = edges[:-1]
er = edges[1:]
sig = W / 2.35482
al = (el - Ecen) / (np.sqrt(2.0) * sig)
ar = (er - Ecen) / (np.sqrt(2.0) * sig)
fxn = A * (erfc(al) - erfc(ar)) / (2.0 * (er - el))
return fxn
[docs]class LowEnergyCutoff(Function):
r"""A Low-Energy cutoff (Multiplicative Component):
:math:`F(E) = \begin{cases}
(E/E_{cut})^{E_{cut}/E_F} e^{(E_{cut}-E)/E_F} & \text{ if } {E \leq E_{cut}} \\
1 & \text{ otherwise }
\end{cases}`
where
* :math:`E_{cut}` is the cutoff energy in keV
* :math:`E_{cen}` is the folding energy in keV
"""
nparams = 2
param_list = [('Ecut', 'keV', 'Cutoff Energy'),
('EF', 'keV', 'Folding Energy')]
default_values = [100.0, 10.0]
delta_abs = [0.1, 0.1]
delta_rel = [0.01, 0.01]
min_values = [0.0, 0.0]
max_values = [np.inf, np.inf]
free = [True, True]
def eval(self, params, x):
Ecut, EF = params
mask = (x <= Ecut)
fxn = np.ones(len(x), dtype=float)
fxn[mask] = (x[mask] / Ecut) ** (Ecut / EF) * np.exp(
(Ecut - x[mask]) / EF)
return fxn
[docs]class HighEnergyCutoff(Function):
r"""A High-Energy cutoff (Multiplicative Component):
:math:`F(E) = \begin{cases}
(E/E_{cut})^{E_{cut}/E_F} e^{(E_{cut}-E)/E_F} & \text{ if } {E > E_{cut}} \\
1 & \text{ otherwise }
\end{cases}`
where
* :math:`E_{cut}` is the cutoff energy in keV
* :math:`E_{cen}` is the folding energy in keV
"""
nparams = 2
param_list = [('Ecut', 'keV', 'Cutoff Energy'),
('EF', 'keV', 'Folding Energy')]
default_values = [1000.0, 100.0]
delta_abs = [0.1, 0.1]
delta_rel = [0.01, 0.01]
min_values = [0.0, 0.0]
max_values = [np.inf, np.inf]
free = [True, True]
def eval(self, params, x):
Ecut, EF = params
mask = (x > Ecut)
fxn = np.ones(len(x), dtype=float)
fxn[mask] = (x[mask] / Ecut) ** (Ecut / EF) * np.exp(
(Ecut - x[mask]) / EF)
return fxn
[docs]class PowerLawMult(Function):
r"""A Power Law (Multiplicative Component):
:math:`F(E) = \bigl(\frac{E}{E_{piv}} \bigr)^\lambda`
where
* :math:`\lambda` is the index
* :math:`E_{piv}` is the pivot energy in keV
"""
nparams = 2
param_list = [('lambda', '', 'Index'),
('Epiv', 'keV', 'Pivot Energy')]
default_values = [-2.0, 100.0]
delta_abs = [0.1, 0.1]
delta_rel = [0.01, 0.01]
min_values = [-10.0, 0.0]
max_values = [10.0, np.inf]
free = [True, False]
def eval(self, params, x):
lam, Epiv = params
fxn = (x / Epiv) ** lam
return fxn
[docs]class GaussLineMult(Function):
r"""A Gaussian Line (Multiplicative Component):
:math:`F(E) = e^{I G}, \\
\text{ where } \\
G = \frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{(E-E_{cen})^2}{2\sigma^2}}; \\
\sigma = W/2.35482`
and where
* :math:`I` is the intensity
* :math:`E_{cen}` is the centroid energy in keV
* :math:`W` is the FWHM in keV
"""
nparams = 3
param_list = [('I', '', 'Intensity'),
('Ecen', 'keV', 'Centroid Energy'),
('W', 'keV', 'FWHM')]
default_values = [1.0, 10.0, 8.0]
delta_abs = [0.1, 0.1, 0.1]
delta_rel = [0.01, 0.01, 0.01]
min_values = [0.0, 0.0, 0.0]
max_values = [np.inf, np.inf, np.inf]
free = [True, True, True]
def eval(self, params, x):
I, Ecen, W = params
sig = W / 2.35482
G = 1.0 / (np.sqrt(2.0 * np.pi) * sig) * np.exp(
-(x - Ecen) ** 2 / (2.0 * sig ** 2))
fxn = np.exp(I * G)
return fxn
[docs]class LorentzLineMult(Function):
r"""A Lorentzian Line (Multiplicative Component):
:math:`F(E) = e^{\frac{I}{W^2 + (E-E_{cen})^2}}`
where
* :math:`I` is the intensity
* :math:`E_{cen}` is the centroid energy in keV
* :math:`W` is the FWHM in keV
"""
nparams = 3
param_list = [('I', '', 'Intensity'),
('Ecen', 'keV', 'Centroid Energy'),
('W', 'keV', 'FWHM')]
default_values = [1.0, 10.0, 8.0]
delta_abs = [0.1, 0.1, 0.1]
delta_rel = [0.01, 0.01, 0.01]
min_values = [0.0, 0.0, 0.0]
max_values = [np.inf, np.inf, np.inf]
free = [True, True, True]
def eval(self, params, x):
I, Ecen, W = params
fxn = np.exp(I / (W ** 2 + (x - Ecen) ** 2))
return fxn