'''
A collection of spectrum functions.
Translate from/to MEGAlib's format
'''
import logging
logger = logging.getLogger(__name__)
from bctools.util import units as u
from bctools.exceptions import UserError,ParsingError
import numpy as np
from numpy import log
from math import isfinite
import sys
from abc import ABC, abstractmethod
from scipy.integrate import quad
[docs]
class Spectrum:
'''
Abstract class representing an spectral model
Child classes need to implement eval(), integrate(), _from_megalib(),
to_megalib(), __str__(), __eq__() and normalize() and property flux
'''
[docs]
@classmethod
def from_name(cls, name, *args, **kwargs):
'''
Instantiate a spectrum from the subclass name. e.g.:
:code:`Spectrum.from_name('PowerLaw', *args, **kwargs) <==> PowerLaw(*args, **kwargs)`
'''
try:
subclass = getattr(sys.modules[__name__], name)
except AttributeError:
raise ValueError("Class 'name' is not a subclass of {}".format(cls.__name__))
if not issubclass(subclass, cls):
raise ValueError("Class 'name' is not a subclass of {}".format(cls.__name__))
return subclass(*args, **kwargs)
[docs]
@abstractmethod
def eval(self, energy):
'''
Evaluate the spectrum F(E)
Args:
energy (iterable): energy :math:`E` values
Returns:
array: flux :math:`F` values
'''
[docs]
@abstractmethod
def integrate(self, bounds):
'''
Evaluate the integral between two energy bounds
Args:
bound (array): Energy bounds. Either an array with 2 elements or 2D
array with 2 columns
Returns:
float or array: Integral flux between bounds
'''
@abstractmethod
def __str__(self):
'''
String with F(E) equation and all parameters
Returns:
str
'''
@abstractmethod
def __eq__(self):
'''
Compare two spectra
Returns:
eq
'''
[docs]
@abstractmethod
def normalize(self, flux, *args, **kwargs):
"""
Change normalization parameters based on an integrated flux.
"""
@property
@abstractmethod
def flux(self):
"""
Total integrated flux between energy bounds
"""
[docs]
@classmethod
def from_megalib(cls, params, flux):
'''
Get object from MEGAlibs format.
Args:
params (list): List of paramters. See `str_megalib_format`
flux (str): Flux value (In MEGAlib the spectrum refers only to the shape)
'''
# Note that MEGAlibs uses same units as bc-tools (keV, s, cm)
if len(params) < 1:
raise ParsingError("Need at least 2 paramter")
stype = params[0]
if stype == "PowerLaw":
return PowerLaw._from_megalib(params, flux)
elif stype == "Mono":
return Mono._from_megalib(params, flux)
elif stype == "BandFunction":
return BandFunction._from_megalib(params, flux)
elif stype == "Comptonized":
return Comptonized._from_megalib(params, flux)
else:
raise ValueError("{} spectrum not yet implemented".format(stype))
@classmethod
@abstractmethod
def _from_megalib(cls, params, flux):
"""
Actual implementation or MEGAlin spectrum parsing. To be overridden by
spectrum type
"""
[docs]
@abstractmethod
def to_megalib(self):
'''
Return equivalent parameters of MEGAlib's Source.Spectrum and Source.Flux
Returns:
list of str: Spectrum type and parameteters for Source.Spectrum
str: Parameter for Source.Flux
'''
[docs]
def set_elims(self, min_energy=None, max_energy=None):
'''
Set energy limits. Function evaluates to 0 outside of these bounds
Limits stay the same if not specified
Args:
min_energy (float): lower energy bound
max_energy (float): uppper energy bound
'''
if min_energy is not None:
self.min_energy = min_energy
if max_energy is not None:
self.max_energy = max_energy
[docs]
class PowerLaw(Spectrum):
r'''
:math:`F(E) = A \bigl(\frac{E}{1 \mathrm{keV}}\bigr)^{-\alpha}`
Args:
norm (float): normalization :math:`A`
index (float): spectral index :math:`\alpha`
Attributes:
norm (float): normalization :math:`A`
index (float): spectral index :math:`\alpha`
min_energy (float): :math:`F=0` if :math:`E<` `min_energy`
max_energy (float): :math:`F=0` if :math:`E>` `max_energy`
'''
def __init__(self, norm, index):
self.min_energy = 0
self.max_energy = u.inf
self.norm = norm
self.index = index
def __str__(self):
return ("F(E) = {:.2E}(E/1keV)^-{:.1f} keV^-1 cm^-2 s^-1. "
"Range: ({:.2E} keV - {:.2E} keV)").format(self.norm,
self.index,
self.min_energy,
self.max_energy)
def __eq__(self, other):
return (self.norm == other.norm and
self.index == other.index and
self.min_energy == other.min_energy and
self.max_energy == other.max_energy)
[docs]
def eval(self, energy):
'''
Evaluate the spectrum F(E)
Args:
energy (iterable): energy :math:`E` values
Returns:
array: flux :math:`F` values. F = 0 outside of the min/max_energy bounds
'''
flux = self.norm * np.array(energy)**(-self.index)
flux[energy < self.min_energy] = 0
flux[energy > self.max_energy] = 0
return flux
[docs]
def integrate(self, bounds):
'''
Evaluate the integral between two energy bounds.
Note that the spectrum evaluates to 0 out of the overall min/max_energy
bounds.
Args:
bound (array): Energy bounds. Either an array with 2 elements or 2D
array with 2 columns
Returns:
float or array: Integral flux between bounds
'''
shape = np.shape(bounds)
if shape == (2,):
# Single bound
min_energy = bounds[0]
max_energy = bounds[1]
if min_energy > max_energy:
raise ValueError("Upper bound must be greater than lower bound")
if max_energy <= self.min_energy or min_energy >= self.max_energy:
return 0
min_energy = max(min_energy, self.min_energy)
max_energy = min(max_energy, self.max_energy)
if self.index == 1:
return self.norm * log(max_energy/min_energy)
else:
return (self.norm
* (max_energy**(1-self.index)
- min_energy**(1-self.index))
/ (1-self.index))
elif len(shape) > 1 and shape[1] == 2:
# Multiple bounds, call recursively
return np.array([self.integrate(erange) for erange in bounds])
else:
raise ValueError("Shape of bounds can either be (2,) or (N,2)")
[docs]
def normalize(self, flux, min_energy = None, max_energy = None):
"""
Change norm based on the total emission between
min_enegy and max_energy
Args:
flux (float): new integral flux between min_energy and max_energy
min_energy (float): Lower integration bound. Defaults to self.min_energy
max_energy (float): Upper integration bound. Defaults to self.max_energy
"""
if min_energy is None:
min_energy = self.min_energy
if max_energy is None:
max_energy = self.max_energy
self.norm = flux / self.integrate([min_energy, max_energy])
@property
def flux(self):
return self.integrate([self.min_energy, self.max_energy])
@classmethod
def _from_megalib(cls, params, flux):
'''
Initialize from MEGAlibs format (PowerLaw spectrum)
Args:
params (list): List of parameters. First parameter must be "PowerLaw"
See `str_megalib_format`
flux (str): Flux value (In MEGAlib the spectrum refers only to the shape)
'''
if len(params) != 4:
raise ParsingError("PowerLaw needs 3 parameters")
if params[0] != "PowerLaw":
raise ParsingError("This doesn't seem like a power law")
flux = float(flux)
index = float(params[3])
min_energy = float(params[1])
max_energy = float(params[2])
obj = PowerLaw(1, index)
obj.set_elims(min_energy, max_energy)
obj.normalize(flux)
return obj
[docs]
def to_megalib(self):
'''
Return equivalent parameters of MEGAlib's Source.Spectrum and Source.Flux
Returns:
list of str: Spectrum type and parameteters for Source.Spectrum \
(i.e. PowerLaw min_energy max_energy photon_index)
str: Parameter for Source.Flux
'''
#MEGAlib flux is the integral between min_energy - max_energy
params = ["PowerLaw"]
params += ["{:.17E}".format(self.min_energy)]
if not isfinite(self.max_energy):
raise UserError("Set finite energy limits by calling set_elims() "
"before calling to_megalib()")
params += ["{:.17E}".format(self.max_energy)]
params += ["{:.17E}".format(self.index)]
# MEGAlib's flux is between min_energy and max_energy
flux = "{:.17E}".format(self.integrate([self.min_energy, self.max_energy])
/u.s/u.cm2)
return params,flux
[docs]
class Mono(Spectrum):
"""
Mono energetic beam
Args:
flux (float): Flux
energy (float): Energy
"""
def __init__(self, flux, energy):
self._flux = flux
self.energy = energy
self.min_energy = 0
self.max_energy = u.inf
def __str__(self):
return f"{self._flux:.2E} cm^-2 s^-1 at {self.energy:.2E} keV"
def __eq__(self, other):
return (isinstance(other, Mono) and
self._flux == other._flux and
self.energy == other.energy)
[docs]
def eval(self, energy):
energy = np.array(energy)
non_zero = (energy >= self.min_enegy and
energy < self.max_energy and
energy == self.energy)
flux = np.zeros(energy.shape)
flux[non_zero] = self._flux
[docs]
def integrate(self, bounds):
'''
Evaluate the integral between two energy bounds
Args:
bound (array): Energy bounds. Either an array with 2 elements or 2D
array with 2 columns
Returns:
float or array: Integral flux between bounds
'''
shape = np.shape(bounds)
if shape == (2,):
# Single bound
min_energy = bounds[0]
max_energy = bounds[1]
if min_energy > max_energy:
raise ValueError("Upper bound must be greater than lower bound")
if (self.min_energy <= self.energy and
self.energy < self.max_energy and
min_energy <= self.energy and
self.energy < max_energy):
return self._flux
else:
return 0
elif len(shape) > 1 and shape[1] == 2:
# Multiple bounds, call recursively
return np.array([self.integrate(erange) for erange in bounds])
else:
raise ValueError("Shape of bounds can either be (2,) or (N,2)")
[docs]
def normalize(self, flux):
"""
Update the flux
Args:
flux (float): new flux
"""
self._flux = flux
@classmethod
def _from_megalib(cls, params, flux):
'''
Initialize from MEGAlibs format (Mono spectrum)
Args:
params (list): List of parameters. First parameter must be "PowerLaw"
See `str_megalib_format`
flux (str): Flux value (In MEGAlib the spectrum refers only to the shape)
'''
if len(params) != 2:
raise ParsingError("Mono needs 2 parameters")
if params[0] != "Mono":
raise ParsingError("This doesn't seem like a mono-energetic beam")
flux = float(flux)
energy = float(params[1])
obj = cls(flux, energy)
return obj
[docs]
def to_megalib(self):
'''
Return equivalent parameters of MEGAlib's Source.Spectrum and Source.Flux
Returns:
list of str: Spectrum type and parameteters for Source.Spectrum \
(i.e. PowerLaw min_energy max_energy photon_index)
str: Parameter for Source.Flux
'''
#MEGAlib flux is the integral between min_energy - max_energy
params = ["Mono"]
params += ["{:.17E}".format(self.energy)]
flux = "{:.17E}".format(self._flux/u.s/u.cm2)
return params,flux
[docs]
class BandFunction(Spectrum):
r'''
:math:`F(E) =
A
\begin{cases}
\left(\frac{E}{1 \mathrm{keV}}\right)^{\alpha_1} \cdot e^{-E/E_{\text{break}}}, & E < E_{\text{trans}} \\
\left[ (\alpha_1-\alpha_2) \cdot \left(\frac{E_{\text{break}}}{1 \mathrm{keV}}\right) \right]^{(\alpha_1-\alpha_2)} \cdot e^{(\alpha_2-\alpha_1)} \cdot
\left(\frac{E}{1 \mathrm{keV}}\right)^{\alpha_2} , & E \geq E_{\text{trans}}
\end{cases}`
:math:`E_{\text{trans}} = (\alpha_1 - \alpha_2) \cdot E_{\text{break}}`
Args:
norm (float): normalization :math:`A`
lower_index (float): spectral index if :math:`E<` `transition_energy` :math:`\alpha_1`
upper_index (float): spectral index if :math:`E>` `transition_energy` :math:`\alpha_2`
break_energy (float): scale energy of the exponential cutoff :math:`E_{\text{break}}`
Attributes:
norm (float): normalization :math:`A`
lower_index (float): spectral index if :math:`E<` `transition_energy` :math:`\alpha_1`
upper_index (float): spectral index if :math:`E>` `transition_energy` :math:`\alpha_2`
break_energy (float): scale energy of the exponential cutoff :math:`E_{\text{break}}`
transition_energy (float): transition value between the low- and high-energy ranges :math:`E_{\text{transition}}`
min_energy (float): :math:`F=0` if :math:`E<` `min_energy`
max_energy (float): :math:`F=0` if :math:`E>` `max_energy`
'''
def __init__(self, norm, lower_index, upper_index, break_energy):
self.min_energy = 0
self.max_energy = u.inf
self.norm = norm
self.lower_index = lower_index
self.upper_index = upper_index
self.break_energy= break_energy
def __str__(self):
return (f"F(E) = {self.norm:.2E} (E/1keV)^{self.lower_index:.1f} e^(-E/{self.break_energy:.2E}) keV^-1 cm^-2 s^-1. "
f"Range: ({self.min_energy:.2E} keV - {self.transition_energy:.2E} keV)\n"
f"F(E) = {self.norm:.2E} [({self.lower_index:.1f}-{self.upper_index:.1f}) ({self.break_energy:.2E}/1keV)]^({self.lower_index:.1f}-{self.upper_index:.1f}) e^({self.upper_index:.1f}-{self.lower_index:.1f}) (E/1keV)^{self.upper_index:.1f} keV^-1 cm^-2 s^-1. "
f"Range: ({self.transition_energy:.2E} keV - {self.max_energy:.2E} keV)\n"
)
def __eq__(self, other):
return (self.norm == other.norm and
self.lower_index == other.lower_index and
self.upper_index == other.upper_index and
self.break_energy == other.break_energy and
self.min_energy == other.min_energy and
self.max_energy == other.max_energy
)
def _eval_shape(self, energy):
'''
Evaluate the spectral shape (no normalization factor).
Args:
energy (iterable): energy :math:`E` values in keV.
Returns:
array: values of the spectral shape evaluation. F = 0 outside of the min/max_energy bounds
'''
value = np.where(energy <= self.transition_energy,
np.power(energy, self.lower_index) * np.exp(- energy / self.break_energy),
np.power(energy, self.upper_index) * np.power(self.break_energy*(self.lower_index-self.upper_index), self.lower_index-self.upper_index) * np.exp(self.upper_index-self.lower_index)
)
value[energy < self.min_energy] = 0
value[energy > self.max_energy] = 0
return value
[docs]
def eval(self, energy):
'''
Evaluate the spectrum F(E)=dN/dE (E) in keV^-1 cm^-2 s^-1.
Args:
energy (iterable): energy :math:`E` values in keV.
Returns:
array: dNdE :math:`F` values. F = 0 outside of the min/max_energy bounds
'''
dNdE = self.norm * self._eval_shape(energy)
return dNdE
[docs]
def integrate(self, bounds, func=None):
'''
Evaluate the integral between two energy bounds.
Note that the spectrum evaluates to 0 out of the overall min/max_energy bounds.
Args:
bound (array): Energy bounds in keV. Either an array with 2 elements or 2D
array with 2 columns
func : function
The function to integrate. If default (func=None), it is self.eval.
Returns:
float or array: Integral photon flux between bounds in cm^-2 s^-1.
'''
shape = np.shape(bounds)
if shape == (2,):
# Single bound
min_energy = bounds[0]
max_energy = bounds[1]
if min_energy > max_energy:
raise ValueError("Upper bound must be greater than lower bound")
if max_energy <= self.min_energy or min_energy >= self.max_energy:
return 0
min_energy = max(min_energy, self.min_energy)
max_energy = min(max_energy, self.max_energy)
# Numerical Integral
if func==None:
func = self.eval # Default is used to compute photon flux
integral, error = quad(func, min_energy, max_energy)
return integral
elif len(shape) > 1 and shape[1] == 2:
# Multiple bounds, call recursively
return np.array([self.integrate(erange) for erange in bounds])
else:
raise ValueError("Shape of bounds can either be (2,) or (N,2)")
[docs]
def normalize(self, flux, min_energy = None, max_energy = None):
"""
Change norm based on the total emission between min_enegy and max_energy.
Args:
flux (float): new integral flux between min_energy and max_energy in cm^-2 s^-1.
min_energy (float): Lower integration bound in keV. Defaults to self.min_energy.
max_energy (float): Upper integration bound in keV. Defaults to self.max_energy.
"""
if min_energy is None:
min_energy = self.min_energy
if max_energy is None:
max_energy = self.max_energy
self.norm = flux / self.integrate([min_energy, max_energy], func=self._eval_shape)
@property
def flux(self):
"""
Returns the photon flux in cm^-2 s^-1.
"""
return self.integrate([self.min_energy, self.max_energy])
@classmethod
def _from_megalib(cls, params, flux):
'''
Initialize from MEGAlibs format (BandFunction spectrum)
Args:
params (list): List of parameters. First parameter must be "BandFunction"
See `str_megalib_format`
flux (str): Flux value (In MEGAlib the spectrum refers only to the shape)
'''
if len(params) != 6:
raise ParsingError("BandFunction needs 5 parameters")
if params[0] != "BandFunction":
raise ParsingError("This doesn't seem like a Band Function")
flux = float(flux)
min_energy = float(params[1])
max_energy = float(params[2])
lower_index = float(params[3])
upper_index = float(params[4])
break_energy= float(params[5])
obj = BandFunction(1, lower_index, upper_index, break_energy)
obj.set_elims(min_energy, max_energy)
obj.normalize(flux)
return obj
[docs]
def to_megalib(self):
'''
Return equivalent parameters of MEGAlib's Source.Spectrum and Source.Flux
Returns:
list of str: Spectrum type and parameteters for Source.Spectrum \
(i.e. BandFunction min_energy max_energy lower_index upper_index break_energy)
str: Parameter for Source.Flux
'''
#MEGAlib flux is the integral between min_energy - max_energy
params = ["BandFunction"]
params += ["{:.17E}".format(self.min_energy)]
if not isfinite(self.max_energy):
raise UserError("Set finite energy limits by calling set_elims() "
"before calling to_megalib()")
params += ["{:.17E}".format(self.max_energy)]
params += ["{:.17E}".format(self.lower_index)]
params += ["{:.17E}".format(self.upper_index)]
params += ["{:.17E}".format(self.break_energy)]
# MEGAlib's flux is between min_energy and max_energy
flux = "{:.17E}".format(self.integrate([self.min_energy, self.max_energy])
/u.s/u.cm2)
return params,flux
def __repr__(self):
return self.__str__()
@property
def transition_energy(self):
"""Energy of the exponential cutoff"""
return (self.lower_index-self.upper_index)*self.break_energy
[docs]
def eval_E2dNdE(self, energy):
"""
Evaluate the energy spectrum E**2 * dN/dE (E) in erg cm^-2 s^-1.
Args:
energy (iterable): energy :math:`E` values in keV.
Returns:
array: dNdE :math:`F` values. F = 0 outside of the min/max_energy bounds.
"""
# Spectrum in keV cm^-2 s^-1.
E2dNdE = np.power(energy,2) * self.eval(energy)
# Convert in erg cm^-2 s^-1.
E2dNdE/=u.erg
return E2dNdE
@property
def eflux(self):
"""
Returns the energy flux in erg cm^-2 s^-1.
"""
# Eflux in keV cm^-2 s^-1.
eflux = self.integrate([self.min_energy, self.max_energy], func=lambda energy: energy * self.eval(energy))
# Convert in erg cm^-2 s^-1.
eflux/=u.erg
return eflux
[docs]
class Comptonized(Spectrum):
r'''
:math:`F(E) = A \bigl(\frac{E}{1 \mathrm{keV}}\bigr)^{\alpha} \cdot e^{- E (2+\alpha)/E_{\text{peak}}}`
:math:`E_{\text{cutoff}} = E_{\text{peak}} / (2+\alpha)`
Args:
norm (float): normalization :math:`A`
alpha (float): spectral index :math:`\alpha`
peak_energy (float): peak energy before exponential cutoff :math:`E_{\text{peak}}`
Attributes:
norm (float): normalization :math:`A`
alpha (float): spectral index :math:`\alpha`
peak_energy (float): peak energy before exponential cutoff :math:`E_{\text{peak}}`
cutoff_energy (float): scale energy of the exponential cutoff :math:`E_{\text{cutoff}}`
min_energy (float): :math:`F=0` if :math:`E<` `min_energy`
max_energy (float): :math:`F=0` if :math:`E>` `max_energy`
'''
def __init__(self, norm, index, peak_energy):
self.min_energy = 0
self.max_energy = u.inf
self.norm = norm
self.index = index
self.peak_energy = peak_energy
def __str__(self):
return (f"F(E) = {self.norm:.2E} (E/1keV)^{self.index:.1f} e^(-E ({2+self.index:.1f})/{self.peak_energy:.2E}) keV^-1 cm^-2 s^-1. "
f"Range: ({self.min_energy:.2E} keV - {self.max_energy:.2E} keV)\n"
f"Cutoff Energy: {self.cutoff_energy:.2E} keV\n"
)
def __eq__(self, other):
return (self.norm == other.norm and
self.index == other.index and
self.peak_energy == other.peak_energy and
self.min_energy == other.min_energy and
self.max_energy == other.max_energy
)
def _eval_shape(self, energy):
'''
Evaluate the spectral shape (no normalization factor).
Args:
energy (iterable): energy :math:`E` values in keV.
Returns:
array: values of the spectral shape evaluation. F = 0 outside of the min/max_energy bounds
'''
energy = np.atleast_1d(energy) # Ensure energy is at least a 1D array, not a float
value = np.power(energy, self.index) * np.exp(- energy * (2+self.index) / self.peak_energy)
value[energy < self.min_energy] = 0
value[energy > self.max_energy] = 0
return value
[docs]
def eval(self, energy):
'''
Evaluate the spectrum F(E)=dN/dE (E) in keV^-1 cm^-2 s^-1.
Args:
energy (iterable): energy :math:`E` values in keV.
Returns:
array: dNdE :math:`F` values. F = 0 outside of the min/max_energy bounds
'''
dNdE = self.norm * self._eval_shape(energy)
return dNdE
[docs]
def integrate(self, bounds, func=None):
'''
Evaluate the integral between two energy bounds.
Note that the spectrum evaluates to 0 out of the overall min/max_energy bounds.
Args:
bound (array): Energy bounds in keV. Either an array with 2 elements or 2D
array with 2 columns
func : function
The function to integrate. If default (func=None), it is self.eval.
Returns:
float or array: Integral photon flux between bounds in cm^-2 s^-1.
'''
shape = np.shape(bounds)
if shape == (2,):
# Single bound
min_energy = bounds[0]
max_energy = bounds[1]
if min_energy > max_energy:
raise ValueError("Upper bound must be greater than lower bound")
if max_energy <= self.min_energy or min_energy >= self.max_energy:
return 0
min_energy = max(min_energy, self.min_energy)
max_energy = min(max_energy, self.max_energy)
# Numerical Integral
if func==None:
func = self.eval # Default is used to compute photon flux
integral, error = quad(func, min_energy, max_energy)
return integral
elif len(shape) > 1 and shape[1] == 2:
# Multiple bounds, call recursively
return np.array([self.integrate(erange) for erange in bounds])
else:
raise ValueError("Shape of bounds can either be (2,) or (N,2)")
[docs]
def normalize(self, flux, min_energy = None, max_energy = None):
"""
Change norm based on the total emission between min_enegy and max_energy.
Args:
flux (float): new integral flux between min_energy and max_energy in cm^-2 s^-1.
min_energy (float): Lower integration bound in keV. Defaults to self.min_energy.
max_energy (float): Upper integration bound in keV. Defaults to self.max_energy.
"""
if min_energy is None:
min_energy = self.min_energy
if max_energy is None:
max_energy = self.max_energy
self.norm = flux / self.integrate([min_energy, max_energy], func=self._eval_shape)
@property
def flux(self):
"""
Returns the photon flux in cm^-2 s^-1.
"""
return self.integrate([self.min_energy, self.max_energy])
@classmethod
def _from_megalib(cls, params, flux):
'''
Initialize from MEGAlibs format (Comptonized spectrum)
Args:
params (list): List of parameters. First parameter must be "Comptonized"
See `str_megalib_format`
flux (str): Flux value (In MEGAlib the spectrum refers only to the shape)
'''
if len(params) != 5:
raise ParsingError("Comptonized needs 5 parameters")
if params[0] != "Comptonized":
raise ParsingError("This doesn't seem like a Comptonized Function")
flux = float(flux)
min_energy = float(params[1])
max_energy = float(params[2])
index = float(params[3])
peak_energy= float(params[4])
obj = Comptonized(1, index, peak_energy)
obj.set_elims(min_energy, max_energy)
obj.normalize(flux)
return obj
[docs]
def to_megalib(self):
'''
Return equivalent parameters of MEGAlib's Source.Spectrum and Source.Flux
Returns:
list of str: Spectrum type and parameteters for Source.Spectrum \
(i.e. Comptonized min_energy max_energy index peak_energy)
str: Parameter for Source.Flux
'''
#MEGAlib flux is the integral between min_energy - max_energy
params = ["Comptonized"]
params += ["{:.17E}".format(self.min_energy)]
if not isfinite(self.max_energy):
raise UserError("Set finite energy limits by calling set_elims() "
"before calling to_megalib()")
params += ["{:.17E}".format(self.max_energy)]
params += ["{:.17E}".format(self.index)]
params += ["{:.17E}".format(self.peak_energy)]
# MEGAlib's flux is between min_energy and max_energy
flux = "{:.17E}".format(self.integrate([self.min_energy, self.max_energy])
/u.s/u.cm2)
return params,flux
def __repr__(self):
return self.__str__()
[docs]
def eval_E2dNdE(self, energy):
"""
Evaluate the energy spectrum E**2 * dN/dE (E) in erg cm^-2 s^-1.
Args:
energy (iterable): energy :math:`E` values in keV.
Returns:
array: dNdE :math:`F` values. F = 0 outside of the min/max_energy bounds.
"""
# Spectrum in keV cm^-2 s^-1.
E2dNdE = np.power(energy,2) * self.eval(energy)
# Convert in erg cm^-2 s^-1.
E2dNdE/=u.erg
return E2dNdE
@property
def eflux(self):
"""
Returns the energy flux in erg cm^-2 s^-1.
"""
# Eflux in keV cm^-2 s^-1.
eflux = self.integrate([self.min_energy, self.max_energy], func=lambda energy: energy * self.eval(energy))
# Convert in erg cm^-2 s^-1.
eflux/=u.erg
return eflux
@property
def cutoff_energy(self):
"""Energy of the exponential cutoff"""
return self.peak_energy/(2+self.index)