# Licensed under a 3-clause BSD style license - see LICENSE.rstimportnumpyasnpfrom.interpolationimportLogScale__all__=["trapz_loglog"]
[docs]deftrapz_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. """fromgammapy.modeling.modelsimportPowerLawSpectralModelaspl# see https://stackoverflow.com/a/56840428x,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 clippinglog=LogScale()withnp.errstate(invalid="ignore",divide="ignore"):index=-log(vals_energy_min/vals_energy_max)/log(energy_min/energy_max)index[np.isnan(index)]=np.infreturnpl.evaluate_integral(energy_min=energy_min,energy_max=energy_max,index=index,reference=energy_min,amplitude=vals_energy_min,)