fixes from PR comments

agnpy/spectra/spectra.py

@@ -199,7 +199,7 @@ def plot(self, gamma=None, gamma_power=0, ax=None, **kwargs):
 
         return ax
 
-    def evaluate_time(self, time, energy_loss_function, subintervals_count=1):
+    def evaluate_time(self, time, energy_loss_function, subintervals_count=10):
         """Performs the time evaluation of energy change for this particle distribution.
 
         Parameters
@@ -226,11 +226,11 @@ def gamma_recalculated_after_loss(gamma):
             new_energy = old_energy - energy_loss
             if np.any(new_energy < 0):
                 raise ValueError(
-                    "Energy loss formula returned value higher then original energy. Use the shorter time ranges.")
+                    "Energy loss formula returned value higher then original energy. Use shorter time ranges.")
             new_gamma = (new_energy / mec2).value
             return new_gamma
 
-        gammas = self.gamma_values()
+        gammas = np.logspace(np.log10(self.gamma_min), np.log10(self.gamma_max), 200)
         n_array = self.__call__(gammas)
 
         # for each gamma point create a narrow bin, calculate the energy loss for start and end of the bin,
@@ -243,12 +243,8 @@ def gamma_recalculated_after_loss(gamma):
             narrowed_bins = bins_end_recalc - gammas
             if np.any(narrowed_bins <= 0):
                 raise ValueError(
-                    "Energy loss formula returned too big value - use shorter time intervals")
+                    "Energy loss formula returned too big value. Use shorter time ranges.")
             density_increase = bins_width / narrowed_bins
-            if np.any(density_increase < 1):
-                # not possible for simple scenarios because higher gammas should lose more energy,
-                # but maybe for more complex scenarios it could happen? then we should revise this assertion
-                raise AssertionError("Unexpected situation, bin has been widened instead of narrowed down")
             n_array = n_array * density_increase
 
         return InterpolatedDistribution(gammas, n_array)
@@ -305,9 +301,6 @@ def evaluate(gamma, k, p, gamma_min, gamma_max):
     def __call__(self, gamma):
         return self.evaluate(gamma, self.k, self.p, self.gamma_min, self.gamma_max)
 
-    def gamma_values(self):
-        return np.logspace(np.log10(self.gamma_min), np.log10(self.gamma_max), 200)
-
     @staticmethod
     def evaluate_SSA_integrand(gamma, k, p, gamma_min, gamma_max):
         r"""Analytical integrand for the synchrotron self-absorption:
@@ -787,7 +780,7 @@ class InterpolatedDistribution(ParticleDistribution):
         function to be used for integration, default is :class:`~numpy.trapz`
     """
 
-    def __init__(self, gamma, n, norm=1, mass=m_e, integrator=np.trapz, gamma_min=None, gamma_max=None):
+    def __init__(self, gamma, n, norm=1, mass=m_e, integrator=np.trapz):
         super().__init__(mass, integrator, "InterpolatedDistribution")
         if n.unit != u.Unit("cm-3"):
             raise ValueError(
@@ -799,9 +792,9 @@ def __init__(self, gamma, n, norm=1, mass=m_e, integrator=np.trapz, gamma_min=No
         # scaling parameter
         self.norm = norm
         # perform the interpolation
-        self.log10_interp = self.log10_interpolation(gamma, n, gamma_min, gamma_max)
+        self.log10_interp = self.log10_interpolation(gamma, n)
 
-    def log10_interpolation(self, gamma, n, gamma_min=None, gamma_max=None):
+    def log10_interpolation(self, gamma, n):
         """Returns the function interpolating in log10 the particle spectrum as
         a function of the Lorentz factor.
         TODO: make possible to pass arguments to CubicSpline.
@@ -815,8 +808,8 @@ def log10_interpolation(self, gamma, n, gamma_min=None, gamma_max=None):
 
         # min and max lorentz factor are now the first gamma values for which
         # the input distribution is not null
-        self.gamma_min = gamma_min if gamma_min else np.min(_gamma)
-        self.gamma_max = gamma_max if gamma_max else np.max(_gamma)
+        self.gamma_min = np.min(_gamma)
+        self.gamma_max = np.max(_gamma)
 
         return interpolator
 

agnpy/synchrotron/synchrotron.py

@@ -96,11 +96,12 @@ def electron_energy_loss_per_time(self, gamma):
         # In case of constant B, the part in brackets is fixed, so as an improvement we could calculate it once
         # and cache, and later we could use the formula (fixed * gamma^2)
         fixed = self._electron_energy_loss_formula_prefix()
-        return (fixed * gamma ** 2).to("erg/s")
+        return fixed * gamma ** 2
 
     def _electron_energy_loss_formula_prefix(self):
+        # using SI units here because of the ambiguity of the different CGS systems in terms of mu0 definition
         magn_energy_dens = (self.blob.B.to("T") ** 2) / (2 * mu0)
-        return (4 / 3) * sigma_T * c * magn_energy_dens
+        return ((4 / 3) * sigma_T * c * magn_energy_dens).to("erg/s")
 
     @staticmethod
     def evaluate_tau_ssa(