Merge pull request #54 from lenasaurenhaus/master

Power law spectrum with subexponential cutoff

Author: cescalara

icecube_tools/source/flux_model.py

@@ -5,6 +5,7 @@
 from .power_law import BoundedPowerLaw
 from .power_law import BrokenBoundedPowerLaw
 from .power_law import BoundedPowerLawExpCutoff
+from .power_law import BoundedPowerLawSubexpCutoff
 
 """
 Module for simple flux models used in 
@@ -477,3 +478,91 @@ def sample(self, N):
         :param N: Number of samples.
         """
         return self.power_law.samples(N)
+    
+
+class PowerLawSubexpCutoffFlux(FluxModel):
+    """
+    Power law flux models with a subexponential cutoff.
+    """
+
+    def __init__(
+            self,
+            normalisation,
+            norm_energy,
+            index1,
+            cutoff_energy,
+            index2,
+            lower_energy=1e2,
+            upper_energy=1e8
+    ):
+        """
+        Power law flux models with a subexponential cutoff.
+
+        :param normalisation: Flux normalisation [GeV^-1 cm^-2 s^-1 sr^-1] or [GeV^-1 cm^-2 s^-1] for point sources.
+        :param norm_energy: Energy at which the spectrum is normalised [GeV].
+        :param index1: Spectral index of the power law.
+        :param cutoff_energy: Cutoff energy [GeV].
+        :param index2: Spectral index in the exponential function. 0 < index2 < 1.
+        :param lower_energy: Lower energy bound [GeV].
+        :param upper_energy: Upper energy bound [GeV].
+        """
+
+        super().__init__()
+
+        self._normalisation = normalisation
+
+        self._norm_energy = norm_energy
+
+        self._index1 = index1
+
+        self._cutoff_energy = cutoff_energy
+
+        self._index2 = index2
+
+        self._lower_energy = lower_energy
+
+        self._upper_energy = upper_energy
+
+        self.power_law = BoundedPowerLawSubexpCutoff(self._index1, self._cutoff_energy, self._index2, self._lower_energy, self._upper_energy)
+
+    def spectrum(self, energy):
+        """
+        dN/dEdAdt or dN/dEdAdtdO.
+        """
+        output = self._normalisation * (energy/self._norm_energy)**(-self._index1) * np.exp(-1 * (energy/self._cutoff_energy)**self._index2)
+
+        if isinstance(energy, np.ndarray):
+            idx = np.logical_or(energy < self._lower_energy, energy > self._upper_energy)
+            output[idx] = np.zeros(len(output[idx]))
+            return output
+        
+        else:
+            if energy < self._lower_energy or energy > self._upper_energy:
+                return 0.0
+            else:
+                return output
+            
+    def integrated_spectrum(self, lower_energy_bound, upper_energy_bound):
+        """
+        Integrates the spectrum with respect to E over a finite energy interval.
+        :param lower_energy_bound: in GeV
+        :param upper_energy_bound: in GeV
+        """
+        norm = self._normalisation
+        E0 = self._norm_energy
+        Ecut = self._cutoff_energy
+        gamma = self._index1
+        lambda_ = self._index2
+
+        E1 = lower_energy_bound
+        E2 = upper_energy_bound
+        incGamma = np.frompyfunc(mp.gammainc, 3, 1)
+
+        if isinstance(lower_energy_bound, np.ndarray) or isinstance(upper_energy_bound, np.ndarray):
+            return norm/E0**(-gamma) * Ecut**(1-gamma)/lambda_ * incGamma((1-gamma)/lambda_, (E1/Ecut)**lambda_, (E2/Ecut)**lambda_).astype('float64')
+        
+        else:
+            return norm/E0**(-gamma) * Ecut**(1-gamma)/lambda_ * float(incGamma((1-gamma)/lambda_, (E1/Ecut)**lambda_, (E2/Ecut)**lambda_))
+        
+    def redshift_factor(self, z: float):
+        return 1.0

icecube_tools/source/power_law.py

@@ -248,7 +248,7 @@ def cdf(self, x):
             if x < self.xmin:
                 return 0.0
             elif x > self.xmax:
-                return 0.0
+                return 1.0
             else:
                 return val
 
@@ -278,3 +278,73 @@ def samples(self, nsamples):
         """
         u = np.random.uniform(0, 1, size=nsamples)
         return self.inv_cdf(u)
+
+
+class BoundedPowerLawSubexpCutoff:
+    """
+    Definition of a bounded power law with a subexponential cutoff.
+    """
+
+    def __init__(self, gamma, xcut, lambda_, xmin, xmax):
+        """
+        Definition of a bounded power law with a subexponential cutoff.
+        """
+        self.gamma = gamma
+        self.xcut = xcut
+        self.lambda_ = lambda_
+        self.xmin = xmin
+        self.xmax = xmax
+
+        # Calculate the normalisation.
+        self.norm = 1. / float(
+            self.xcut ** (1 - self.gamma) / self.lambda_ 
+            * mp.gammainc((1 - self.gamma) / self.lambda_, 
+                          (self.xmin / self.xcut) ** self.lambda_, 
+                          (self.xmax/self.xcut) ** self.lambda_))
+        
+    def pdf(self, x):
+        """
+        Evaluates the probability density function at x.
+        """
+        val = self.norm * x ** (-self.gamma) * np.exp(-1 * (x / self.xcut) ** self.lambda_)
+
+        if not isinstance(x, np.ndarray):
+            if x < self.xmin or x > self.xmax:
+                return 0.0
+            else:
+                return val
+            
+        else:
+            idx = np.logical_or(x < self.xmin, x > self.xmax)
+            val[idx] = np.zeros(len(val[idx]))
+            return val
+        
+    def cdf(self, x):
+        """
+        Evaluated the cumulative distribution function at x.
+        """
+        incGamma = np.frompyfunc(mp.gammainc, 3, 1)
+        val = self.norm * self.xcut ** (1 - self.gamma) / self.lambda_ * incGamma((1 - self.gamma) / self.lambda_, 
+                                                                                  (self.xmin / self.xcut) ** self.lambda_, 
+                                                                                  (x / self.xcut) ** self.lambda_).astype('float64')
+        
+        if not isinstance(x, np.ndarray):
+            if x < self.xmin:
+                return 0.0
+            elif x > self.xmax:
+                return 1.0
+            else: 
+                return val
+            
+        else:
+            idx = x < self.xmin
+            val[idx] = np.zeros(len(val[idx]))
+            idx = x > self.xmax
+            val[idx] = np.ones(len(val[idx]))
+            return val
+        
+    def inv_cdf(self, x):
+        raise NotImplementedError
+    
+    def samples(self, nsamples):
+        raise NotImplementedError

tests/test_simulation.py

@@ -6,7 +6,7 @@
 from icecube_tools.detector.angular_resolution import AngularResolution
 from icecube_tools.detector.r2021 import R2021IRF
 from icecube_tools.detector.detector import IceCube
-from icecube_tools.source.flux_model import PowerLawFlux, BrokenPowerLawFlux, PowerLawExpCutoffFlux
+from icecube_tools.source.flux_model import PowerLawFlux, BrokenPowerLawFlux, PowerLawExpCutoffFlux, PowerLawSubexpCutoffFlux
 from icecube_tools.source.source_model import DiffuseSource, PointSource
 from icecube_tools.neutrino_calculator import NeutrinoCalculator, PhiSolver
 from icecube_tools.simulator import Simulator, TimeDependentSimulator
@@ -33,12 +33,18 @@
 point_source_plec = PointSource(plec_flux, z=0.0)
 diffuse_source_plec = DiffuseSource(plec_flux)
 
-sources = [point_source, diffuse_source, point_source_bpl, diffuse_source_bpl, point_source_plec, diffuse_source_plec]
+pl_subexp_cutoff_params = (1e-13, 1e3, 2., 1e3, 0.5, 1e4, 1e8)
+pl_subexp_cutoff_flux = PowerLawSubexpCutoffFlux(*pl_subexp_cutoff_params)
+point_source_pl_subexp_cutoff = PointSource(pl_subexp_cutoff_flux, z=0.0)
+diffuse_source_pl_subexp_cutoff = DiffuseSource(pl_subexp_cutoff_flux)
+
+sources = [point_source, diffuse_source, point_source_bpl, diffuse_source_bpl, point_source_plec, 
+           diffuse_source_plec, point_source_pl_subexp_cutoff, diffuse_source_pl_subexp_cutoff]
 
 
 def test_nu_calc():
 
-    nu_calc = NeutrinoCalculator([point_source, point_source_bpl, point_source_plec], aeff)
+    nu_calc = NeutrinoCalculator([point_source, point_source_bpl, point_source_plec, point_source_pl_subexp_cutoff], aeff)
 
     Nnu = nu_calc(
         time=1,

tests/test_source_interface.py

@@ -1,15 +1,16 @@
 import numpy as np
 
-from icecube_tools.source.flux_model import PowerLawFlux, BrokenPowerLawFlux, PowerLawExpCutoffFlux
+from icecube_tools.source.flux_model import PowerLawFlux, BrokenPowerLawFlux, PowerLawExpCutoffFlux, PowerLawSubexpCutoffFlux
 from icecube_tools.source.source_model import PointSource, DiffuseSource
 
 
 pl_params = (1e-18, 1e5, 2.2, 1e4, 1e8)
 bpl_params = (1e-18, 1e5, 2.2, 3.0, 1e4, 1e8)
 plec_params = (1e-13, 1e3, -2., 1e3, 1e4, 1e8)
+pl_subexp_cutoff_params = (1e-13, 1e3, 2., 1e3, 0.5, 1e4, 1e8)
 
-flux_models = [PowerLawFlux, BrokenPowerLawFlux, PowerLawExpCutoffFlux]
-params_list = [pl_params, bpl_params, plec_params]
+flux_models = [PowerLawFlux, BrokenPowerLawFlux, PowerLawExpCutoffFlux, PowerLawSubexpCutoffFlux]
+params_list = [pl_params, bpl_params, plec_params, pl_subexp_cutoff_params]
 
 
 def test_flux_models():
@@ -33,7 +34,7 @@ def test_flux_models():
 
             assert int_low > int_high
 
-        elif flux_model == PowerLawExpCutoffFlux:
+        elif flux_model == PowerLawExpCutoffFlux or flux_model == PowerLawSubexpCutoffFlux:
             integral = flux.integrated_spectrum(1e4, 1e6)
 
             assert integral > 0
@@ -43,11 +44,12 @@ def test_flux_models():
             assert flux.total_flux_density() > 0
 
         # Check sampling
-        samples = flux.sample(1000)
+        if not flux_model == PowerLawSubexpCutoffFlux:
+            samples = flux.sample(1000)
 
-        assert np.all(samples >= 1e4)
+            assert np.all(samples >= 1e4)
 
-        assert np.all(samples <= 1e8)
+            assert np.all(samples <= 1e8)
 
 
 def test_source_definition():