# Licensed under a 3-clause BSD style license - see LICENSE.rst
import numpy as np
from astropy.units import Quantity
__all__ = ["integrate_spectrum"]
[docs]def integrate_spectrum(func, xmin, xmax, ndecade=100, intervals=False):
"""Integrate 1d function using the log-log trapezoidal rule.
If scalar values for xmin and xmax are passed an oversampled grid is generated using the
``ndecade`` keyword argument. If xmin and xmax arrays are passed, no
oversampling is performed and the integral is computed in the provided
grid.
Parameters
----------
func : callable
Function to integrate.
xmin : `~astropy.units.Quantity` or array-like
Integration range minimum
xmax : `~astropy.units.Quantity` or array-like
Integration range minimum
ndecade : int, optional
Number of grid points per decade used for the integration.
Default : 100.
intervals : bool, optional
Return integrals in the grid not the sum, default: False
"""
is_quantity = False
if isinstance(xmin, Quantity):
unit = xmin.unit
xmin = xmin.value
xmax = xmax.to_value(unit)
is_quantity = True
if np.isscalar(xmin):
logmin = np.log10(xmin)
logmax = np.log10(xmax)
n = int((logmax - logmin) * ndecade)
x = np.logspace(logmin, logmax, n)
else:
x = np.append(xmin, xmax[-1])
if is_quantity:
x = x * unit
y = func(x)
val = _trapz_loglog(y, x, intervals=intervals)
return val
# This function is copied over from https://github.com/zblz/naima/blob/master/naima/utils.py#L261
# and slightly modified to allow use with the uncertainties package
def _trapz_loglog(y, x, axis=-1, intervals=False):
"""Integrate using the composite trapezoidal rule in log-log space.
Integrate `y` (`x`) along given axis in loglog space.
Parameters
----------
y : array_like
Input array to integrate.
x : array_like, optional
Independent variable to integrate over.
axis : int, optional
Specify the axis.
intervals : bool, optional
Return array of shape x not the total integral, default: False
Returns
-------
trapz : float
Definite integral as approximated by trapezoidal rule in loglog space.
"""
log10 = np.log10
try:
y_unit = y.unit
y = y.value
except AttributeError:
y_unit = 1.0
try:
x_unit = x.unit
x = x.value
except AttributeError:
x_unit = 1.0
y = np.asanyarray(y)
x = np.asanyarray(x)
slice1 = [slice(None)] * y.ndim
slice2 = [slice(None)] * y.ndim
slice1[axis] = slice(None, -1)
slice2[axis] = slice(1, None)
slice1, slice2 = tuple(slice1), tuple(slice2)
# arrays with uncertainties contain objects
if y.dtype == "O":
from uncertainties.unumpy import log10
# uncertainties.unumpy.log10 can't deal with tiny values see
# https://github.com/gammapy/gammapy/issues/687, so we filter out the values
# here. As the values are so small it doesn't affect the final result.
# the sqrt is taken to create a margin, because of the later division
# y[slice2] / y[slice1]
valid = y > np.sqrt(np.finfo(float).tiny)
x, y = x[valid], y[valid]
if x.ndim == 1:
shape = [1] * y.ndim
shape[axis] = x.shape[0]
x = x.reshape(shape)
with np.errstate(invalid="ignore", divide="ignore"):
# Compute the power law indices in each integration bin
b = log10(y[slice2] / y[slice1]) / log10(x[slice2] / x[slice1])
# if local powerlaw index is -1, use \int 1/x = log(x); otherwise use normal
# powerlaw integration
trapzs = np.where(
np.abs(b + 1.0) > 1e-10,
(y[slice1] * (x[slice2] * (x[slice2] / x[slice1]) ** b - x[slice1]))
/ (b + 1),
x[slice1] * y[slice1] * np.log(x[slice2] / x[slice1]),
)
tozero = (y[slice1] == 0.0) + (y[slice2] == 0.0) + (x[slice1] == x[slice2])
trapzs[tozero] = 0.0
if intervals:
return trapzs * x_unit * y_unit
ret = np.add.reduce(trapzs, axis) * x_unit * y_unit
return ret