Source code for gammapy.utils.integrate

# Licensed under a 3-clause BSD style license - see LICENSE.rst
import numpy as np
from .interpolation import LogScale

__all__ = ["evaluate_integral_pwl", "trapz_loglog"]


[docs]def evaluate_integral_pwl(emin, emax, index, amplitude, reference): """Evaluate pwl integral (static function).""" val = -1 * index + 1 prefactor = amplitude * reference / val upper = np.power((emax / reference), val) lower = np.power((emin / reference), val) integral = prefactor * (upper - lower) mask = np.isclose(val, 0) if mask.any(): integral[mask] = (amplitude * reference * np.log(emax / emin))[mask] return integral
[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 : array_like Input array to integrate. x : array_like, optional Independent variable to integrate over. axis : int, optional Specify the axis. Returns ------- trapz : float Definite integral as approximated by trapezoidal rule in loglog space. """ # see https://stackoverflow.com/a/56840428 x, y = np.moveaxis(x, axis, 0), np.moveaxis(y, axis, 0) emin, emax = x[:-1], x[1:] vals_emin, vals_emax = y[:-1], y[1:] # log scale has the build-in zero clipping log = LogScale() index = -log(vals_emin / vals_emax) / log(emin / emax) index[np.isnan(index)] = np.inf return evaluate_integral_pwl( emin=emin, emax=emax, index=index, reference=emin, amplitude=vals_emin )