# Licensed under a 3-clause BSD style license - see LICENSE.rst
import numpy as np
from .interpolation import LogScale
__all__ = ["trapz_loglog"]
[docs]def trapz_loglog(y, x, axis=-1):
"""Integrate using the composite trapezoidal rule in log-log space.
Integrate `y` (`x`) along given axis in loglog space.
Parameters
----------
y : `~numpy.ndarray`
Input array to integrate.
x : `~numpy.ndarray`
Independent variable to integrate over.
axis : int, optional
Specify the axis. Default is -1.
Returns
-------
trapz : float
Definite integral as approximated by trapezoidal rule in loglog space.
"""
from gammapy.modeling.models import PowerLawSpectralModel as pl
# see https://stackoverflow.com/a/56840428
x, y = np.moveaxis(x, axis, 0), np.moveaxis(y, axis, 0)
energy_min, energy_max = x[:-1], x[1:]
vals_energy_min, vals_energy_max = y[:-1], y[1:]
# log scale has the built-in zero clipping
log = LogScale()
with np.errstate(invalid="ignore", divide="ignore"):
index = -log(vals_energy_min / vals_energy_max) / log(energy_min / energy_max)
index[np.isnan(index)] = np.inf
return pl.evaluate_integral(
energy_min=energy_min,
energy_max=energy_max,
index=index,
reference=energy_min,
amplitude=vals_energy_min,
)