Skip to content

API Reference

crflux.models

:mod:models --- models of the cosmic ray flux at the top of Earth's atmosphere

This module is a collection of models of the high-energy primary cosmic ray flux, found in the literature of the last decades. The base class :class:PrimaryFlux has various methods, which can be employed in lepton flux calculations, cosmic ray physics, radiation physics etc.

The numbering scheme for nuclei is adapted from the Cosmic Ray Air-Shower Monte Carlo CORSIKA <https://web.ikp.kit.edu/corsika/>_. Protons have the ID 14. The composite ID for nuclei follows the formula :math:ID=100 \cdot A + Z, where :math:A is the mass number and :math:Z the charge. With this scheme one can easily obtain charge and mass from the ID and vice-versa (see :func:PrimaryFlux.Z_A).

The physics underlying each of the models can be found following the references in this documentation.

.. note::

As always, if you use one of the models for your work, please cite the corresponding
publication!

Example::

import crflux.models as mods
import numpy as np

model = mods.HillasGaisser2012("H3a")
E = np.logspace(1, 11, 100)
pfrac, p, n = model.p_and_n_flux(E)

# With geomagnetic cutoff of 7 GV
model_cut = mods.HillasGaisser2012("H3a", geomagnetic_cutoff=7.0)
pfrac_cut, p_cut, n_cut = model_cut.p_and_n_flux(E)

PrimaryFlux

Base class for primary cosmic ray flux models.

Parameters:

Name Type Description Default
geomagnetic_cutoff float

rigidity cutoff in GV. If set, the flux is zeroed below the corresponding energy for each nucleus. For a nucleus with charge Z, the energy cutoff is E_cut = Z * R_cut.

 None
Source code in crflux/models.py 
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
class PrimaryFlux(metaclass=ABCMeta):
    """Base class for primary cosmic ray flux models.

    Args:
      geomagnetic_cutoff (float, optional): rigidity cutoff in GV. If set,
        the flux is zeroed below the corresponding energy for each nucleus.
        For a nucleus with charge Z, the energy cutoff is E_cut = Z * R_cut.
    """

    def __init__(self, geomagnetic_cutoff=None):
        self.nucleus_ids = None
        self.geomagnetic_cutoff = geomagnetic_cutoff

    @abstractmethod
    def _nucleus_flux(self, corsika_id, E):
        """Returns the flux of nuclei corresponding to
        the ``corsika_id`` at energy ``E``.

        Args:
          corsika_id (int): see :mod:`crflux` for description.
          E (numpy.ndarray): laboratory energy of nucleus in GeV
        Returns:
          (numpy.ndarray): flux of single nucleus type :math:`\\Phi_{nucleus}`
          in :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
        """
        raise NotImplementedError(
            self.__class__.__name__ + "::_nucleus_flux(): Base class method called."
        )

    def nucleus_flux(self, corsika_id, E):
        """Returns the flux of nuclei corresponding to
        the ``corsika_id`` at energy ``E``, with optional geomagnetic cutoff.

        Args:
          corsika_id (int): see :mod:`crflux` for description.
          E (float or numpy.ndarray): laboratory energy of nucleus in GeV
        Returns:
          (numpy.ndarray): flux of single nucleus type :math:`\\Phi_{nucleus}`
          in :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
        """
        E = np.atleast_1d(np.asarray(E, dtype=float))
        flux = self._nucleus_flux(corsika_id, E)
        if self.geomagnetic_cutoff is not None:
            Z, A = self.Z_A(corsika_id)
            E_cut = Z * self.geomagnetic_cutoff
            flux = np.where(E >= E_cut, flux, 0.0)
        return flux

    def total_flux(self, E):
        """Returns total flux of nuclei, the "all-particle-flux".

        Args:
          E (float or numpy.ndarray): laboratory energy of particles in GeV
        Returns:
          (numpy.ndarray): particle flux in :math:`\\Phi_{particles}` in
          :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
        """
        return sum(self.nucleus_flux(corsika_id, E) for corsika_id in self.nucleus_ids)

    def tot_nucleon_flux(self, E):
        """Returns total flux of nucleons, the "all-nucleon-flux".

        Args:
          E (float or numpy.ndarray): laboratory energy of nucleons in GeV
        Returns:
          (numpy.ndarray): nucleon flux :math:`\\Phi_{nucleons}` in
          :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
        """
        return sum(
            self.Z_A(corsika_id)[1] ** 2.0
            * self.nucleus_flux(corsika_id, E * self.Z_A(corsika_id)[1])
            for corsika_id in self.nucleus_ids
        )

    def nucleon_gamma(self, E, rel_delta=0.01):
        """Returns differential spectral index of all-nucleon-flux
        obtained from a numerical derivative.

        Args:
          E (float or numpy.ndarray): laboratory energy of nucleons in GeV
          rel_delta (float): range of derivative relative to log10(E)
        Returns:
          (numpy.ndarray): spectral index :math:`\\gamma` of nucleons
        """
        delta = rel_delta * E
        fl = self.tot_nucleon_flux
        return np.log10(fl(E + delta) / fl(E - delta)) / np.log10(
            (E + delta) / (E - delta)
        )

    def nucleus_gamma(self, E, corsika_id, rel_delta=0.01):
        """Returns differential spectral index of nuclei
        obtained from a numerical derivative.

        Args:
          E (float or numpy.ndarray): laboratory energy of nuclei in GeV
          corsika_id (int): corsika id of nucleus/mass group
          rel_delta (float): range of derivative relative to log10(E)
        Returns:
          (numpy.ndarray): spectral index :math:`\\gamma` of nuclei
        """
        delta = rel_delta * E
        return np.log10(
            self.nucleus_flux(corsika_id, E + delta)
            / self.nucleus_flux(corsika_id, E - delta)
        ) / np.log10((E + delta) / (E - delta))

    def delta_0(self, E):
        """Returns proton excess.

        The proton excess is defined as
        :math:`\\delta_0 = \\frac{\\Phi_p - \\Phi_n}{\\Phi_p + \\Phi_n}`.

        Args:
          E (float or numpy.ndarray): laboratory energy of nucleons in GeV
        Returns:
          (numpy.ndarray): proton excess :math:`\\delta_0`
        """
        p_0 = 0.0
        n_0 = 0.0

        p_0 += self.nucleus_flux(14, E)

        p_0 += 2.0**2 * self.nucleus_flux(402, E * 4.0)
        n_0 += 2.0**2 * self.nucleus_flux(402, E * 4.0)

        p_0 += 6.0**2 * self.nucleus_flux(1206, E * 12.0)
        n_0 += 6.0**2 * self.nucleus_flux(1206, E * 12.0)

        p_0 += 14.0**2 * self.nucleus_flux(2814, E * 28.0)
        n_0 += 14.0**2 * self.nucleus_flux(2814, E * 28.0)

        a_fe = self.Z_A(5426)[1]
        p_0 += 26.0**2 * self.nucleus_flux(5426, E * a_fe)
        n_0 += 26.0**2 * self.nucleus_flux(5426, E * a_fe)

        return (p_0 - n_0) / (p_0 + n_0)

    def p_and_n_flux(self, E):
        """Returns tuple with proton fraction, proton flux and neutron flux.

        The proton fraction is defined as :math:`\\frac{\\Phi_p}{\\Phi_p + \\Phi_n}`.

        Args:
          E (float or numpy.ndarray): laboratory energy of nucleons in GeV
        Returns:
          (float,float,float): proton fraction, proton flux, neutron flux
        """
        za = self.Z_A

        p_flux = sum(
            za(corsika_id)[0]
            * za(corsika_id)[1]
            * self.nucleus_flux(corsika_id, E * za(corsika_id)[1])
            for corsika_id in self.nucleus_ids
        )

        n_flux = sum(
            (za(corsika_id)[1] - za(corsika_id)[0])
            * za(corsika_id)[1]
            * self.nucleus_flux(corsika_id, E * za(corsika_id)[1])
            for corsika_id in self.nucleus_ids
        )

        total = p_flux + n_flux
        p_frac = np.where(total > 0, p_flux / total, 0.0)
        return p_frac, p_flux, n_flux

    def lnA(self, E):
        """Returns mean logarithmic mass <ln A>.

        Args:
          E (float or numpy.ndarray): laboratory energy of particles in GeV
        Returns:
          (numpy.ndarray): mean (natural) logarithmic mass
        """
        sum_weight = 0.0

        for cid in self.nucleus_ids:
            if cid == 14:
                continue  # p has lnA = 0
            sum_weight += np.log(self.Z_A(cid)[1]) * self.nucleus_flux(cid, E)

        return sum_weight / self.total_flux(E)

    def _find_nearby_id(self, corsika_id, delta_A=3):
        """Looks in :attr:`self.params` for a nucleus with same ``corsika_id``
        and returns the corsika_id if these parameters exist. If not, it will
        look for nuclei of mass number +- delta_A around the requested nucleus
        and return its corsika_id if it exists.

        Args:
          corsika_id (int): corsika id of nucleus/mass group
          delta_A (int): maximum mass number difference to accept
        Returns:
          (int): corsika_id of requested or similar nucleus
        Raises:
          Exception: if no nucleus with mass number closer than delta_A can
            be found in parameters
        """
        if corsika_id in self.nucleus_ids:
            return corsika_id

        A_in = (corsika_id - corsika_id % 100) // 100
        closest_id = _get_closest(corsika_id, np.array(self.nucleus_ids))[1]
        A_close = (closest_id - closest_id % 100) // 100
        if np.abs(A_in - A_close) > delta_A:
            e = (
                "{0}::_find_nearby_id(): No similar nucleus found with "
                + "delta_A <= {1} for A_in = {2}. Closest is {3}."
            )
            raise Exception(e.format(self.__class__.__name__, delta_A, A_in, A_close))
        else:
            return closest_id

    def Z_A(self, corsika_id):
        """Returns charge :math:`Z` and mass number :math:`A` corresponding
        to ``corsika_id``

        Args:
          corsika_id (int): corsika id of nucleus/mass group
        Returns:
          (int,int): (Z, A) tuple
        """
        Z, A = 1, 1
        if corsika_id > 14:
            Z = corsika_id % 100
            A = (corsika_id - Z) // 100
        return Z, A

nucleus_flux(corsika_id, E)

Returns the flux of nuclei corresponding to the corsika_id at energy E, with optional geomagnetic cutoff.

Parameters:

Name Type Description Default
corsika_id int

see :mod:crflux for description.

 required
E float or ndarray

laboratory energy of nucleus in GeV

 required

Returns: (numpy.ndarray): flux of single nucleus type :math:\Phi_{nucleus} in :math:(\text{m}^2 \text{s sr GeV})^{-1}

Source code in crflux/models.py 
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
def nucleus_flux(self, corsika_id, E):
    """Returns the flux of nuclei corresponding to
    the ``corsika_id`` at energy ``E``, with optional geomagnetic cutoff.

    Args:
      corsika_id (int): see :mod:`crflux` for description.
      E (float or numpy.ndarray): laboratory energy of nucleus in GeV
    Returns:
      (numpy.ndarray): flux of single nucleus type :math:`\\Phi_{nucleus}`
      in :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
    """
    E = np.atleast_1d(np.asarray(E, dtype=float))
    flux = self._nucleus_flux(corsika_id, E)
    if self.geomagnetic_cutoff is not None:
        Z, A = self.Z_A(corsika_id)
        E_cut = Z * self.geomagnetic_cutoff
        flux = np.where(E >= E_cut, flux, 0.0)
    return flux

total_flux(E)

Returns total flux of nuclei, the "all-particle-flux".

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of particles in GeV

 required

Returns: (numpy.ndarray): particle flux in :math:\Phi_{particles} in :math:(\text{m}^2 \text{s sr GeV})^{-1}

Source code in crflux/models.py 
107
108
109
110
111
112
113
114
115
116
def total_flux(self, E):
    """Returns total flux of nuclei, the "all-particle-flux".

    Args:
      E (float or numpy.ndarray): laboratory energy of particles in GeV
    Returns:
      (numpy.ndarray): particle flux in :math:`\\Phi_{particles}` in
      :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
    """
    return sum(self.nucleus_flux(corsika_id, E) for corsika_id in self.nucleus_ids)

tot_nucleon_flux(E)

Returns total flux of nucleons, the "all-nucleon-flux".

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required

Returns: (numpy.ndarray): nucleon flux :math:\Phi_{nucleons} in :math:(\text{m}^2 \text{s sr GeV})^{-1}

Source code in crflux/models.py 
118
119
120
121
122
123
124
125
126
127
128
129
130
131
def tot_nucleon_flux(self, E):
    """Returns total flux of nucleons, the "all-nucleon-flux".

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
    Returns:
      (numpy.ndarray): nucleon flux :math:`\\Phi_{nucleons}` in
      :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
    """
    return sum(
        self.Z_A(corsika_id)[1] ** 2.0
        * self.nucleus_flux(corsika_id, E * self.Z_A(corsika_id)[1])
        for corsika_id in self.nucleus_ids
    )

nucleon_gamma(E, rel_delta=0.01)

Returns differential spectral index of all-nucleon-flux obtained from a numerical derivative.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required
rel_delta float

range of derivative relative to log10(E)

 0.01

Returns: (numpy.ndarray): spectral index :math:\gamma of nucleons

Source code in crflux/models.py 
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
def nucleon_gamma(self, E, rel_delta=0.01):
    """Returns differential spectral index of all-nucleon-flux
    obtained from a numerical derivative.

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
      rel_delta (float): range of derivative relative to log10(E)
    Returns:
      (numpy.ndarray): spectral index :math:`\\gamma` of nucleons
    """
    delta = rel_delta * E
    fl = self.tot_nucleon_flux
    return np.log10(fl(E + delta) / fl(E - delta)) / np.log10(
        (E + delta) / (E - delta)
    )

nucleus_gamma(E, corsika_id, rel_delta=0.01)

Returns differential spectral index of nuclei obtained from a numerical derivative.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nuclei in GeV

 required
corsika_id int

corsika id of nucleus/mass group

 required
rel_delta float

range of derivative relative to log10(E)

 0.01

Returns: (numpy.ndarray): spectral index :math:\gamma of nuclei

Source code in crflux/models.py 
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
def nucleus_gamma(self, E, corsika_id, rel_delta=0.01):
    """Returns differential spectral index of nuclei
    obtained from a numerical derivative.

    Args:
      E (float or numpy.ndarray): laboratory energy of nuclei in GeV
      corsika_id (int): corsika id of nucleus/mass group
      rel_delta (float): range of derivative relative to log10(E)
    Returns:
      (numpy.ndarray): spectral index :math:`\\gamma` of nuclei
    """
    delta = rel_delta * E
    return np.log10(
        self.nucleus_flux(corsika_id, E + delta)
        / self.nucleus_flux(corsika_id, E - delta)
    ) / np.log10((E + delta) / (E - delta))

delta_0(E)

Returns proton excess.

The proton excess is defined as :math:\delta_0 = \frac{\Phi_p - \Phi_n}{\Phi_p + \Phi_n}.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required

Returns: (numpy.ndarray): proton excess :math:\delta_0

Source code in crflux/models.py 
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
def delta_0(self, E):
    """Returns proton excess.

    The proton excess is defined as
    :math:`\\delta_0 = \\frac{\\Phi_p - \\Phi_n}{\\Phi_p + \\Phi_n}`.

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
    Returns:
      (numpy.ndarray): proton excess :math:`\\delta_0`
    """
    p_0 = 0.0
    n_0 = 0.0

    p_0 += self.nucleus_flux(14, E)

    p_0 += 2.0**2 * self.nucleus_flux(402, E * 4.0)
    n_0 += 2.0**2 * self.nucleus_flux(402, E * 4.0)

    p_0 += 6.0**2 * self.nucleus_flux(1206, E * 12.0)
    n_0 += 6.0**2 * self.nucleus_flux(1206, E * 12.0)

    p_0 += 14.0**2 * self.nucleus_flux(2814, E * 28.0)
    n_0 += 14.0**2 * self.nucleus_flux(2814, E * 28.0)

    a_fe = self.Z_A(5426)[1]
    p_0 += 26.0**2 * self.nucleus_flux(5426, E * a_fe)
    n_0 += 26.0**2 * self.nucleus_flux(5426, E * a_fe)

    return (p_0 - n_0) / (p_0 + n_0)

p_and_n_flux(E)

Returns tuple with proton fraction, proton flux and neutron flux.

The proton fraction is defined as :math:\frac{\Phi_p}{\Phi_p + \Phi_n}.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required

Returns: (float,float,float): proton fraction, proton flux, neutron flux

Source code in crflux/models.py 
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
def p_and_n_flux(self, E):
    """Returns tuple with proton fraction, proton flux and neutron flux.

    The proton fraction is defined as :math:`\\frac{\\Phi_p}{\\Phi_p + \\Phi_n}`.

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
    Returns:
      (float,float,float): proton fraction, proton flux, neutron flux
    """
    za = self.Z_A

    p_flux = sum(
        za(corsika_id)[0]
        * za(corsika_id)[1]
        * self.nucleus_flux(corsika_id, E * za(corsika_id)[1])
        for corsika_id in self.nucleus_ids
    )

    n_flux = sum(
        (za(corsika_id)[1] - za(corsika_id)[0])
        * za(corsika_id)[1]
        * self.nucleus_flux(corsika_id, E * za(corsika_id)[1])
        for corsika_id in self.nucleus_ids
    )

    total = p_flux + n_flux
    p_frac = np.where(total > 0, p_flux / total, 0.0)
    return p_frac, p_flux, n_flux

lnA(E)

Returns mean logarithmic mass .

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of particles in GeV

 required

Returns: (numpy.ndarray): mean (natural) logarithmic mass

Source code in crflux/models.py 
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
def lnA(self, E):
    """Returns mean logarithmic mass <ln A>.

    Args:
      E (float or numpy.ndarray): laboratory energy of particles in GeV
    Returns:
      (numpy.ndarray): mean (natural) logarithmic mass
    """
    sum_weight = 0.0

    for cid in self.nucleus_ids:
        if cid == 14:
            continue  # p has lnA = 0
        sum_weight += np.log(self.Z_A(cid)[1]) * self.nucleus_flux(cid, E)

    return sum_weight / self.total_flux(E)

Z_A(corsika_id)

Returns charge :math:Z and mass number :math:A corresponding to corsika_id

Parameters:

Name Type Description Default
corsika_id int

corsika id of nucleus/mass group

 required

Returns: (int,int): (Z, A) tuple

Source code in crflux/models.py 
274
275
276
277
278
279
280
281
282
283
284
285
286
287
def Z_A(self, corsika_id):
    """Returns charge :math:`Z` and mass number :math:`A` corresponding
    to ``corsika_id``

    Args:
      corsika_id (int): corsika id of nucleus/mass group
    Returns:
      (int,int): (Z, A) tuple
    """
    Z, A = 1, 1
    if corsika_id > 14:
        Z = corsika_id % 100
        A = (corsika_id - Z) // 100
    return Z, A

PolyGonato

Bases: PrimaryFlux

J. R. Hoerandel, Astroparticle Physics 19, 193 (2003).

Source code in crflux/models.py 
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
class PolyGonato(PrimaryFlux):
    """J. R. Hoerandel, Astroparticle Physics 19, 193 (2003)."""

    def __init__(self, constdelta=False, **kwargs):
        super().__init__(**kwargs)

        self.name = "poly-gonato"
        self.sname = "pg"
        self.constdelta = constdelta
        self.E_p = 4.51e6
        self.gamma_c = -4.68
        self.epsilon_c = 1.87
        self.delta_gamma = 2.10
        if self.constdelta:
            self.epsilon_c = 1.90
            self.E_p = 4.49e6

        self.params = {}

        self.params[14] = (8.73e-2, 2.71, 1)  # H
        self.params[402] = (5.71e-2, 2.64, 2)  # He
        self.params[1206] = (1.06e-2, 2.66, 6)  # C
        self.params[1407] = (2.35e-3, 2.72, 7)  # N
        self.params[1608] = (1.57e-2, 2.68, 8)  # O
        self.params[2412] = (8.01e-3, 2.64, 12)  # Mg
        self.params[2613] = (1.15e-3, 2.66, 13)  # Al
        self.params[2814] = (7.96e-3, 2.75, 14)  # Si
        self.params[5025] = (1.35e-3, 2.46, 25)  # Mn
        self.params[5426] = (2.04e-2, 2.59, 26)  # Fe
        self.params[5427] = (7.51e-5, 2.72, 27)  # Co

        self.nucleus_ids = list(self.params.keys())

    def _nucleus_flux(self, corsika_id, E):
        corsika_id = self._find_nearby_id(corsika_id)
        return self._polygonato_formula(corsika_id, E)

    def _polygonato_formula(self, corsika_id, E):
        p = self.params[corsika_id]
        gam = (self.gamma_c + p[1]) if self.constdelta else -self.delta_gamma

        return (
            p[0]
            / 1000.0
            * (E / 1000.0) ** (-p[1])
            * (1 + (E / p[2] / self.E_p) ** self.epsilon_c) ** (gam / self.epsilon_c)
        )

nucleus_flux(corsika_id, E)

Returns the flux of nuclei corresponding to the corsika_id at energy E, with optional geomagnetic cutoff.

Parameters:

Name Type Description Default
corsika_id int

see :mod:crflux for description.

 required
E float or ndarray

laboratory energy of nucleus in GeV

 required

Returns: (numpy.ndarray): flux of single nucleus type :math:\Phi_{nucleus} in :math:(\text{m}^2 \text{s sr GeV})^{-1}

Source code in crflux/models.py 
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
def nucleus_flux(self, corsika_id, E):
    """Returns the flux of nuclei corresponding to
    the ``corsika_id`` at energy ``E``, with optional geomagnetic cutoff.

    Args:
      corsika_id (int): see :mod:`crflux` for description.
      E (float or numpy.ndarray): laboratory energy of nucleus in GeV
    Returns:
      (numpy.ndarray): flux of single nucleus type :math:`\\Phi_{nucleus}`
      in :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
    """
    E = np.atleast_1d(np.asarray(E, dtype=float))
    flux = self._nucleus_flux(corsika_id, E)
    if self.geomagnetic_cutoff is not None:
        Z, A = self.Z_A(corsika_id)
        E_cut = Z * self.geomagnetic_cutoff
        flux = np.where(E >= E_cut, flux, 0.0)
    return flux

total_flux(E)

Returns total flux of nuclei, the "all-particle-flux".

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of particles in GeV

 required

Returns: (numpy.ndarray): particle flux in :math:\Phi_{particles} in :math:(\text{m}^2 \text{s sr GeV})^{-1}

Source code in crflux/models.py 
107
108
109
110
111
112
113
114
115
116
def total_flux(self, E):
    """Returns total flux of nuclei, the "all-particle-flux".

    Args:
      E (float or numpy.ndarray): laboratory energy of particles in GeV
    Returns:
      (numpy.ndarray): particle flux in :math:`\\Phi_{particles}` in
      :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
    """
    return sum(self.nucleus_flux(corsika_id, E) for corsika_id in self.nucleus_ids)

tot_nucleon_flux(E)

Returns total flux of nucleons, the "all-nucleon-flux".

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required

Returns: (numpy.ndarray): nucleon flux :math:\Phi_{nucleons} in :math:(\text{m}^2 \text{s sr GeV})^{-1}

Source code in crflux/models.py 
118
119
120
121
122
123
124
125
126
127
128
129
130
131
def tot_nucleon_flux(self, E):
    """Returns total flux of nucleons, the "all-nucleon-flux".

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
    Returns:
      (numpy.ndarray): nucleon flux :math:`\\Phi_{nucleons}` in
      :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
    """
    return sum(
        self.Z_A(corsika_id)[1] ** 2.0
        * self.nucleus_flux(corsika_id, E * self.Z_A(corsika_id)[1])
        for corsika_id in self.nucleus_ids
    )

nucleon_gamma(E, rel_delta=0.01)

Returns differential spectral index of all-nucleon-flux obtained from a numerical derivative.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required
rel_delta float

range of derivative relative to log10(E)

 0.01

Returns: (numpy.ndarray): spectral index :math:\gamma of nucleons

Source code in crflux/models.py 
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
def nucleon_gamma(self, E, rel_delta=0.01):
    """Returns differential spectral index of all-nucleon-flux
    obtained from a numerical derivative.

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
      rel_delta (float): range of derivative relative to log10(E)
    Returns:
      (numpy.ndarray): spectral index :math:`\\gamma` of nucleons
    """
    delta = rel_delta * E
    fl = self.tot_nucleon_flux
    return np.log10(fl(E + delta) / fl(E - delta)) / np.log10(
        (E + delta) / (E - delta)
    )

nucleus_gamma(E, corsika_id, rel_delta=0.01)

Returns differential spectral index of nuclei obtained from a numerical derivative.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nuclei in GeV

 required
corsika_id int

corsika id of nucleus/mass group

 required
rel_delta float

range of derivative relative to log10(E)

 0.01

Returns: (numpy.ndarray): spectral index :math:\gamma of nuclei

Source code in crflux/models.py 
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
def nucleus_gamma(self, E, corsika_id, rel_delta=0.01):
    """Returns differential spectral index of nuclei
    obtained from a numerical derivative.

    Args:
      E (float or numpy.ndarray): laboratory energy of nuclei in GeV
      corsika_id (int): corsika id of nucleus/mass group
      rel_delta (float): range of derivative relative to log10(E)
    Returns:
      (numpy.ndarray): spectral index :math:`\\gamma` of nuclei
    """
    delta = rel_delta * E
    return np.log10(
        self.nucleus_flux(corsika_id, E + delta)
        / self.nucleus_flux(corsika_id, E - delta)
    ) / np.log10((E + delta) / (E - delta))

delta_0(E)

Returns proton excess.

The proton excess is defined as :math:\delta_0 = \frac{\Phi_p - \Phi_n}{\Phi_p + \Phi_n}.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required

Returns: (numpy.ndarray): proton excess :math:\delta_0

Source code in crflux/models.py 
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
def delta_0(self, E):
    """Returns proton excess.

    The proton excess is defined as
    :math:`\\delta_0 = \\frac{\\Phi_p - \\Phi_n}{\\Phi_p + \\Phi_n}`.

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
    Returns:
      (numpy.ndarray): proton excess :math:`\\delta_0`
    """
    p_0 = 0.0
    n_0 = 0.0

    p_0 += self.nucleus_flux(14, E)

    p_0 += 2.0**2 * self.nucleus_flux(402, E * 4.0)
    n_0 += 2.0**2 * self.nucleus_flux(402, E * 4.0)

    p_0 += 6.0**2 * self.nucleus_flux(1206, E * 12.0)
    n_0 += 6.0**2 * self.nucleus_flux(1206, E * 12.0)

    p_0 += 14.0**2 * self.nucleus_flux(2814, E * 28.0)
    n_0 += 14.0**2 * self.nucleus_flux(2814, E * 28.0)

    a_fe = self.Z_A(5426)[1]
    p_0 += 26.0**2 * self.nucleus_flux(5426, E * a_fe)
    n_0 += 26.0**2 * self.nucleus_flux(5426, E * a_fe)

    return (p_0 - n_0) / (p_0 + n_0)

p_and_n_flux(E)

Returns tuple with proton fraction, proton flux and neutron flux.

The proton fraction is defined as :math:\frac{\Phi_p}{\Phi_p + \Phi_n}.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required

Returns: (float,float,float): proton fraction, proton flux, neutron flux

Source code in crflux/models.py 
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
def p_and_n_flux(self, E):
    """Returns tuple with proton fraction, proton flux and neutron flux.

    The proton fraction is defined as :math:`\\frac{\\Phi_p}{\\Phi_p + \\Phi_n}`.

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
    Returns:
      (float,float,float): proton fraction, proton flux, neutron flux
    """
    za = self.Z_A

    p_flux = sum(
        za(corsika_id)[0]
        * za(corsika_id)[1]
        * self.nucleus_flux(corsika_id, E * za(corsika_id)[1])
        for corsika_id in self.nucleus_ids
    )

    n_flux = sum(
        (za(corsika_id)[1] - za(corsika_id)[0])
        * za(corsika_id)[1]
        * self.nucleus_flux(corsika_id, E * za(corsika_id)[1])
        for corsika_id in self.nucleus_ids
    )

    total = p_flux + n_flux
    p_frac = np.where(total > 0, p_flux / total, 0.0)
    return p_frac, p_flux, n_flux

lnA(E)

Returns mean logarithmic mass .

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of particles in GeV

 required

Returns: (numpy.ndarray): mean (natural) logarithmic mass

Source code in crflux/models.py 
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
def lnA(self, E):
    """Returns mean logarithmic mass <ln A>.

    Args:
      E (float or numpy.ndarray): laboratory energy of particles in GeV
    Returns:
      (numpy.ndarray): mean (natural) logarithmic mass
    """
    sum_weight = 0.0

    for cid in self.nucleus_ids:
        if cid == 14:
            continue  # p has lnA = 0
        sum_weight += np.log(self.Z_A(cid)[1]) * self.nucleus_flux(cid, E)

    return sum_weight / self.total_flux(E)

Z_A(corsika_id)

Returns charge :math:Z and mass number :math:A corresponding to corsika_id

Parameters:

Name Type Description Default
corsika_id int

corsika id of nucleus/mass group

 required

Returns: (int,int): (Z, A) tuple

Source code in crflux/models.py 
274
275
276
277
278
279
280
281
282
283
284
285
286
287
def Z_A(self, corsika_id):
    """Returns charge :math:`Z` and mass number :math:`A` corresponding
    to ``corsika_id``

    Args:
      corsika_id (int): corsika id of nucleus/mass group
    Returns:
      (int,int): (Z, A) tuple
    """
    Z, A = 1, 1
    if corsika_id > 14:
        Z = corsika_id % 100
        A = (corsika_id - Z) // 100
    return Z, A

HillasGaisser2012

Bases: PrimaryFlux

Gaisser, T.K., Astroparticle Physics 35, 801 (2012).

Model is based on Hillas ideas and eye-ball fits by T.K. Gaisser. H3a is a 3-'peters cycle' 5 mass group model with mixed composition above the ankle. H4a has protons only in the 3. component.

Parameters:

Name Type Description Default
model str

can be either H3a or H4a.

 'H4a'
Source code in crflux/models.py 
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
class HillasGaisser2012(PrimaryFlux):
    """Gaisser, T.K., Astroparticle Physics 35, 801 (2012).

    Model is based on Hillas ideas and eye-ball fits by T.K. Gaisser.
    H3a is a 3-'peters cycle' 5 mass group model with mixed
    composition above the ankle. H4a has protons only in the
    3. component.

    Args:
      model (str): can be either H3a or H4a.
    """

    def __init__(self, model="H4a", **kwargs):
        super().__init__(**kwargs)

        self.name = "Hillas-Gaisser (" + model + ")"
        self.sname = model
        self.model = model
        self.params = {}
        self.rid_cutoff = {}

        mass_comp = [14, 402, 1206, 2814, 5426]
        for mcomp in mass_comp:
            self.params[mcomp] = {}

        self.rid_cutoff[1] = 4e6
        self.rid_cutoff[2] = 30e6
        self.rid_cutoff[3] = 2e9
        self.params[14][1] = (7860, 1.66, 1)  # H
        self.params[402][1] = (3550, 1.58, 2)  # He
        self.params[1206][1] = (2200, 1.63, 6)  # CNO
        self.params[2814][1] = (1430, 1.67, 14)  # MgAlSi
        self.params[5426][1] = (2120, 1.63, 26)  # Fe

        self.params[14][2] = (20, 1.4, 1)  # H
        self.params[402][2] = (20, 1.4, 2)  # He
        self.params[1206][2] = (13.4, 1.4, 6)  # CNO
        self.params[2814][2] = (13.4, 1.4, 14)  # MgAlSi
        self.params[5426][2] = (13.4, 1.4, 26)  # Fe

        if self.model == "H3a":
            self.rid_cutoff[3] = 2e9
            self.params[14][3] = (1.7, 1.4, 1)  # H
            self.params[402][3] = (1.7, 1.4, 2)  # He
            self.params[1206][3] = (1.14, 1.4, 6)  # CNO
            self.params[2814][3] = (1.14, 1.4, 14)  # MgAlSi
            self.params[5426][3] = (1.14, 1.4, 26)  # Fe
        elif self.model == "H4a":
            self.rid_cutoff[3] = 60e9
            self.params[14][3] = (200.0, 1.6, 1)  # H
            self.params[402][3] = (0, 1.4, 2)  # He
            self.params[1206][3] = (0, 1.4, 6)  # CNO
            self.params[2814][3] = (0, 1.4, 14)  # MgAlSi
            self.params[5426][3] = (0, 1.4, 26)  # Fe
        else:
            raise Exception("HillasGaisser2012(): Unknown model version requested.")

        self.nucleus_ids = list(self.params.keys())

    def _nucleus_flux(self, corsika_id, E):
        corsika_id = self._find_nearby_id(corsika_id)

        flux = 0.0
        for i in range(1, 4):
            p = self.params[corsika_id][i]
            flux += p[0] * E ** (-p[1] - 1.0) * np.exp(-E / p[2] / self.rid_cutoff[i])
        return flux

nucleus_flux(corsika_id, E)

Returns the flux of nuclei corresponding to the corsika_id at energy E, with optional geomagnetic cutoff.

Parameters:

Name Type Description Default
corsika_id int

see :mod:crflux for description.

 required
E float or ndarray

laboratory energy of nucleus in GeV

 required

Returns: (numpy.ndarray): flux of single nucleus type :math:\Phi_{nucleus} in :math:(\text{m}^2 \text{s sr GeV})^{-1}

Source code in crflux/models.py 
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
def nucleus_flux(self, corsika_id, E):
    """Returns the flux of nuclei corresponding to
    the ``corsika_id`` at energy ``E``, with optional geomagnetic cutoff.

    Args:
      corsika_id (int): see :mod:`crflux` for description.
      E (float or numpy.ndarray): laboratory energy of nucleus in GeV
    Returns:
      (numpy.ndarray): flux of single nucleus type :math:`\\Phi_{nucleus}`
      in :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
    """
    E = np.atleast_1d(np.asarray(E, dtype=float))
    flux = self._nucleus_flux(corsika_id, E)
    if self.geomagnetic_cutoff is not None:
        Z, A = self.Z_A(corsika_id)
        E_cut = Z * self.geomagnetic_cutoff
        flux = np.where(E >= E_cut, flux, 0.0)
    return flux

total_flux(E)

Returns total flux of nuclei, the "all-particle-flux".

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of particles in GeV

 required

Returns: (numpy.ndarray): particle flux in :math:\Phi_{particles} in :math:(\text{m}^2 \text{s sr GeV})^{-1}

Source code in crflux/models.py 
107
108
109
110
111
112
113
114
115
116
def total_flux(self, E):
    """Returns total flux of nuclei, the "all-particle-flux".

    Args:
      E (float or numpy.ndarray): laboratory energy of particles in GeV
    Returns:
      (numpy.ndarray): particle flux in :math:`\\Phi_{particles}` in
      :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
    """
    return sum(self.nucleus_flux(corsika_id, E) for corsika_id in self.nucleus_ids)

tot_nucleon_flux(E)

Returns total flux of nucleons, the "all-nucleon-flux".

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required

Returns: (numpy.ndarray): nucleon flux :math:\Phi_{nucleons} in :math:(\text{m}^2 \text{s sr GeV})^{-1}

Source code in crflux/models.py 
118
119
120
121
122
123
124
125
126
127
128
129
130
131
def tot_nucleon_flux(self, E):
    """Returns total flux of nucleons, the "all-nucleon-flux".

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
    Returns:
      (numpy.ndarray): nucleon flux :math:`\\Phi_{nucleons}` in
      :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
    """
    return sum(
        self.Z_A(corsika_id)[1] ** 2.0
        * self.nucleus_flux(corsika_id, E * self.Z_A(corsika_id)[1])
        for corsika_id in self.nucleus_ids
    )

nucleon_gamma(E, rel_delta=0.01)

Returns differential spectral index of all-nucleon-flux obtained from a numerical derivative.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required
rel_delta float

range of derivative relative to log10(E)

 0.01

Returns: (numpy.ndarray): spectral index :math:\gamma of nucleons

Source code in crflux/models.py 
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
def nucleon_gamma(self, E, rel_delta=0.01):
    """Returns differential spectral index of all-nucleon-flux
    obtained from a numerical derivative.

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
      rel_delta (float): range of derivative relative to log10(E)
    Returns:
      (numpy.ndarray): spectral index :math:`\\gamma` of nucleons
    """
    delta = rel_delta * E
    fl = self.tot_nucleon_flux
    return np.log10(fl(E + delta) / fl(E - delta)) / np.log10(
        (E + delta) / (E - delta)
    )

nucleus_gamma(E, corsika_id, rel_delta=0.01)

Returns differential spectral index of nuclei obtained from a numerical derivative.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nuclei in GeV

 required
corsika_id int

corsika id of nucleus/mass group

 required
rel_delta float

range of derivative relative to log10(E)

 0.01

Returns: (numpy.ndarray): spectral index :math:\gamma of nuclei

Source code in crflux/models.py 
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
def nucleus_gamma(self, E, corsika_id, rel_delta=0.01):
    """Returns differential spectral index of nuclei
    obtained from a numerical derivative.

    Args:
      E (float or numpy.ndarray): laboratory energy of nuclei in GeV
      corsika_id (int): corsika id of nucleus/mass group
      rel_delta (float): range of derivative relative to log10(E)
    Returns:
      (numpy.ndarray): spectral index :math:`\\gamma` of nuclei
    """
    delta = rel_delta * E
    return np.log10(
        self.nucleus_flux(corsika_id, E + delta)
        / self.nucleus_flux(corsika_id, E - delta)
    ) / np.log10((E + delta) / (E - delta))

delta_0(E)

Returns proton excess.

The proton excess is defined as :math:\delta_0 = \frac{\Phi_p - \Phi_n}{\Phi_p + \Phi_n}.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required

Returns: (numpy.ndarray): proton excess :math:\delta_0

Source code in crflux/models.py 
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
def delta_0(self, E):
    """Returns proton excess.

    The proton excess is defined as
    :math:`\\delta_0 = \\frac{\\Phi_p - \\Phi_n}{\\Phi_p + \\Phi_n}`.

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
    Returns:
      (numpy.ndarray): proton excess :math:`\\delta_0`
    """
    p_0 = 0.0
    n_0 = 0.0

    p_0 += self.nucleus_flux(14, E)

    p_0 += 2.0**2 * self.nucleus_flux(402, E * 4.0)
    n_0 += 2.0**2 * self.nucleus_flux(402, E * 4.0)

    p_0 += 6.0**2 * self.nucleus_flux(1206, E * 12.0)
    n_0 += 6.0**2 * self.nucleus_flux(1206, E * 12.0)

    p_0 += 14.0**2 * self.nucleus_flux(2814, E * 28.0)
    n_0 += 14.0**2 * self.nucleus_flux(2814, E * 28.0)

    a_fe = self.Z_A(5426)[1]
    p_0 += 26.0**2 * self.nucleus_flux(5426, E * a_fe)
    n_0 += 26.0**2 * self.nucleus_flux(5426, E * a_fe)

    return (p_0 - n_0) / (p_0 + n_0)

p_and_n_flux(E)

Returns tuple with proton fraction, proton flux and neutron flux.

The proton fraction is defined as :math:\frac{\Phi_p}{\Phi_p + \Phi_n}.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required

Returns: (float,float,float): proton fraction, proton flux, neutron flux

Source code in crflux/models.py 
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
def p_and_n_flux(self, E):
    """Returns tuple with proton fraction, proton flux and neutron flux.

    The proton fraction is defined as :math:`\\frac{\\Phi_p}{\\Phi_p + \\Phi_n}`.

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
    Returns:
      (float,float,float): proton fraction, proton flux, neutron flux
    """
    za = self.Z_A

    p_flux = sum(
        za(corsika_id)[0]
        * za(corsika_id)[1]
        * self.nucleus_flux(corsika_id, E * za(corsika_id)[1])
        for corsika_id in self.nucleus_ids
    )

    n_flux = sum(
        (za(corsika_id)[1] - za(corsika_id)[0])
        * za(corsika_id)[1]
        * self.nucleus_flux(corsika_id, E * za(corsika_id)[1])
        for corsika_id in self.nucleus_ids
    )

    total = p_flux + n_flux
    p_frac = np.where(total > 0, p_flux / total, 0.0)
    return p_frac, p_flux, n_flux

lnA(E)

Returns mean logarithmic mass .

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of particles in GeV

 required

Returns: (numpy.ndarray): mean (natural) logarithmic mass

Source code in crflux/models.py 
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
def lnA(self, E):
    """Returns mean logarithmic mass <ln A>.

    Args:
      E (float or numpy.ndarray): laboratory energy of particles in GeV
    Returns:
      (numpy.ndarray): mean (natural) logarithmic mass
    """
    sum_weight = 0.0

    for cid in self.nucleus_ids:
        if cid == 14:
            continue  # p has lnA = 0
        sum_weight += np.log(self.Z_A(cid)[1]) * self.nucleus_flux(cid, E)

    return sum_weight / self.total_flux(E)

Z_A(corsika_id)

Returns charge :math:Z and mass number :math:A corresponding to corsika_id

Parameters:

Name Type Description Default
corsika_id int

corsika id of nucleus/mass group

 required

Returns: (int,int): (Z, A) tuple

Source code in crflux/models.py 
274
275
276
277
278
279
280
281
282
283
284
285
286
287
def Z_A(self, corsika_id):
    """Returns charge :math:`Z` and mass number :math:`A` corresponding
    to ``corsika_id``

    Args:
      corsika_id (int): corsika id of nucleus/mass group
    Returns:
      (int,int): (Z, A) tuple
    """
    Z, A = 1, 1
    if corsika_id > 14:
        Z = corsika_id % 100
        A = (corsika_id - Z) // 100
    return Z, A

H3a_polygonato

Bases: HillasGaisser2012

Modified version of Gaisser, T.K., Astroparticle Physics 35, 801 (2012).

Model is based on Hillas ideas and eye-ball fits by T.K. Gaisser. H3a is a 3-'peters cycle' 5 mass group model with mixed composition above the ankle. H4a has protons only in the 3. component.

This version has a five component poly-gonato at lower energies and resembles H3/4a at energies above the knee.

Parameters:

Name Type Description Default
model str

can be either H3a or H4a.

 'H3a'
Source code in crflux/models.py 
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
class H3a_polygonato(HillasGaisser2012):
    """Modified version of Gaisser, T.K., Astroparticle
    Physics 35, 801 (2012).

    Model is based on Hillas ideas and eye-ball fits by T.K. Gaisser.
    H3a is a 3-'peters cycle' 5 mass group model with mixed
    composition above the ankle. H4a has protons only in the
    3. component.

    This version has a five component poly-gonato at lower energies
    and resembles H3/4a at energies above the knee.

    Args:
      model (str): can be either H3a or H4a.
    """

    def __init__(self, model="H3a", **kwargs):
        super().__init__(model, **kwargs)
        self.rid_cutoff[1] = 4.49e6
        self.rid_cutoff[2] = 30e6
        self.rid_cutoff[3] = 2e9
        self.params[14][1] = (11800, 1.71, 1)  # H
        self.params[402][1] = (4750, 1.64, 2)  # He
        self.params[1206][1] = (3860, 1.67, 6)  # CNO
        self.params[2814][1] = (3120, 1.70, 14)  # MgAlSi
        self.params[5426][1] = (1080, 1.55, 26)  # Fe

        self.params[14][2] = (11.8, 1.4, 1)  # H
        self.params[402][2] = (11.8, 1.4, 2)  # He
        self.params[1206][2] = (7.88, 1.4, 6)  # CNO
        self.params[2814][2] = (7.88, 1.4, 14)  # MgAlSi
        self.params[5426][2] = (7.88, 1.4, 26)  # Fe

nucleus_flux(corsika_id, E)

Returns the flux of nuclei corresponding to the corsika_id at energy E, with optional geomagnetic cutoff.

Parameters:

Name Type Description Default
corsika_id int

see :mod:crflux for description.

 required
E float or ndarray

laboratory energy of nucleus in GeV

 required

Returns: (numpy.ndarray): flux of single nucleus type :math:\Phi_{nucleus} in :math:(\text{m}^2 \text{s sr GeV})^{-1}

Source code in crflux/models.py 
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
def nucleus_flux(self, corsika_id, E):
    """Returns the flux of nuclei corresponding to
    the ``corsika_id`` at energy ``E``, with optional geomagnetic cutoff.

    Args:
      corsika_id (int): see :mod:`crflux` for description.
      E (float or numpy.ndarray): laboratory energy of nucleus in GeV
    Returns:
      (numpy.ndarray): flux of single nucleus type :math:`\\Phi_{nucleus}`
      in :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
    """
    E = np.atleast_1d(np.asarray(E, dtype=float))
    flux = self._nucleus_flux(corsika_id, E)
    if self.geomagnetic_cutoff is not None:
        Z, A = self.Z_A(corsika_id)
        E_cut = Z * self.geomagnetic_cutoff
        flux = np.where(E >= E_cut, flux, 0.0)
    return flux

total_flux(E)

Returns total flux of nuclei, the "all-particle-flux".

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of particles in GeV

 required

Returns: (numpy.ndarray): particle flux in :math:\Phi_{particles} in :math:(\text{m}^2 \text{s sr GeV})^{-1}

Source code in crflux/models.py 
107
108
109
110
111
112
113
114
115
116
def total_flux(self, E):
    """Returns total flux of nuclei, the "all-particle-flux".

    Args:
      E (float or numpy.ndarray): laboratory energy of particles in GeV
    Returns:
      (numpy.ndarray): particle flux in :math:`\\Phi_{particles}` in
      :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
    """
    return sum(self.nucleus_flux(corsika_id, E) for corsika_id in self.nucleus_ids)

tot_nucleon_flux(E)

Returns total flux of nucleons, the "all-nucleon-flux".

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required

Returns: (numpy.ndarray): nucleon flux :math:\Phi_{nucleons} in :math:(\text{m}^2 \text{s sr GeV})^{-1}

Source code in crflux/models.py 
118
119
120
121
122
123
124
125
126
127
128
129
130
131
def tot_nucleon_flux(self, E):
    """Returns total flux of nucleons, the "all-nucleon-flux".

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
    Returns:
      (numpy.ndarray): nucleon flux :math:`\\Phi_{nucleons}` in
      :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
    """
    return sum(
        self.Z_A(corsika_id)[1] ** 2.0
        * self.nucleus_flux(corsika_id, E * self.Z_A(corsika_id)[1])
        for corsika_id in self.nucleus_ids
    )

nucleon_gamma(E, rel_delta=0.01)

Returns differential spectral index of all-nucleon-flux obtained from a numerical derivative.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required
rel_delta float

range of derivative relative to log10(E)

 0.01

Returns: (numpy.ndarray): spectral index :math:\gamma of nucleons

Source code in crflux/models.py 
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
def nucleon_gamma(self, E, rel_delta=0.01):
    """Returns differential spectral index of all-nucleon-flux
    obtained from a numerical derivative.

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
      rel_delta (float): range of derivative relative to log10(E)
    Returns:
      (numpy.ndarray): spectral index :math:`\\gamma` of nucleons
    """
    delta = rel_delta * E
    fl = self.tot_nucleon_flux
    return np.log10(fl(E + delta) / fl(E - delta)) / np.log10(
        (E + delta) / (E - delta)
    )

nucleus_gamma(E, corsika_id, rel_delta=0.01)

Returns differential spectral index of nuclei obtained from a numerical derivative.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nuclei in GeV

 required
corsika_id int

corsika id of nucleus/mass group

 required
rel_delta float

range of derivative relative to log10(E)

 0.01

Returns: (numpy.ndarray): spectral index :math:\gamma of nuclei

Source code in crflux/models.py 
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
def nucleus_gamma(self, E, corsika_id, rel_delta=0.01):
    """Returns differential spectral index of nuclei
    obtained from a numerical derivative.

    Args:
      E (float or numpy.ndarray): laboratory energy of nuclei in GeV
      corsika_id (int): corsika id of nucleus/mass group
      rel_delta (float): range of derivative relative to log10(E)
    Returns:
      (numpy.ndarray): spectral index :math:`\\gamma` of nuclei
    """
    delta = rel_delta * E
    return np.log10(
        self.nucleus_flux(corsika_id, E + delta)
        / self.nucleus_flux(corsika_id, E - delta)
    ) / np.log10((E + delta) / (E - delta))

delta_0(E)

Returns proton excess.

The proton excess is defined as :math:\delta_0 = \frac{\Phi_p - \Phi_n}{\Phi_p + \Phi_n}.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required

Returns: (numpy.ndarray): proton excess :math:\delta_0

Source code in crflux/models.py 
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
def delta_0(self, E):
    """Returns proton excess.

    The proton excess is defined as
    :math:`\\delta_0 = \\frac{\\Phi_p - \\Phi_n}{\\Phi_p + \\Phi_n}`.

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
    Returns:
      (numpy.ndarray): proton excess :math:`\\delta_0`
    """
    p_0 = 0.0
    n_0 = 0.0

    p_0 += self.nucleus_flux(14, E)

    p_0 += 2.0**2 * self.nucleus_flux(402, E * 4.0)
    n_0 += 2.0**2 * self.nucleus_flux(402, E * 4.0)

    p_0 += 6.0**2 * self.nucleus_flux(1206, E * 12.0)
    n_0 += 6.0**2 * self.nucleus_flux(1206, E * 12.0)

    p_0 += 14.0**2 * self.nucleus_flux(2814, E * 28.0)
    n_0 += 14.0**2 * self.nucleus_flux(2814, E * 28.0)

    a_fe = self.Z_A(5426)[1]
    p_0 += 26.0**2 * self.nucleus_flux(5426, E * a_fe)
    n_0 += 26.0**2 * self.nucleus_flux(5426, E * a_fe)

    return (p_0 - n_0) / (p_0 + n_0)

p_and_n_flux(E)

Returns tuple with proton fraction, proton flux and neutron flux.

The proton fraction is defined as :math:\frac{\Phi_p}{\Phi_p + \Phi_n}.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required

Returns: (float,float,float): proton fraction, proton flux, neutron flux

Source code in crflux/models.py 
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
def p_and_n_flux(self, E):
    """Returns tuple with proton fraction, proton flux and neutron flux.

    The proton fraction is defined as :math:`\\frac{\\Phi_p}{\\Phi_p + \\Phi_n}`.

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
    Returns:
      (float,float,float): proton fraction, proton flux, neutron flux
    """
    za = self.Z_A

    p_flux = sum(
        za(corsika_id)[0]
        * za(corsika_id)[1]
        * self.nucleus_flux(corsika_id, E * za(corsika_id)[1])
        for corsika_id in self.nucleus_ids
    )

    n_flux = sum(
        (za(corsika_id)[1] - za(corsika_id)[0])
        * za(corsika_id)[1]
        * self.nucleus_flux(corsika_id, E * za(corsika_id)[1])
        for corsika_id in self.nucleus_ids
    )

    total = p_flux + n_flux
    p_frac = np.where(total > 0, p_flux / total, 0.0)
    return p_frac, p_flux, n_flux

lnA(E)

Returns mean logarithmic mass .

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of particles in GeV

 required

Returns: (numpy.ndarray): mean (natural) logarithmic mass

Source code in crflux/models.py 
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
def lnA(self, E):
    """Returns mean logarithmic mass <ln A>.

    Args:
      E (float or numpy.ndarray): laboratory energy of particles in GeV
    Returns:
      (numpy.ndarray): mean (natural) logarithmic mass
    """
    sum_weight = 0.0

    for cid in self.nucleus_ids:
        if cid == 14:
            continue  # p has lnA = 0
        sum_weight += np.log(self.Z_A(cid)[1]) * self.nucleus_flux(cid, E)

    return sum_weight / self.total_flux(E)

Z_A(corsika_id)

Returns charge :math:Z and mass number :math:A corresponding to corsika_id

Parameters:

Name Type Description Default
corsika_id int

corsika id of nucleus/mass group

 required

Returns: (int,int): (Z, A) tuple

Source code in crflux/models.py 
274
275
276
277
278
279
280
281
282
283
284
285
286
287
def Z_A(self, corsika_id):
    """Returns charge :math:`Z` and mass number :math:`A` corresponding
    to ``corsika_id``

    Args:
      corsika_id (int): corsika id of nucleus/mass group
    Returns:
      (int,int): (Z, A) tuple
    """
    Z, A = 1, 1
    if corsika_id > 14:
        Z = corsika_id % 100
        A = (corsika_id - Z) // 100
    return Z, A

GaisserStanevTilav

Bases: PrimaryFlux

T. K. Gaisser, T. Stanev, and S. Tilav, arXiv:1303.3565, (2013).

Parameters:

Name Type Description Default
model str

3-gen or 4-gen

 '3-gen'
include_heavy_in_total bool

if True, the total flux will include the heavy component groups (Te, Hg) as well. Default is False for legacy compatibility.

 False

Raises:

Type Description
Exception

if model not properly specified.
Source code in crflux/models.py 
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
class GaisserStanevTilav(PrimaryFlux):
    """T. K. Gaisser, T. Stanev, and S. Tilav, arXiv:1303.3565, (2013).

    Args:
      model (str): 3-gen or 4-gen
      include_heavy_in_total (bool, optional): if True, the total flux will
        include the heavy component groups (Te, Hg) as well. Default is
        False for legacy compatibility.

    Raises:
      Exception: if ``model`` not properly specified.
    """

    def __init__(self, model="3-gen", include_heavy_in_total=False, **kwargs):
        super().__init__(**kwargs)

        self.name = "GST (" + model + ")"
        self.sname = "GST" + model[0]
        self.model = model
        self.params = {}
        self.rid_cutoff = {}

        self.rid_cutoff[1] = 120e3
        self.rid_cutoff[2] = 4e6

        mass_comp = [14, 402, 1206, 1608, 5426, 12852, 20180]
        for mcomp in mass_comp:
            self.params[mcomp] = {}

        self.params[14][1] = (7000, 1.66, 1)  # H
        self.params[402][1] = (3200, 1.58, 2)  # He
        self.params[1206][1] = (100, 1.4, 6)  # C
        self.params[1608][1] = (130, 1.4, 8)  # O
        self.params[5426][1] = (60, 1.3, 26)  # Fe
        self.params[12852][1] = (0, 1.0, 52)  # Te
        self.params[20180][1] = (0, 1.0, 80)  # Hg

        self.params[14][2] = (150, 1.4, 1)  # H
        self.params[402][2] = (65, 1.3, 2)  # He
        self.params[1206][2] = (6, 1.3, 6)  # C
        self.params[1608][2] = (7, 1.3, 8)  # O

        if self.model == "3-gen":
            self.params[5426][2] = (2.3, 1.2, 26)  # Fe
            self.params[12852][2] = (0.1, 1.2, 52)  # Te
            self.params[20180][2] = (0.4, 1.2, 80)  # Hg
            self.rid_cutoff[3] = 1.3e9

            self.params[14][3] = (14, 1.4, 1)  # H
            self.params[402][3] = (0, 1.4, 2)  # He
            self.params[1206][3] = (0, 1.4, 6)  # CNO
            self.params[1608][3] = (0, 1.3, 8)  # O
            self.params[5426][3] = (0.025, 1.2, 26)  # Fe
            self.params[12852][3] = (0, 1.0, 52)  # Te
            self.params[20180][3] = (0, 1.0, 80)  # Hg

        elif self.model == "4-gen":
            self.params[5426][2] = (2.1, 1.2, 26)  # Fe
            self.params[12852][2] = (0.1, 1.2, 52)  # Te
            self.params[20180][2] = (0.53, 1.2, 80)  # Hg

            self.rid_cutoff[3] = 1.5e9
            self.params[14][3] = (12.0, 1.4, 1)  # H
            self.params[402][3] = (0, 1.4, 2)  # He
            self.params[1206][3] = (0, 1.4, 6)  # CNO
            self.params[1608][3] = (0, 1.3, 8)  # O
            self.params[5426][3] = (0.011, 1.2, 26)  # Fe
            self.params[12852][3] = (0, 1.0, 52)  # Te
            self.params[20180][3] = (0, 1.0, 80)  # Hg

            self.rid_cutoff[4] = 40e9
            self.params[14][4] = (1.2, 1.4, 1)  # H
            self.params[402][4] = (0, 0, 2)  # He
            self.params[1206][4] = (0, 0, 6)  # CNO
            self.params[1608][4] = (0, 0, 8)  # O
            self.params[5426][4] = (0, 0, 26)  # Fe
            self.params[12852][4] = (0, 1.0, 52)  # Te
            self.params[20180][4] = (0, 1.0, 80)  # Hg
        else:
            raise Exception("GaisserStanevTilav(): Unknown model version.")
        if include_heavy_in_total:
            self.nucleus_ids = list(self.params.keys())
        else:
            self.nucleus_ids = [
                k for k in self.params.keys() if k not in [12852, 20180]
            ]

    def _nucleus_flux(self, corsika_id, E):
        corsika_id = self._find_nearby_id(corsika_id)

        flux = 0.0
        ngen = 0

        if self.model == "3-gen":
            ngen = 4
        elif self.model == "4-gen":
            ngen = 5
        else:
            raise Exception("GaisserStanevTilav(): Unknown model type.")

        for i in range(1, ngen):
            p = self.params[corsika_id][i]
            flux += p[0] * E ** (-p[1] - 1.0) * np.exp(-E / p[2] / self.rid_cutoff[i])
        return flux

nucleus_flux(corsika_id, E)

Returns the flux of nuclei corresponding to the corsika_id at energy E, with optional geomagnetic cutoff.

Parameters:

Name Type Description Default
corsika_id int

see :mod:crflux for description.

 required
E float or ndarray

laboratory energy of nucleus in GeV

 required

Returns: (numpy.ndarray): flux of single nucleus type :math:\Phi_{nucleus} in :math:(\text{m}^2 \text{s sr GeV})^{-1}

Source code in crflux/models.py 
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
def nucleus_flux(self, corsika_id, E):
    """Returns the flux of nuclei corresponding to
    the ``corsika_id`` at energy ``E``, with optional geomagnetic cutoff.

    Args:
      corsika_id (int): see :mod:`crflux` for description.
      E (float or numpy.ndarray): laboratory energy of nucleus in GeV
    Returns:
      (numpy.ndarray): flux of single nucleus type :math:`\\Phi_{nucleus}`
      in :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
    """
    E = np.atleast_1d(np.asarray(E, dtype=float))
    flux = self._nucleus_flux(corsika_id, E)
    if self.geomagnetic_cutoff is not None:
        Z, A = self.Z_A(corsika_id)
        E_cut = Z * self.geomagnetic_cutoff
        flux = np.where(E >= E_cut, flux, 0.0)
    return flux

total_flux(E)

Returns total flux of nuclei, the "all-particle-flux".

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of particles in GeV

 required

Returns: (numpy.ndarray): particle flux in :math:\Phi_{particles} in :math:(\text{m}^2 \text{s sr GeV})^{-1}

Source code in crflux/models.py 
107
108
109
110
111
112
113
114
115
116
def total_flux(self, E):
    """Returns total flux of nuclei, the "all-particle-flux".

    Args:
      E (float or numpy.ndarray): laboratory energy of particles in GeV
    Returns:
      (numpy.ndarray): particle flux in :math:`\\Phi_{particles}` in
      :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
    """
    return sum(self.nucleus_flux(corsika_id, E) for corsika_id in self.nucleus_ids)

tot_nucleon_flux(E)

Returns total flux of nucleons, the "all-nucleon-flux".

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required

Returns: (numpy.ndarray): nucleon flux :math:\Phi_{nucleons} in :math:(\text{m}^2 \text{s sr GeV})^{-1}

Source code in crflux/models.py 
118
119
120
121
122
123
124
125
126
127
128
129
130
131
def tot_nucleon_flux(self, E):
    """Returns total flux of nucleons, the "all-nucleon-flux".

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
    Returns:
      (numpy.ndarray): nucleon flux :math:`\\Phi_{nucleons}` in
      :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
    """
    return sum(
        self.Z_A(corsika_id)[1] ** 2.0
        * self.nucleus_flux(corsika_id, E * self.Z_A(corsika_id)[1])
        for corsika_id in self.nucleus_ids
    )

nucleon_gamma(E, rel_delta=0.01)

Returns differential spectral index of all-nucleon-flux obtained from a numerical derivative.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required
rel_delta float

range of derivative relative to log10(E)

 0.01

Returns: (numpy.ndarray): spectral index :math:\gamma of nucleons

Source code in crflux/models.py 
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
def nucleon_gamma(self, E, rel_delta=0.01):
    """Returns differential spectral index of all-nucleon-flux
    obtained from a numerical derivative.

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
      rel_delta (float): range of derivative relative to log10(E)
    Returns:
      (numpy.ndarray): spectral index :math:`\\gamma` of nucleons
    """
    delta = rel_delta * E
    fl = self.tot_nucleon_flux
    return np.log10(fl(E + delta) / fl(E - delta)) / np.log10(
        (E + delta) / (E - delta)
    )

nucleus_gamma(E, corsika_id, rel_delta=0.01)

Returns differential spectral index of nuclei obtained from a numerical derivative.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nuclei in GeV

 required
corsika_id int

corsika id of nucleus/mass group

 required
rel_delta float

range of derivative relative to log10(E)

 0.01

Returns: (numpy.ndarray): spectral index :math:\gamma of nuclei

Source code in crflux/models.py 
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
def nucleus_gamma(self, E, corsika_id, rel_delta=0.01):
    """Returns differential spectral index of nuclei
    obtained from a numerical derivative.

    Args:
      E (float or numpy.ndarray): laboratory energy of nuclei in GeV
      corsika_id (int): corsika id of nucleus/mass group
      rel_delta (float): range of derivative relative to log10(E)
    Returns:
      (numpy.ndarray): spectral index :math:`\\gamma` of nuclei
    """
    delta = rel_delta * E
    return np.log10(
        self.nucleus_flux(corsika_id, E + delta)
        / self.nucleus_flux(corsika_id, E - delta)
    ) / np.log10((E + delta) / (E - delta))

delta_0(E)

Returns proton excess.

The proton excess is defined as :math:\delta_0 = \frac{\Phi_p - \Phi_n}{\Phi_p + \Phi_n}.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required

Returns: (numpy.ndarray): proton excess :math:\delta_0

Source code in crflux/models.py 
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
def delta_0(self, E):
    """Returns proton excess.

    The proton excess is defined as
    :math:`\\delta_0 = \\frac{\\Phi_p - \\Phi_n}{\\Phi_p + \\Phi_n}`.

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
    Returns:
      (numpy.ndarray): proton excess :math:`\\delta_0`
    """
    p_0 = 0.0
    n_0 = 0.0

    p_0 += self.nucleus_flux(14, E)

    p_0 += 2.0**2 * self.nucleus_flux(402, E * 4.0)
    n_0 += 2.0**2 * self.nucleus_flux(402, E * 4.0)

    p_0 += 6.0**2 * self.nucleus_flux(1206, E * 12.0)
    n_0 += 6.0**2 * self.nucleus_flux(1206, E * 12.0)

    p_0 += 14.0**2 * self.nucleus_flux(2814, E * 28.0)
    n_0 += 14.0**2 * self.nucleus_flux(2814, E * 28.0)

    a_fe = self.Z_A(5426)[1]
    p_0 += 26.0**2 * self.nucleus_flux(5426, E * a_fe)
    n_0 += 26.0**2 * self.nucleus_flux(5426, E * a_fe)

    return (p_0 - n_0) / (p_0 + n_0)

p_and_n_flux(E)

Returns tuple with proton fraction, proton flux and neutron flux.

The proton fraction is defined as :math:\frac{\Phi_p}{\Phi_p + \Phi_n}.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required

Returns: (float,float,float): proton fraction, proton flux, neutron flux

Source code in crflux/models.py 
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
def p_and_n_flux(self, E):
    """Returns tuple with proton fraction, proton flux and neutron flux.

    The proton fraction is defined as :math:`\\frac{\\Phi_p}{\\Phi_p + \\Phi_n}`.

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
    Returns:
      (float,float,float): proton fraction, proton flux, neutron flux
    """
    za = self.Z_A

    p_flux = sum(
        za(corsika_id)[0]
        * za(corsika_id)[1]
        * self.nucleus_flux(corsika_id, E * za(corsika_id)[1])
        for corsika_id in self.nucleus_ids
    )

    n_flux = sum(
        (za(corsika_id)[1] - za(corsika_id)[0])
        * za(corsika_id)[1]
        * self.nucleus_flux(corsika_id, E * za(corsika_id)[1])
        for corsika_id in self.nucleus_ids
    )

    total = p_flux + n_flux
    p_frac = np.where(total > 0, p_flux / total, 0.0)
    return p_frac, p_flux, n_flux

lnA(E)

Returns mean logarithmic mass .

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of particles in GeV

 required

Returns: (numpy.ndarray): mean (natural) logarithmic mass

Source code in crflux/models.py 
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
def lnA(self, E):
    """Returns mean logarithmic mass <ln A>.

    Args:
      E (float or numpy.ndarray): laboratory energy of particles in GeV
    Returns:
      (numpy.ndarray): mean (natural) logarithmic mass
    """
    sum_weight = 0.0

    for cid in self.nucleus_ids:
        if cid == 14:
            continue  # p has lnA = 0
        sum_weight += np.log(self.Z_A(cid)[1]) * self.nucleus_flux(cid, E)

    return sum_weight / self.total_flux(E)

Z_A(corsika_id)

Returns charge :math:Z and mass number :math:A corresponding to corsika_id

Parameters:

Name Type Description Default
corsika_id int

corsika id of nucleus/mass group

 required

Returns: (int,int): (Z, A) tuple

Source code in crflux/models.py 
274
275
276
277
278
279
280
281
282
283
284
285
286
287
def Z_A(self, corsika_id):
    """Returns charge :math:`Z` and mass number :math:`A` corresponding
    to ``corsika_id``

    Args:
      corsika_id (int): corsika id of nucleus/mass group
    Returns:
      (int,int): (Z, A) tuple
    """
    Z, A = 1, 1
    if corsika_id > 14:
        Z = corsika_id % 100
        A = (corsika_id - Z) // 100
    return Z, A

CombinedGHandHG

Bases: PrimaryFlux

A. Fedynitch, J. Becker Tjus, and P. Desiati, Phys. Rev. D 86, 114024 (2012).

In the paper the high energy models were called cHGm for GH+H3a and cHGp for GH+H4a. The names have been changed to use the quite unintuitive names H3a and H4a in ongoing literature.

Source code in crflux/models.py 
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
class CombinedGHandHG(PrimaryFlux):
    """A. Fedynitch, J. Becker Tjus, and P. Desiati,
    Phys. Rev. D 86, 114024 (2012).

    In the paper the high energy models were called cHGm for GH+H3a
    and cHGp for GH+H4a. The names have been changed to use the quite
    unintuitive names H3a and H4a in ongoing literature.
    """

    def __init__(self, model="H3a", **kwargs):
        super().__init__(**kwargs)
        self.name = "comb. GH and " + model
        self.sname = "c" + model
        self.params = {}
        self.leModel = GaisserHonda(**kwargs)
        self.heModel = HillasGaisser2012(model, **kwargs)
        self.heCutOff = 1e5
        cid_list = [14, 402, 1206, 2814, 5426]

        # Store low- to high-energy model transitions in params
        for cid in cid_list:
            self.params[cid] = self.find_transition(cid)
        self.nucleus_ids = list(self.params.keys())

    def find_transition(self, corsika_id):
        from scipy.optimize import fsolve

        def func(logE):
            return self.leModel._nucleus_flux(
                corsika_id, np.atleast_1d(10**logE)
            ) - self.heModel._nucleus_flux(corsika_id, np.atleast_1d(10**logE))

        result = fsolve(func, 3.1)
        return 10 ** result[0]

    def _nucleus_flux(self, corsika_id, E):
        corsika_id = self._find_nearby_id(corsika_id)
        E = np.atleast_1d(np.asarray(E, dtype=float))
        result = np.empty_like(E)
        lo = E < self.params[corsika_id]
        hi = ~lo
        if np.any(lo):
            result[lo] = self.leModel._nucleus_flux(corsika_id, E[lo])
        if np.any(hi):
            result[hi] = self.heModel._nucleus_flux(corsika_id, E[hi])
        return result

nucleus_flux(corsika_id, E)

Returns the flux of nuclei corresponding to the corsika_id at energy E, with optional geomagnetic cutoff.

Parameters:

Name Type Description Default
corsika_id int

see :mod:crflux for description.

 required
E float or ndarray

laboratory energy of nucleus in GeV

 required

Returns: (numpy.ndarray): flux of single nucleus type :math:\Phi_{nucleus} in :math:(\text{m}^2 \text{s sr GeV})^{-1}

Source code in crflux/models.py 
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
def nucleus_flux(self, corsika_id, E):
    """Returns the flux of nuclei corresponding to
    the ``corsika_id`` at energy ``E``, with optional geomagnetic cutoff.

    Args:
      corsika_id (int): see :mod:`crflux` for description.
      E (float or numpy.ndarray): laboratory energy of nucleus in GeV
    Returns:
      (numpy.ndarray): flux of single nucleus type :math:`\\Phi_{nucleus}`
      in :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
    """
    E = np.atleast_1d(np.asarray(E, dtype=float))
    flux = self._nucleus_flux(corsika_id, E)
    if self.geomagnetic_cutoff is not None:
        Z, A = self.Z_A(corsika_id)
        E_cut = Z * self.geomagnetic_cutoff
        flux = np.where(E >= E_cut, flux, 0.0)
    return flux

total_flux(E)

Returns total flux of nuclei, the "all-particle-flux".

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of particles in GeV

 required

Returns: (numpy.ndarray): particle flux in :math:\Phi_{particles} in :math:(\text{m}^2 \text{s sr GeV})^{-1}

Source code in crflux/models.py 
107
108
109
110
111
112
113
114
115
116
def total_flux(self, E):
    """Returns total flux of nuclei, the "all-particle-flux".

    Args:
      E (float or numpy.ndarray): laboratory energy of particles in GeV
    Returns:
      (numpy.ndarray): particle flux in :math:`\\Phi_{particles}` in
      :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
    """
    return sum(self.nucleus_flux(corsika_id, E) for corsika_id in self.nucleus_ids)

tot_nucleon_flux(E)

Returns total flux of nucleons, the "all-nucleon-flux".

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required

Returns: (numpy.ndarray): nucleon flux :math:\Phi_{nucleons} in :math:(\text{m}^2 \text{s sr GeV})^{-1}

Source code in crflux/models.py 
118
119
120
121
122
123
124
125
126
127
128
129
130
131
def tot_nucleon_flux(self, E):
    """Returns total flux of nucleons, the "all-nucleon-flux".

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
    Returns:
      (numpy.ndarray): nucleon flux :math:`\\Phi_{nucleons}` in
      :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
    """
    return sum(
        self.Z_A(corsika_id)[1] ** 2.0
        * self.nucleus_flux(corsika_id, E * self.Z_A(corsika_id)[1])
        for corsika_id in self.nucleus_ids
    )

nucleon_gamma(E, rel_delta=0.01)

Returns differential spectral index of all-nucleon-flux obtained from a numerical derivative.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required
rel_delta float

range of derivative relative to log10(E)

 0.01

Returns: (numpy.ndarray): spectral index :math:\gamma of nucleons

Source code in crflux/models.py 
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
def nucleon_gamma(self, E, rel_delta=0.01):
    """Returns differential spectral index of all-nucleon-flux
    obtained from a numerical derivative.

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
      rel_delta (float): range of derivative relative to log10(E)
    Returns:
      (numpy.ndarray): spectral index :math:`\\gamma` of nucleons
    """
    delta = rel_delta * E
    fl = self.tot_nucleon_flux
    return np.log10(fl(E + delta) / fl(E - delta)) / np.log10(
        (E + delta) / (E - delta)
    )

nucleus_gamma(E, corsika_id, rel_delta=0.01)

Returns differential spectral index of nuclei obtained from a numerical derivative.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nuclei in GeV

 required
corsika_id int

corsika id of nucleus/mass group

 required
rel_delta float

range of derivative relative to log10(E)

 0.01

Returns: (numpy.ndarray): spectral index :math:\gamma of nuclei

Source code in crflux/models.py 
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
def nucleus_gamma(self, E, corsika_id, rel_delta=0.01):
    """Returns differential spectral index of nuclei
    obtained from a numerical derivative.

    Args:
      E (float or numpy.ndarray): laboratory energy of nuclei in GeV
      corsika_id (int): corsika id of nucleus/mass group
      rel_delta (float): range of derivative relative to log10(E)
    Returns:
      (numpy.ndarray): spectral index :math:`\\gamma` of nuclei
    """
    delta = rel_delta * E
    return np.log10(
        self.nucleus_flux(corsika_id, E + delta)
        / self.nucleus_flux(corsika_id, E - delta)
    ) / np.log10((E + delta) / (E - delta))

delta_0(E)

Returns proton excess.

The proton excess is defined as :math:\delta_0 = \frac{\Phi_p - \Phi_n}{\Phi_p + \Phi_n}.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required

Returns: (numpy.ndarray): proton excess :math:\delta_0

Source code in crflux/models.py 
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
def delta_0(self, E):
    """Returns proton excess.

    The proton excess is defined as
    :math:`\\delta_0 = \\frac{\\Phi_p - \\Phi_n}{\\Phi_p + \\Phi_n}`.

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
    Returns:
      (numpy.ndarray): proton excess :math:`\\delta_0`
    """
    p_0 = 0.0
    n_0 = 0.0

    p_0 += self.nucleus_flux(14, E)

    p_0 += 2.0**2 * self.nucleus_flux(402, E * 4.0)
    n_0 += 2.0**2 * self.nucleus_flux(402, E * 4.0)

    p_0 += 6.0**2 * self.nucleus_flux(1206, E * 12.0)
    n_0 += 6.0**2 * self.nucleus_flux(1206, E * 12.0)

    p_0 += 14.0**2 * self.nucleus_flux(2814, E * 28.0)
    n_0 += 14.0**2 * self.nucleus_flux(2814, E * 28.0)

    a_fe = self.Z_A(5426)[1]
    p_0 += 26.0**2 * self.nucleus_flux(5426, E * a_fe)
    n_0 += 26.0**2 * self.nucleus_flux(5426, E * a_fe)

    return (p_0 - n_0) / (p_0 + n_0)

p_and_n_flux(E)

Returns tuple with proton fraction, proton flux and neutron flux.

The proton fraction is defined as :math:\frac{\Phi_p}{\Phi_p + \Phi_n}.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required

Returns: (float,float,float): proton fraction, proton flux, neutron flux

Source code in crflux/models.py 
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
def p_and_n_flux(self, E):
    """Returns tuple with proton fraction, proton flux and neutron flux.

    The proton fraction is defined as :math:`\\frac{\\Phi_p}{\\Phi_p + \\Phi_n}`.

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
    Returns:
      (float,float,float): proton fraction, proton flux, neutron flux
    """
    za = self.Z_A

    p_flux = sum(
        za(corsika_id)[0]
        * za(corsika_id)[1]
        * self.nucleus_flux(corsika_id, E * za(corsika_id)[1])
        for corsika_id in self.nucleus_ids
    )

    n_flux = sum(
        (za(corsika_id)[1] - za(corsika_id)[0])
        * za(corsika_id)[1]
        * self.nucleus_flux(corsika_id, E * za(corsika_id)[1])
        for corsika_id in self.nucleus_ids
    )

    total = p_flux + n_flux
    p_frac = np.where(total > 0, p_flux / total, 0.0)
    return p_frac, p_flux, n_flux

lnA(E)

Returns mean logarithmic mass .

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of particles in GeV

 required

Returns: (numpy.ndarray): mean (natural) logarithmic mass

Source code in crflux/models.py 
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
def lnA(self, E):
    """Returns mean logarithmic mass <ln A>.

    Args:
      E (float or numpy.ndarray): laboratory energy of particles in GeV
    Returns:
      (numpy.ndarray): mean (natural) logarithmic mass
    """
    sum_weight = 0.0

    for cid in self.nucleus_ids:
        if cid == 14:
            continue  # p has lnA = 0
        sum_weight += np.log(self.Z_A(cid)[1]) * self.nucleus_flux(cid, E)

    return sum_weight / self.total_flux(E)

Z_A(corsika_id)

Returns charge :math:Z and mass number :math:A corresponding to corsika_id

Parameters:

Name Type Description Default
corsika_id int

corsika id of nucleus/mass group

 required

Returns: (int,int): (Z, A) tuple

Source code in crflux/models.py 
274
275
276
277
278
279
280
281
282
283
284
285
286
287
def Z_A(self, corsika_id):
    """Returns charge :math:`Z` and mass number :math:`A` corresponding
    to ``corsika_id``

    Args:
      corsika_id (int): corsika id of nucleus/mass group
    Returns:
      (int,int): (Z, A) tuple
    """
    Z, A = 1, 1
    if corsika_id > 14:
        Z = corsika_id % 100
        A = (corsika_id - Z) // 100
    return Z, A

ZatsepinSokolskaya

Bases: PrimaryFlux

The model was first released in V. I. Zatsepin and N. V. Sokolskaya, Astronomy and Astrophysics 458, 1 (2006).

Later, the PAMELA experiment has fitted the parameters of this model to their data in PAMELA Collaboration, O. Adriani et al., Science 332, 69 (2011). Both versions of parameters can be accessed here.

The model does not describe the flux above the knee. Therefore, the highest energies should not exceed 1-10 PeV.

Parameters:

Name Type Description Default
model str

'default' for original or 'pamela' for PAMELA parameters

 'pamela'
Source code in crflux/models.py 
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
class ZatsepinSokolskaya(PrimaryFlux):
    """The model was first released in V. I. Zatsepin and N. V. Sokolskaya,
    Astronomy and Astrophysics 458, 1 (2006).

    Later, the PAMELA experiment has fitted the parameters of this model
    to their data in PAMELA Collaboration, O. Adriani et al.,
    Science 332, 69 (2011). Both versions of parameters can be accessed here.

    The model does not describe the flux above the knee. Therefore, the
    highest energies should not exceed 1-10 PeV.

    Args:
      model (str): 'default' for original or 'pamela' for PAMELA parameters
    """

    def __init__(self, model="pamela", **kwargs):
        super().__init__(**kwargs)
        if model == "pamela":
            self.name = "Zatsepin-Sokolskaya/Pamela"
            self.sname = "ZSP"
            self.R_0 = 5.5
            self.alpha = (2.3, 2.1, 2.57)
            self.R_max = (8e4, 4e6, 2e2)
            self.gamma = (2.63, 2.43, 2.9)
            self.gamma_k = (8, 4.5, 4.5)
            self.f_norm = {}
            self.f_norm[14] = (7.1e3, 6.25e3, 3.0, 74.0, 1)
            self.f_norm[402] = (9.5e3, 8.5e3, 0.74, 18.0, 2)
            self.f_norm[1206] = (6.75e3, 1.8e3, 30, 5.8, 7)
            self.f_norm[2814] = (5.5e3, 1.5e3, 110, 3.5, 12)
            self.f_norm[5426] = (3.5e3, 1.2e3, 750, 2.4, 20)
            self.m_p = 0.983
        elif model == "default":
            self.name = "Zatsepin-Sokolskaya"
            self.sname = "ZS"
            self.R_0 = 5.5
            self.alpha = (2.3, 2.1, 2.57)
            self.R_max = (8e4, 4e6, 2e2)
            self.gamma = (2.63, 2.43, 2.9)
            self.gamma_k = (8, 4.5, 4.5)
            self.f_norm = {}
            self.f_norm[14] = (1.36e4, 6.25e3, 2.1, 74.0, 1)
            self.f_norm[402] = (8.75e3, 8.5e3, 3.0, 18.0, 2)
            self.f_norm[1206] = (6.75e3, 1.8e3, 30, 5.8, 7)
            self.f_norm[2814] = (5.5e3, 1.5e3, 110, 3.5, 12)
            self.f_norm[5426] = (3.5e3, 1.2e3, 750, 2.4, 20)
            self.m_p = 0.983
        else:
            raise Exception(
                f"{self.__class__.__name__}():: Unknown model selection '{model}'."
            )
        self.nucleus_ids = list(self.f_norm.keys())

    def lamba_esc(self, R):
        return (
            4.2 * (R / self.R_0) ** (-1.0 / 3.0) * (1 + (R / self.R_0) ** (-2.0 / 3.0))
        )

    def Q(self, R, gen):
        return R ** (-self.alpha[gen]) * self.phi(R, gen)

    def phi(self, R, gen):
        return (1 + (R / self.R_max[gen]) ** 2) ** (
            (self.gamma[gen] - self.gamma_k[gen]) / 2.0
        )

    def dR_dE(self, E, corsika_id):
        Z, A = self.Z_A(corsika_id)
        return 1.0 / Z * (E + self.m_p * A) / np.sqrt(E**2 + 2 * self.m_p * A * E)

    def R(self, E, corsika_id):
        Z, A = self.Z_A(corsika_id)
        return 1.0 / Z * np.sqrt(E**2 + 2 * self.m_p * A * E)

    def f_mod(self, E, corsika_id):
        Z = self.Z_A(corsika_id)[0]

        P = (E**2 + 2 * self.m_p * E) / (
            (E + Z * 0.511e-3 * 0.6) ** 2 + 2 * self.m_p * (E + Z * 0.511e-3 * 0.6)
        )
        return self.flux(E + Z * 0.511e-3 * 0.6, corsika_id) * P

    def dN_dR(self, R, corsika_id, gen):
        return (
            self.Q(R, gen)
            * self.lamba_esc(R)
            / (1 + self.lamba_esc(R) / self.f_norm[corsika_id][3])
        )

    def dN_dE(self, E, corsika_id, gen):
        return self.dR_dE(E, corsika_id) * self.dN_dR(
            self.R(E, corsika_id), corsika_id, gen
        )

    def flux(self, E, corsika_id):
        flux = 0.0
        for gen in range(3):
            flux += (
                self.f_norm[corsika_id][gen]
                * 1e4 ** (-2.75)
                / self.dN_dE(1e4, corsika_id, gen)
                * self.dN_dE(E, corsika_id, gen)
            )
        return flux

    def _nucleus_flux(self, corsika_id, E):
        corsika_id = self._find_nearby_id(corsika_id)
        E = np.atleast_1d(np.asarray(E, dtype=float))
        result = np.empty_like(E)
        lo = E < 300
        hi = ~lo
        if np.any(lo):
            result[lo] = self.f_mod(E[lo], corsika_id)
        if np.any(hi):
            result[hi] = self.flux(E[hi], corsika_id)
        return result

nucleus_flux(corsika_id, E)

Returns the flux of nuclei corresponding to the corsika_id at energy E, with optional geomagnetic cutoff.

Parameters:

Name Type Description Default
corsika_id int

see :mod:crflux for description.

 required
E float or ndarray

laboratory energy of nucleus in GeV

 required

Returns: (numpy.ndarray): flux of single nucleus type :math:\Phi_{nucleus} in :math:(\text{m}^2 \text{s sr GeV})^{-1}

Source code in crflux/models.py 
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
def nucleus_flux(self, corsika_id, E):
    """Returns the flux of nuclei corresponding to
    the ``corsika_id`` at energy ``E``, with optional geomagnetic cutoff.

    Args:
      corsika_id (int): see :mod:`crflux` for description.
      E (float or numpy.ndarray): laboratory energy of nucleus in GeV
    Returns:
      (numpy.ndarray): flux of single nucleus type :math:`\\Phi_{nucleus}`
      in :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
    """
    E = np.atleast_1d(np.asarray(E, dtype=float))
    flux = self._nucleus_flux(corsika_id, E)
    if self.geomagnetic_cutoff is not None:
        Z, A = self.Z_A(corsika_id)
        E_cut = Z * self.geomagnetic_cutoff
        flux = np.where(E >= E_cut, flux, 0.0)
    return flux

total_flux(E)

Returns total flux of nuclei, the "all-particle-flux".

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of particles in GeV

 required

Returns: (numpy.ndarray): particle flux in :math:\Phi_{particles} in :math:(\text{m}^2 \text{s sr GeV})^{-1}

Source code in crflux/models.py 
107
108
109
110
111
112
113
114
115
116
def total_flux(self, E):
    """Returns total flux of nuclei, the "all-particle-flux".

    Args:
      E (float or numpy.ndarray): laboratory energy of particles in GeV
    Returns:
      (numpy.ndarray): particle flux in :math:`\\Phi_{particles}` in
      :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
    """
    return sum(self.nucleus_flux(corsika_id, E) for corsika_id in self.nucleus_ids)

tot_nucleon_flux(E)

Returns total flux of nucleons, the "all-nucleon-flux".

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required

Returns: (numpy.ndarray): nucleon flux :math:\Phi_{nucleons} in :math:(\text{m}^2 \text{s sr GeV})^{-1}

Source code in crflux/models.py 
118
119
120
121
122
123
124
125
126
127
128
129
130
131
def tot_nucleon_flux(self, E):
    """Returns total flux of nucleons, the "all-nucleon-flux".

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
    Returns:
      (numpy.ndarray): nucleon flux :math:`\\Phi_{nucleons}` in
      :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
    """
    return sum(
        self.Z_A(corsika_id)[1] ** 2.0
        * self.nucleus_flux(corsika_id, E * self.Z_A(corsika_id)[1])
        for corsika_id in self.nucleus_ids
    )

nucleon_gamma(E, rel_delta=0.01)

Returns differential spectral index of all-nucleon-flux obtained from a numerical derivative.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required
rel_delta float

range of derivative relative to log10(E)

 0.01

Returns: (numpy.ndarray): spectral index :math:\gamma of nucleons

Source code in crflux/models.py 
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
def nucleon_gamma(self, E, rel_delta=0.01):
    """Returns differential spectral index of all-nucleon-flux
    obtained from a numerical derivative.

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
      rel_delta (float): range of derivative relative to log10(E)
    Returns:
      (numpy.ndarray): spectral index :math:`\\gamma` of nucleons
    """
    delta = rel_delta * E
    fl = self.tot_nucleon_flux
    return np.log10(fl(E + delta) / fl(E - delta)) / np.log10(
        (E + delta) / (E - delta)
    )

nucleus_gamma(E, corsika_id, rel_delta=0.01)

Returns differential spectral index of nuclei obtained from a numerical derivative.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nuclei in GeV

 required
corsika_id int

corsika id of nucleus/mass group

 required
rel_delta float

range of derivative relative to log10(E)

 0.01

Returns: (numpy.ndarray): spectral index :math:\gamma of nuclei

Source code in crflux/models.py 
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
def nucleus_gamma(self, E, corsika_id, rel_delta=0.01):
    """Returns differential spectral index of nuclei
    obtained from a numerical derivative.

    Args:
      E (float or numpy.ndarray): laboratory energy of nuclei in GeV
      corsika_id (int): corsika id of nucleus/mass group
      rel_delta (float): range of derivative relative to log10(E)
    Returns:
      (numpy.ndarray): spectral index :math:`\\gamma` of nuclei
    """
    delta = rel_delta * E
    return np.log10(
        self.nucleus_flux(corsika_id, E + delta)
        / self.nucleus_flux(corsika_id, E - delta)
    ) / np.log10((E + delta) / (E - delta))

delta_0(E)

Returns proton excess.

The proton excess is defined as :math:\delta_0 = \frac{\Phi_p - \Phi_n}{\Phi_p + \Phi_n}.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required

Returns: (numpy.ndarray): proton excess :math:\delta_0

Source code in crflux/models.py 
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
def delta_0(self, E):
    """Returns proton excess.

    The proton excess is defined as
    :math:`\\delta_0 = \\frac{\\Phi_p - \\Phi_n}{\\Phi_p + \\Phi_n}`.

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
    Returns:
      (numpy.ndarray): proton excess :math:`\\delta_0`
    """
    p_0 = 0.0
    n_0 = 0.0

    p_0 += self.nucleus_flux(14, E)

    p_0 += 2.0**2 * self.nucleus_flux(402, E * 4.0)
    n_0 += 2.0**2 * self.nucleus_flux(402, E * 4.0)

    p_0 += 6.0**2 * self.nucleus_flux(1206, E * 12.0)
    n_0 += 6.0**2 * self.nucleus_flux(1206, E * 12.0)

    p_0 += 14.0**2 * self.nucleus_flux(2814, E * 28.0)
    n_0 += 14.0**2 * self.nucleus_flux(2814, E * 28.0)

    a_fe = self.Z_A(5426)[1]
    p_0 += 26.0**2 * self.nucleus_flux(5426, E * a_fe)
    n_0 += 26.0**2 * self.nucleus_flux(5426, E * a_fe)

    return (p_0 - n_0) / (p_0 + n_0)

p_and_n_flux(E)

Returns tuple with proton fraction, proton flux and neutron flux.

The proton fraction is defined as :math:\frac{\Phi_p}{\Phi_p + \Phi_n}.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required

Returns: (float,float,float): proton fraction, proton flux, neutron flux

Source code in crflux/models.py 
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
def p_and_n_flux(self, E):
    """Returns tuple with proton fraction, proton flux and neutron flux.

    The proton fraction is defined as :math:`\\frac{\\Phi_p}{\\Phi_p + \\Phi_n}`.

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
    Returns:
      (float,float,float): proton fraction, proton flux, neutron flux
    """
    za = self.Z_A

    p_flux = sum(
        za(corsika_id)[0]
        * za(corsika_id)[1]
        * self.nucleus_flux(corsika_id, E * za(corsika_id)[1])
        for corsika_id in self.nucleus_ids
    )

    n_flux = sum(
        (za(corsika_id)[1] - za(corsika_id)[0])
        * za(corsika_id)[1]
        * self.nucleus_flux(corsika_id, E * za(corsika_id)[1])
        for corsika_id in self.nucleus_ids
    )

    total = p_flux + n_flux
    p_frac = np.where(total > 0, p_flux / total, 0.0)
    return p_frac, p_flux, n_flux

lnA(E)

Returns mean logarithmic mass .

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of particles in GeV

 required

Returns: (numpy.ndarray): mean (natural) logarithmic mass

Source code in crflux/models.py 
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
def lnA(self, E):
    """Returns mean logarithmic mass <ln A>.

    Args:
      E (float or numpy.ndarray): laboratory energy of particles in GeV
    Returns:
      (numpy.ndarray): mean (natural) logarithmic mass
    """
    sum_weight = 0.0

    for cid in self.nucleus_ids:
        if cid == 14:
            continue  # p has lnA = 0
        sum_weight += np.log(self.Z_A(cid)[1]) * self.nucleus_flux(cid, E)

    return sum_weight / self.total_flux(E)

Z_A(corsika_id)

Returns charge :math:Z and mass number :math:A corresponding to corsika_id

Parameters:

Name Type Description Default
corsika_id int

corsika id of nucleus/mass group

 required

Returns: (int,int): (Z, A) tuple

Source code in crflux/models.py 
274
275
276
277
278
279
280
281
282
283
284
285
286
287
def Z_A(self, corsika_id):
    """Returns charge :math:`Z` and mass number :math:`A` corresponding
    to ``corsika_id``

    Args:
      corsika_id (int): corsika id of nucleus/mass group
    Returns:
      (int,int): (Z, A) tuple
    """
    Z, A = 1, 1
    if corsika_id > 14:
        Z = corsika_id % 100
        A = (corsika_id - Z) // 100
    return Z, A

GaisserHonda

Bases: PrimaryFlux

5 mass group single power-law model from T.K. Gaisser and M. Honda, Annual Review of Nuclear and Particle Science 52, 153 (2002)

This model was tuned to lower energy baloon data. It fails to describe the flux at and above the knee. A safe range for using it is < 100 TeV/nucleon.

Source code in crflux/models.py 
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
class GaisserHonda(PrimaryFlux):
    """5 mass group single power-law model from T.K. Gaisser and M. Honda,
    Annual Review of Nuclear and Particle Science 52, 153 (2002)

    This model was tuned to lower energy baloon data. It fails to
    describe the flux at and above the knee. A safe range for using
    it is < 100 TeV/nucleon.
    """

    def __init__(self, *args, **kwargs):
        super().__init__(geomagnetic_cutoff=kwargs.pop("geomagnetic_cutoff", None))
        self.name = "Gaisser-Honda"
        self.sname = "GH"
        self.params = {}
        self.params[14] = (2.74, 14900, 2.15, 0.21)
        self.params[402] = (2.64, 600, 1.25, 0.14)
        self.params[1206] = (2.60, 33.2, 0.97, 0.01)
        self.params[2814] = (2.79 + 0.08, 34.2 - 6.0, 2.14, 0.01)
        self.params[5426] = (2.68, 4.45, 3.07, 0.41)
        self.nucleus_ids = list(self.params.keys())

    def _nucleus_flux(self, corsika_id, E):
        corsika_id = self._find_nearby_id(corsika_id)
        A = self.Z_A(corsika_id)[1]

        alpha = self.params[corsika_id][0]
        K = self.params[corsika_id][1]
        b = self.params[corsika_id][2]
        c = self.params[corsika_id][3]

        return K / A * (E / A + b * np.exp(-c * np.sqrt(E / A))) ** -alpha

nucleus_flux(corsika_id, E)

Returns the flux of nuclei corresponding to the corsika_id at energy E, with optional geomagnetic cutoff.

Parameters:

Name Type Description Default
corsika_id int

see :mod:crflux for description.

 required
E float or ndarray

laboratory energy of nucleus in GeV

 required

Returns: (numpy.ndarray): flux of single nucleus type :math:\Phi_{nucleus} in :math:(\text{m}^2 \text{s sr GeV})^{-1}

Source code in crflux/models.py 
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
def nucleus_flux(self, corsika_id, E):
    """Returns the flux of nuclei corresponding to
    the ``corsika_id`` at energy ``E``, with optional geomagnetic cutoff.

    Args:
      corsika_id (int): see :mod:`crflux` for description.
      E (float or numpy.ndarray): laboratory energy of nucleus in GeV
    Returns:
      (numpy.ndarray): flux of single nucleus type :math:`\\Phi_{nucleus}`
      in :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
    """
    E = np.atleast_1d(np.asarray(E, dtype=float))
    flux = self._nucleus_flux(corsika_id, E)
    if self.geomagnetic_cutoff is not None:
        Z, A = self.Z_A(corsika_id)
        E_cut = Z * self.geomagnetic_cutoff
        flux = np.where(E >= E_cut, flux, 0.0)
    return flux

total_flux(E)

Returns total flux of nuclei, the "all-particle-flux".

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of particles in GeV

 required

Returns: (numpy.ndarray): particle flux in :math:\Phi_{particles} in :math:(\text{m}^2 \text{s sr GeV})^{-1}

Source code in crflux/models.py 
107
108
109
110
111
112
113
114
115
116
def total_flux(self, E):
    """Returns total flux of nuclei, the "all-particle-flux".

    Args:
      E (float or numpy.ndarray): laboratory energy of particles in GeV
    Returns:
      (numpy.ndarray): particle flux in :math:`\\Phi_{particles}` in
      :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
    """
    return sum(self.nucleus_flux(corsika_id, E) for corsika_id in self.nucleus_ids)

tot_nucleon_flux(E)

Returns total flux of nucleons, the "all-nucleon-flux".

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required

Returns: (numpy.ndarray): nucleon flux :math:\Phi_{nucleons} in :math:(\text{m}^2 \text{s sr GeV})^{-1}

Source code in crflux/models.py 
118
119
120
121
122
123
124
125
126
127
128
129
130
131
def tot_nucleon_flux(self, E):
    """Returns total flux of nucleons, the "all-nucleon-flux".

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
    Returns:
      (numpy.ndarray): nucleon flux :math:`\\Phi_{nucleons}` in
      :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
    """
    return sum(
        self.Z_A(corsika_id)[1] ** 2.0
        * self.nucleus_flux(corsika_id, E * self.Z_A(corsika_id)[1])
        for corsika_id in self.nucleus_ids
    )

nucleon_gamma(E, rel_delta=0.01)

Returns differential spectral index of all-nucleon-flux obtained from a numerical derivative.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required
rel_delta float

range of derivative relative to log10(E)

 0.01

Returns: (numpy.ndarray): spectral index :math:\gamma of nucleons

Source code in crflux/models.py 
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
def nucleon_gamma(self, E, rel_delta=0.01):
    """Returns differential spectral index of all-nucleon-flux
    obtained from a numerical derivative.

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
      rel_delta (float): range of derivative relative to log10(E)
    Returns:
      (numpy.ndarray): spectral index :math:`\\gamma` of nucleons
    """
    delta = rel_delta * E
    fl = self.tot_nucleon_flux
    return np.log10(fl(E + delta) / fl(E - delta)) / np.log10(
        (E + delta) / (E - delta)
    )

nucleus_gamma(E, corsika_id, rel_delta=0.01)

Returns differential spectral index of nuclei obtained from a numerical derivative.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nuclei in GeV

 required
corsika_id int

corsika id of nucleus/mass group

 required
rel_delta float

range of derivative relative to log10(E)

 0.01

Returns: (numpy.ndarray): spectral index :math:\gamma of nuclei

Source code in crflux/models.py 
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
def nucleus_gamma(self, E, corsika_id, rel_delta=0.01):
    """Returns differential spectral index of nuclei
    obtained from a numerical derivative.

    Args:
      E (float or numpy.ndarray): laboratory energy of nuclei in GeV
      corsika_id (int): corsika id of nucleus/mass group
      rel_delta (float): range of derivative relative to log10(E)
    Returns:
      (numpy.ndarray): spectral index :math:`\\gamma` of nuclei
    """
    delta = rel_delta * E
    return np.log10(
        self.nucleus_flux(corsika_id, E + delta)
        / self.nucleus_flux(corsika_id, E - delta)
    ) / np.log10((E + delta) / (E - delta))

delta_0(E)

Returns proton excess.

The proton excess is defined as :math:\delta_0 = \frac{\Phi_p - \Phi_n}{\Phi_p + \Phi_n}.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required

Returns: (numpy.ndarray): proton excess :math:\delta_0

Source code in crflux/models.py 
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
def delta_0(self, E):
    """Returns proton excess.

    The proton excess is defined as
    :math:`\\delta_0 = \\frac{\\Phi_p - \\Phi_n}{\\Phi_p + \\Phi_n}`.

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
    Returns:
      (numpy.ndarray): proton excess :math:`\\delta_0`
    """
    p_0 = 0.0
    n_0 = 0.0

    p_0 += self.nucleus_flux(14, E)

    p_0 += 2.0**2 * self.nucleus_flux(402, E * 4.0)
    n_0 += 2.0**2 * self.nucleus_flux(402, E * 4.0)

    p_0 += 6.0**2 * self.nucleus_flux(1206, E * 12.0)
    n_0 += 6.0**2 * self.nucleus_flux(1206, E * 12.0)

    p_0 += 14.0**2 * self.nucleus_flux(2814, E * 28.0)
    n_0 += 14.0**2 * self.nucleus_flux(2814, E * 28.0)

    a_fe = self.Z_A(5426)[1]
    p_0 += 26.0**2 * self.nucleus_flux(5426, E * a_fe)
    n_0 += 26.0**2 * self.nucleus_flux(5426, E * a_fe)

    return (p_0 - n_0) / (p_0 + n_0)

p_and_n_flux(E)

Returns tuple with proton fraction, proton flux and neutron flux.

The proton fraction is defined as :math:\frac{\Phi_p}{\Phi_p + \Phi_n}.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required

Returns: (float,float,float): proton fraction, proton flux, neutron flux

Source code in crflux/models.py 
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
def p_and_n_flux(self, E):
    """Returns tuple with proton fraction, proton flux and neutron flux.

    The proton fraction is defined as :math:`\\frac{\\Phi_p}{\\Phi_p + \\Phi_n}`.

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
    Returns:
      (float,float,float): proton fraction, proton flux, neutron flux
    """
    za = self.Z_A

    p_flux = sum(
        za(corsika_id)[0]
        * za(corsika_id)[1]
        * self.nucleus_flux(corsika_id, E * za(corsika_id)[1])
        for corsika_id in self.nucleus_ids
    )

    n_flux = sum(
        (za(corsika_id)[1] - za(corsika_id)[0])
        * za(corsika_id)[1]
        * self.nucleus_flux(corsika_id, E * za(corsika_id)[1])
        for corsika_id in self.nucleus_ids
    )

    total = p_flux + n_flux
    p_frac = np.where(total > 0, p_flux / total, 0.0)
    return p_frac, p_flux, n_flux

lnA(E)

Returns mean logarithmic mass .

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of particles in GeV

 required

Returns: (numpy.ndarray): mean (natural) logarithmic mass

Source code in crflux/models.py 
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
def lnA(self, E):
    """Returns mean logarithmic mass <ln A>.

    Args:
      E (float or numpy.ndarray): laboratory energy of particles in GeV
    Returns:
      (numpy.ndarray): mean (natural) logarithmic mass
    """
    sum_weight = 0.0

    for cid in self.nucleus_ids:
        if cid == 14:
            continue  # p has lnA = 0
        sum_weight += np.log(self.Z_A(cid)[1]) * self.nucleus_flux(cid, E)

    return sum_weight / self.total_flux(E)

Z_A(corsika_id)

Returns charge :math:Z and mass number :math:A corresponding to corsika_id

Parameters:

Name Type Description Default
corsika_id int

corsika id of nucleus/mass group

 required

Returns: (int,int): (Z, A) tuple

Source code in crflux/models.py 
274
275
276
277
278
279
280
281
282
283
284
285
286
287
def Z_A(self, corsika_id):
    """Returns charge :math:`Z` and mass number :math:`A` corresponding
    to ``corsika_id``

    Args:
      corsika_id (int): corsika id of nucleus/mass group
    Returns:
      (int,int): (Z, A) tuple
    """
    Z, A = 1, 1
    if corsika_id > 14:
        Z = corsika_id % 100
        A = (corsika_id - Z) // 100
    return Z, A

Thunman

Bases: PrimaryFlux

Popular broken power-law flux model.

The parameters of this model are taken from the prompt flux calculation paper by M. Thunman, G. Ingelman, and P. Gondolo, Astroparticle Physics 5, 309 (1996). The model contains only protons with a power-law index of -2.7 below the knee, located at 5 PeV, and -3.0 for energies higher than that.

Source code in crflux/models.py 
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
class Thunman(PrimaryFlux):
    """Popular broken power-law flux model.

    The parameters of this model are taken from the prompt flux calculation
    paper by M. Thunman, G. Ingelman, and P. Gondolo,
    Astroparticle Physics 5, 309 (1996).
    The model contains only protons with a power-law index of -2.7 below
    the knee, located at 5 PeV, and -3.0 for energies higher than that.
    """

    name = "Thunman et al. ('96)"
    sname = "TIG"

    def __init__(self, *args, **kwargs):
        super().__init__(geomagnetic_cutoff=kwargs.pop("geomagnetic_cutoff", None))
        self.params = {}
        self.params["low_e"] = (1e4 * 1.7, -2.7)
        self.params["high_e"] = (1e4 * 174, -3.0)
        self.params["trans"] = 5e6

        self.nucleus_ids = [14]

    def _nucleus_flux(self, corsika_id, E):
        """Broken power law spectrum for protons."""
        E = np.atleast_1d(np.asarray(E, dtype=float))
        if corsika_id != 14:
            return np.zeros_like(E)

        le = E < self.params["trans"]
        he = ~le
        result = np.empty_like(E)
        result[le] = self.params["low_e"][0] * E[le] ** self.params["low_e"][1]
        result[he] = self.params["high_e"][0] * E[he] ** self.params["high_e"][1]
        return result

nucleus_flux(corsika_id, E)

Returns the flux of nuclei corresponding to the corsika_id at energy E, with optional geomagnetic cutoff.

Parameters:

Name Type Description Default
corsika_id int

see :mod:crflux for description.

 required
E float or ndarray

laboratory energy of nucleus in GeV

 required

Returns: (numpy.ndarray): flux of single nucleus type :math:\Phi_{nucleus} in :math:(\text{m}^2 \text{s sr GeV})^{-1}

Source code in crflux/models.py 
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
def nucleus_flux(self, corsika_id, E):
    """Returns the flux of nuclei corresponding to
    the ``corsika_id`` at energy ``E``, with optional geomagnetic cutoff.

    Args:
      corsika_id (int): see :mod:`crflux` for description.
      E (float or numpy.ndarray): laboratory energy of nucleus in GeV
    Returns:
      (numpy.ndarray): flux of single nucleus type :math:`\\Phi_{nucleus}`
      in :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
    """
    E = np.atleast_1d(np.asarray(E, dtype=float))
    flux = self._nucleus_flux(corsika_id, E)
    if self.geomagnetic_cutoff is not None:
        Z, A = self.Z_A(corsika_id)
        E_cut = Z * self.geomagnetic_cutoff
        flux = np.where(E >= E_cut, flux, 0.0)
    return flux

total_flux(E)

Returns total flux of nuclei, the "all-particle-flux".

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of particles in GeV

 required

Returns: (numpy.ndarray): particle flux in :math:\Phi_{particles} in :math:(\text{m}^2 \text{s sr GeV})^{-1}

Source code in crflux/models.py 
107
108
109
110
111
112
113
114
115
116
def total_flux(self, E):
    """Returns total flux of nuclei, the "all-particle-flux".

    Args:
      E (float or numpy.ndarray): laboratory energy of particles in GeV
    Returns:
      (numpy.ndarray): particle flux in :math:`\\Phi_{particles}` in
      :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
    """
    return sum(self.nucleus_flux(corsika_id, E) for corsika_id in self.nucleus_ids)

tot_nucleon_flux(E)

Returns total flux of nucleons, the "all-nucleon-flux".

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required

Returns: (numpy.ndarray): nucleon flux :math:\Phi_{nucleons} in :math:(\text{m}^2 \text{s sr GeV})^{-1}

Source code in crflux/models.py 
118
119
120
121
122
123
124
125
126
127
128
129
130
131
def tot_nucleon_flux(self, E):
    """Returns total flux of nucleons, the "all-nucleon-flux".

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
    Returns:
      (numpy.ndarray): nucleon flux :math:`\\Phi_{nucleons}` in
      :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
    """
    return sum(
        self.Z_A(corsika_id)[1] ** 2.0
        * self.nucleus_flux(corsika_id, E * self.Z_A(corsika_id)[1])
        for corsika_id in self.nucleus_ids
    )

nucleon_gamma(E, rel_delta=0.01)

Returns differential spectral index of all-nucleon-flux obtained from a numerical derivative.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required
rel_delta float

range of derivative relative to log10(E)

 0.01

Returns: (numpy.ndarray): spectral index :math:\gamma of nucleons

Source code in crflux/models.py 
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
def nucleon_gamma(self, E, rel_delta=0.01):
    """Returns differential spectral index of all-nucleon-flux
    obtained from a numerical derivative.

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
      rel_delta (float): range of derivative relative to log10(E)
    Returns:
      (numpy.ndarray): spectral index :math:`\\gamma` of nucleons
    """
    delta = rel_delta * E
    fl = self.tot_nucleon_flux
    return np.log10(fl(E + delta) / fl(E - delta)) / np.log10(
        (E + delta) / (E - delta)
    )

nucleus_gamma(E, corsika_id, rel_delta=0.01)

Returns differential spectral index of nuclei obtained from a numerical derivative.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nuclei in GeV

 required
corsika_id int

corsika id of nucleus/mass group

 required
rel_delta float

range of derivative relative to log10(E)

 0.01

Returns: (numpy.ndarray): spectral index :math:\gamma of nuclei

Source code in crflux/models.py 
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
def nucleus_gamma(self, E, corsika_id, rel_delta=0.01):
    """Returns differential spectral index of nuclei
    obtained from a numerical derivative.

    Args:
      E (float or numpy.ndarray): laboratory energy of nuclei in GeV
      corsika_id (int): corsika id of nucleus/mass group
      rel_delta (float): range of derivative relative to log10(E)
    Returns:
      (numpy.ndarray): spectral index :math:`\\gamma` of nuclei
    """
    delta = rel_delta * E
    return np.log10(
        self.nucleus_flux(corsika_id, E + delta)
        / self.nucleus_flux(corsika_id, E - delta)
    ) / np.log10((E + delta) / (E - delta))

delta_0(E)

Returns proton excess.

The proton excess is defined as :math:\delta_0 = \frac{\Phi_p - \Phi_n}{\Phi_p + \Phi_n}.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required

Returns: (numpy.ndarray): proton excess :math:\delta_0

Source code in crflux/models.py 
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
def delta_0(self, E):
    """Returns proton excess.

    The proton excess is defined as
    :math:`\\delta_0 = \\frac{\\Phi_p - \\Phi_n}{\\Phi_p + \\Phi_n}`.

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
    Returns:
      (numpy.ndarray): proton excess :math:`\\delta_0`
    """
    p_0 = 0.0
    n_0 = 0.0

    p_0 += self.nucleus_flux(14, E)

    p_0 += 2.0**2 * self.nucleus_flux(402, E * 4.0)
    n_0 += 2.0**2 * self.nucleus_flux(402, E * 4.0)

    p_0 += 6.0**2 * self.nucleus_flux(1206, E * 12.0)
    n_0 += 6.0**2 * self.nucleus_flux(1206, E * 12.0)

    p_0 += 14.0**2 * self.nucleus_flux(2814, E * 28.0)
    n_0 += 14.0**2 * self.nucleus_flux(2814, E * 28.0)

    a_fe = self.Z_A(5426)[1]
    p_0 += 26.0**2 * self.nucleus_flux(5426, E * a_fe)
    n_0 += 26.0**2 * self.nucleus_flux(5426, E * a_fe)

    return (p_0 - n_0) / (p_0 + n_0)

p_and_n_flux(E)

Returns tuple with proton fraction, proton flux and neutron flux.

The proton fraction is defined as :math:\frac{\Phi_p}{\Phi_p + \Phi_n}.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required

Returns: (float,float,float): proton fraction, proton flux, neutron flux

Source code in crflux/models.py 
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
def p_and_n_flux(self, E):
    """Returns tuple with proton fraction, proton flux and neutron flux.

    The proton fraction is defined as :math:`\\frac{\\Phi_p}{\\Phi_p + \\Phi_n}`.

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
    Returns:
      (float,float,float): proton fraction, proton flux, neutron flux
    """
    za = self.Z_A

    p_flux = sum(
        za(corsika_id)[0]
        * za(corsika_id)[1]
        * self.nucleus_flux(corsika_id, E * za(corsika_id)[1])
        for corsika_id in self.nucleus_ids
    )

    n_flux = sum(
        (za(corsika_id)[1] - za(corsika_id)[0])
        * za(corsika_id)[1]
        * self.nucleus_flux(corsika_id, E * za(corsika_id)[1])
        for corsika_id in self.nucleus_ids
    )

    total = p_flux + n_flux
    p_frac = np.where(total > 0, p_flux / total, 0.0)
    return p_frac, p_flux, n_flux

lnA(E)

Returns mean logarithmic mass .

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of particles in GeV

 required

Returns: (numpy.ndarray): mean (natural) logarithmic mass

Source code in crflux/models.py 
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
def lnA(self, E):
    """Returns mean logarithmic mass <ln A>.

    Args:
      E (float or numpy.ndarray): laboratory energy of particles in GeV
    Returns:
      (numpy.ndarray): mean (natural) logarithmic mass
    """
    sum_weight = 0.0

    for cid in self.nucleus_ids:
        if cid == 14:
            continue  # p has lnA = 0
        sum_weight += np.log(self.Z_A(cid)[1]) * self.nucleus_flux(cid, E)

    return sum_weight / self.total_flux(E)

Z_A(corsika_id)

Returns charge :math:Z and mass number :math:A corresponding to corsika_id

Parameters:

Name Type Description Default
corsika_id int

corsika id of nucleus/mass group

 required

Returns: (int,int): (Z, A) tuple

Source code in crflux/models.py 
274
275
276
277
278
279
280
281
282
283
284
285
286
287
def Z_A(self, corsika_id):
    """Returns charge :math:`Z` and mass number :math:`A` corresponding
    to ``corsika_id``

    Args:
      corsika_id (int): corsika id of nucleus/mass group
    Returns:
      (int,int): (Z, A) tuple
    """
    Z, A = 1, 1
    if corsika_id > 14:
        Z = corsika_id % 100
        A = (corsika_id - Z) // 100
    return Z, A

SimplePowerlaw27

Bases: PrimaryFlux

Simple E**-2.7 parametrization based on values from :class:Thunman below knee.

Source code in crflux/models.py 
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
class SimplePowerlaw27(PrimaryFlux):
    """Simple E**-2.7 parametrization based on values from
    :class:`Thunman` below knee.
    """

    name = r"$E^{-2.7}$"
    sname = r"$E^{-2.7}$"

    def __init__(self, *args, **kwargs):
        super().__init__(geomagnetic_cutoff=kwargs.pop("geomagnetic_cutoff", None))
        self.params = (1e4 * 1.7, -2.7)
        self.nucleus_ids = [14]

    def _nucleus_flux(self, corsika_id, E):
        E = np.atleast_1d(np.asarray(E, dtype=float))
        if corsika_id != 14:
            return np.zeros_like(E)

        return self.params[0] * E ** (self.params[1])

nucleus_flux(corsika_id, E)

Returns the flux of nuclei corresponding to the corsika_id at energy E, with optional geomagnetic cutoff.

Parameters:

Name Type Description Default
corsika_id int

see :mod:crflux for description.

 required
E float or ndarray

laboratory energy of nucleus in GeV

 required

Returns: (numpy.ndarray): flux of single nucleus type :math:\Phi_{nucleus} in :math:(\text{m}^2 \text{s sr GeV})^{-1}

Source code in crflux/models.py 
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
def nucleus_flux(self, corsika_id, E):
    """Returns the flux of nuclei corresponding to
    the ``corsika_id`` at energy ``E``, with optional geomagnetic cutoff.

    Args:
      corsika_id (int): see :mod:`crflux` for description.
      E (float or numpy.ndarray): laboratory energy of nucleus in GeV
    Returns:
      (numpy.ndarray): flux of single nucleus type :math:`\\Phi_{nucleus}`
      in :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
    """
    E = np.atleast_1d(np.asarray(E, dtype=float))
    flux = self._nucleus_flux(corsika_id, E)
    if self.geomagnetic_cutoff is not None:
        Z, A = self.Z_A(corsika_id)
        E_cut = Z * self.geomagnetic_cutoff
        flux = np.where(E >= E_cut, flux, 0.0)
    return flux

total_flux(E)

Returns total flux of nuclei, the "all-particle-flux".

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of particles in GeV

 required

Returns: (numpy.ndarray): particle flux in :math:\Phi_{particles} in :math:(\text{m}^2 \text{s sr GeV})^{-1}

Source code in crflux/models.py 
107
108
109
110
111
112
113
114
115
116
def total_flux(self, E):
    """Returns total flux of nuclei, the "all-particle-flux".

    Args:
      E (float or numpy.ndarray): laboratory energy of particles in GeV
    Returns:
      (numpy.ndarray): particle flux in :math:`\\Phi_{particles}` in
      :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
    """
    return sum(self.nucleus_flux(corsika_id, E) for corsika_id in self.nucleus_ids)

tot_nucleon_flux(E)

Returns total flux of nucleons, the "all-nucleon-flux".

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required

Returns: (numpy.ndarray): nucleon flux :math:\Phi_{nucleons} in :math:(\text{m}^2 \text{s sr GeV})^{-1}

Source code in crflux/models.py 
118
119
120
121
122
123
124
125
126
127
128
129
130
131
def tot_nucleon_flux(self, E):
    """Returns total flux of nucleons, the "all-nucleon-flux".

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
    Returns:
      (numpy.ndarray): nucleon flux :math:`\\Phi_{nucleons}` in
      :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
    """
    return sum(
        self.Z_A(corsika_id)[1] ** 2.0
        * self.nucleus_flux(corsika_id, E * self.Z_A(corsika_id)[1])
        for corsika_id in self.nucleus_ids
    )

nucleon_gamma(E, rel_delta=0.01)

Returns differential spectral index of all-nucleon-flux obtained from a numerical derivative.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required
rel_delta float

range of derivative relative to log10(E)

 0.01

Returns: (numpy.ndarray): spectral index :math:\gamma of nucleons

Source code in crflux/models.py 
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
def nucleon_gamma(self, E, rel_delta=0.01):
    """Returns differential spectral index of all-nucleon-flux
    obtained from a numerical derivative.

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
      rel_delta (float): range of derivative relative to log10(E)
    Returns:
      (numpy.ndarray): spectral index :math:`\\gamma` of nucleons
    """
    delta = rel_delta * E
    fl = self.tot_nucleon_flux
    return np.log10(fl(E + delta) / fl(E - delta)) / np.log10(
        (E + delta) / (E - delta)
    )

nucleus_gamma(E, corsika_id, rel_delta=0.01)

Returns differential spectral index of nuclei obtained from a numerical derivative.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nuclei in GeV

 required
corsika_id int

corsika id of nucleus/mass group

 required
rel_delta float

range of derivative relative to log10(E)

 0.01

Returns: (numpy.ndarray): spectral index :math:\gamma of nuclei

Source code in crflux/models.py 
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
def nucleus_gamma(self, E, corsika_id, rel_delta=0.01):
    """Returns differential spectral index of nuclei
    obtained from a numerical derivative.

    Args:
      E (float or numpy.ndarray): laboratory energy of nuclei in GeV
      corsika_id (int): corsika id of nucleus/mass group
      rel_delta (float): range of derivative relative to log10(E)
    Returns:
      (numpy.ndarray): spectral index :math:`\\gamma` of nuclei
    """
    delta = rel_delta * E
    return np.log10(
        self.nucleus_flux(corsika_id, E + delta)
        / self.nucleus_flux(corsika_id, E - delta)
    ) / np.log10((E + delta) / (E - delta))

delta_0(E)

Returns proton excess.

The proton excess is defined as :math:\delta_0 = \frac{\Phi_p - \Phi_n}{\Phi_p + \Phi_n}.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required

Returns: (numpy.ndarray): proton excess :math:\delta_0

Source code in crflux/models.py 
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
def delta_0(self, E):
    """Returns proton excess.

    The proton excess is defined as
    :math:`\\delta_0 = \\frac{\\Phi_p - \\Phi_n}{\\Phi_p + \\Phi_n}`.

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
    Returns:
      (numpy.ndarray): proton excess :math:`\\delta_0`
    """
    p_0 = 0.0
    n_0 = 0.0

    p_0 += self.nucleus_flux(14, E)

    p_0 += 2.0**2 * self.nucleus_flux(402, E * 4.0)
    n_0 += 2.0**2 * self.nucleus_flux(402, E * 4.0)

    p_0 += 6.0**2 * self.nucleus_flux(1206, E * 12.0)
    n_0 += 6.0**2 * self.nucleus_flux(1206, E * 12.0)

    p_0 += 14.0**2 * self.nucleus_flux(2814, E * 28.0)
    n_0 += 14.0**2 * self.nucleus_flux(2814, E * 28.0)

    a_fe = self.Z_A(5426)[1]
    p_0 += 26.0**2 * self.nucleus_flux(5426, E * a_fe)
    n_0 += 26.0**2 * self.nucleus_flux(5426, E * a_fe)

    return (p_0 - n_0) / (p_0 + n_0)

p_and_n_flux(E)

Returns tuple with proton fraction, proton flux and neutron flux.

The proton fraction is defined as :math:\frac{\Phi_p}{\Phi_p + \Phi_n}.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required

Returns: (float,float,float): proton fraction, proton flux, neutron flux

Source code in crflux/models.py 
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
def p_and_n_flux(self, E):
    """Returns tuple with proton fraction, proton flux and neutron flux.

    The proton fraction is defined as :math:`\\frac{\\Phi_p}{\\Phi_p + \\Phi_n}`.

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
    Returns:
      (float,float,float): proton fraction, proton flux, neutron flux
    """
    za = self.Z_A

    p_flux = sum(
        za(corsika_id)[0]
        * za(corsika_id)[1]
        * self.nucleus_flux(corsika_id, E * za(corsika_id)[1])
        for corsika_id in self.nucleus_ids
    )

    n_flux = sum(
        (za(corsika_id)[1] - za(corsika_id)[0])
        * za(corsika_id)[1]
        * self.nucleus_flux(corsika_id, E * za(corsika_id)[1])
        for corsika_id in self.nucleus_ids
    )

    total = p_flux + n_flux
    p_frac = np.where(total > 0, p_flux / total, 0.0)
    return p_frac, p_flux, n_flux

lnA(E)

Returns mean logarithmic mass .

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of particles in GeV

 required

Returns: (numpy.ndarray): mean (natural) logarithmic mass

Source code in crflux/models.py 
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
def lnA(self, E):
    """Returns mean logarithmic mass <ln A>.

    Args:
      E (float or numpy.ndarray): laboratory energy of particles in GeV
    Returns:
      (numpy.ndarray): mean (natural) logarithmic mass
    """
    sum_weight = 0.0

    for cid in self.nucleus_ids:
        if cid == 14:
            continue  # p has lnA = 0
        sum_weight += np.log(self.Z_A(cid)[1]) * self.nucleus_flux(cid, E)

    return sum_weight / self.total_flux(E)

Z_A(corsika_id)

Returns charge :math:Z and mass number :math:A corresponding to corsika_id

Parameters:

Name Type Description Default
corsika_id int

corsika id of nucleus/mass group

 required

Returns: (int,int): (Z, A) tuple

Source code in crflux/models.py 
274
275
276
277
278
279
280
281
282
283
284
285
286
287
def Z_A(self, corsika_id):
    """Returns charge :math:`Z` and mass number :math:`A` corresponding
    to ``corsika_id``

    Args:
      corsika_id (int): corsika id of nucleus/mass group
    Returns:
      (int,int): (Z, A) tuple
    """
    Z, A = 1, 1
    if corsika_id > 14:
        Z = corsika_id % 100
        A = (corsika_id - Z) // 100
    return Z, A

GlobalSplineFit

Bases: PrimaryFlux

Global Spline Fit (GSF) model through the globalsplinefit package.

Per-nucleus fluxes come from GSFEnergy (total energy per nucleus) and nucleon fluxes from GSFEnergyPerNucleon.p_and_n_flux (total energy per nucleon), both in (m² s sr GeV)^-1. Requires the optional gsf extra.

Use for reference: Fedynitch et al., arXiv:2609.32649 (2026).

Parameters

version : str globalsplinefit model version, e.g. "2026.1", "2025", "2019", "2017". time_interval : tuple or str, optional Solar-modulation interval (YYYYMM, YYYYMM), or "LIS". None uses the Solar Cycle 24 average. geomagnetic_cutoff : float, optional Rigidity cutoff in GV, applied by globalsplinefit per species in rigidity (not the base-class E >= Z R_cut step).

Notes

nucleus_ids lists every species of the version as CORSIKA ids (100 A + Z with A the rounded mass, 14 for protons, 201 for deuterium where the version carries it). gsf is the underlying GSFEnergy model for uncertainties, covariances and composition.

Source code in crflux/models.py 
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
class GlobalSplineFit(PrimaryFlux):
    """Global Spline Fit (GSF) model through the globalsplinefit package.

    Per-nucleus fluxes come from ``GSFEnergy`` (total energy per nucleus) and
    nucleon fluxes from ``GSFEnergyPerNucleon.p_and_n_flux`` (total energy per
    nucleon), both in (m² s sr GeV)^-1. Requires the optional ``gsf`` extra.

    Use for reference: Fedynitch et al., arXiv:2609.32649 (2026).

    Parameters
    ----------
    version : str
        globalsplinefit model version, e.g. "2026.1", "2025", "2019", "2017".
    time_interval : tuple or str, optional
        Solar-modulation interval (YYYYMM, YYYYMM), or "LIS". None uses the
        Solar Cycle 24 average.
    geomagnetic_cutoff : float, optional
        Rigidity cutoff in GV, applied by globalsplinefit per species in
        rigidity (not the base-class E >= Z R_cut step).

    Notes
    -----
    ``nucleus_ids`` lists every species of the version as CORSIKA ids
    (100 A + Z with A the rounded mass, 14 for protons, 201 for deuterium
    where the version carries it). ``gsf`` is the underlying
    ``GSFEnergy`` model for uncertainties, covariances and composition.
    """

    def __init__(self, version="2026.1", time_interval=None, geomagnetic_cutoff=None):
        super().__init__(geomagnetic_cutoff=geomagnetic_cutoff)
        from globalsplinefit import GSFEnergy, GSFEnergyPerNucleon

        opts = {
            "version": version,
            "default_time_interval": time_interval,
            "default_rigidity_cutoff": geomagnetic_cutoff,
        }
        self.gsf = GSFEnergy(**opts)
        self._gsf_nucleon = GSFEnergyPerNucleon(**opts)
        self._species = {self._corsika_id(*sid): sid for sid in self.gsf.species}
        self.nucleus_ids = sorted(self._species)
        self.name = f"Global Spline Fit {self.gsf.version}"
        self.sname = f"GSF{self.gsf.version}"

    @staticmethod
    def _corsika_id(z, a):
        a = int(round(a))
        return 14 if (z, a) == (1, 1) else 100 * a + z

    def nucleus_flux(self, corsika_id, E):
        # the cutoff is already applied by globalsplinefit
        return self._nucleus_flux(corsika_id, np.atleast_1d(np.asarray(E, dtype=float)))

    def _nucleus_flux(self, corsika_id, E):
        if corsika_id not in self._species:
            return np.zeros_like(E)
        return self.gsf.flux(E, self._species[corsika_id])

    def p_and_n_flux(self, E):
        energy = np.atleast_1d(np.asarray(E, dtype=float))
        model = self._gsf_nucleon
        proton, neutron = sum(
            model.p_and_n_flux(energy, g) for g in model.active_groups
        )
        total = proton + neutron
        fraction = np.divide(proton, total, out=np.zeros_like(total), where=total > 0)
        return fraction, proton, neutron

    def tot_nucleon_flux(self, E):
        return np.sum(self.p_and_n_flux(E)[1:], axis=0)

total_flux(E)

Returns total flux of nuclei, the "all-particle-flux".

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of particles in GeV

 required

Returns: (numpy.ndarray): particle flux in :math:\Phi_{particles} in :math:(\text{m}^2 \text{s sr GeV})^{-1}

Source code in crflux/models.py 
107
108
109
110
111
112
113
114
115
116
def total_flux(self, E):
    """Returns total flux of nuclei, the "all-particle-flux".

    Args:
      E (float or numpy.ndarray): laboratory energy of particles in GeV
    Returns:
      (numpy.ndarray): particle flux in :math:`\\Phi_{particles}` in
      :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
    """
    return sum(self.nucleus_flux(corsika_id, E) for corsika_id in self.nucleus_ids)

nucleon_gamma(E, rel_delta=0.01)

Returns differential spectral index of all-nucleon-flux obtained from a numerical derivative.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required
rel_delta float

range of derivative relative to log10(E)

 0.01

Returns: (numpy.ndarray): spectral index :math:\gamma of nucleons

Source code in crflux/models.py 
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
def nucleon_gamma(self, E, rel_delta=0.01):
    """Returns differential spectral index of all-nucleon-flux
    obtained from a numerical derivative.

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
      rel_delta (float): range of derivative relative to log10(E)
    Returns:
      (numpy.ndarray): spectral index :math:`\\gamma` of nucleons
    """
    delta = rel_delta * E
    fl = self.tot_nucleon_flux
    return np.log10(fl(E + delta) / fl(E - delta)) / np.log10(
        (E + delta) / (E - delta)
    )

nucleus_gamma(E, corsika_id, rel_delta=0.01)

Returns differential spectral index of nuclei obtained from a numerical derivative.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nuclei in GeV

 required
corsika_id int

corsika id of nucleus/mass group

 required
rel_delta float

range of derivative relative to log10(E)

 0.01

Returns: (numpy.ndarray): spectral index :math:\gamma of nuclei

Source code in crflux/models.py 
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
def nucleus_gamma(self, E, corsika_id, rel_delta=0.01):
    """Returns differential spectral index of nuclei
    obtained from a numerical derivative.

    Args:
      E (float or numpy.ndarray): laboratory energy of nuclei in GeV
      corsika_id (int): corsika id of nucleus/mass group
      rel_delta (float): range of derivative relative to log10(E)
    Returns:
      (numpy.ndarray): spectral index :math:`\\gamma` of nuclei
    """
    delta = rel_delta * E
    return np.log10(
        self.nucleus_flux(corsika_id, E + delta)
        / self.nucleus_flux(corsika_id, E - delta)
    ) / np.log10((E + delta) / (E - delta))

delta_0(E)

Returns proton excess.

The proton excess is defined as :math:\delta_0 = \frac{\Phi_p - \Phi_n}{\Phi_p + \Phi_n}.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required

Returns: (numpy.ndarray): proton excess :math:\delta_0

Source code in crflux/models.py 
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
def delta_0(self, E):
    """Returns proton excess.

    The proton excess is defined as
    :math:`\\delta_0 = \\frac{\\Phi_p - \\Phi_n}{\\Phi_p + \\Phi_n}`.

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
    Returns:
      (numpy.ndarray): proton excess :math:`\\delta_0`
    """
    p_0 = 0.0
    n_0 = 0.0

    p_0 += self.nucleus_flux(14, E)

    p_0 += 2.0**2 * self.nucleus_flux(402, E * 4.0)
    n_0 += 2.0**2 * self.nucleus_flux(402, E * 4.0)

    p_0 += 6.0**2 * self.nucleus_flux(1206, E * 12.0)
    n_0 += 6.0**2 * self.nucleus_flux(1206, E * 12.0)

    p_0 += 14.0**2 * self.nucleus_flux(2814, E * 28.0)
    n_0 += 14.0**2 * self.nucleus_flux(2814, E * 28.0)

    a_fe = self.Z_A(5426)[1]
    p_0 += 26.0**2 * self.nucleus_flux(5426, E * a_fe)
    n_0 += 26.0**2 * self.nucleus_flux(5426, E * a_fe)

    return (p_0 - n_0) / (p_0 + n_0)

lnA(E)

Returns mean logarithmic mass .

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of particles in GeV

 required

Returns: (numpy.ndarray): mean (natural) logarithmic mass

Source code in crflux/models.py 
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
def lnA(self, E):
    """Returns mean logarithmic mass <ln A>.

    Args:
      E (float or numpy.ndarray): laboratory energy of particles in GeV
    Returns:
      (numpy.ndarray): mean (natural) logarithmic mass
    """
    sum_weight = 0.0

    for cid in self.nucleus_ids:
        if cid == 14:
            continue  # p has lnA = 0
        sum_weight += np.log(self.Z_A(cid)[1]) * self.nucleus_flux(cid, E)

    return sum_weight / self.total_flux(E)

Z_A(corsika_id)

Returns charge :math:Z and mass number :math:A corresponding to corsika_id

Parameters:

Name Type Description Default
corsika_id int

corsika id of nucleus/mass group

 required

Returns: (int,int): (Z, A) tuple

Source code in crflux/models.py 
274
275
276
277
278
279
280
281
282
283
284
285
286
287
def Z_A(self, corsika_id):
    """Returns charge :math:`Z` and mass number :math:`A` corresponding
    to ``corsika_id``

    Args:
      corsika_id (int): corsika id of nucleus/mass group
    Returns:
      (int,int): (Z, A) tuple
    """
    Z, A = 1, 1
    if corsika_id > 14:
        Z = corsika_id % 100
        A = (corsika_id - Z) // 100
    return Z, A

GlobalSplineFitReduced

Bases: GlobalSplineFit, PrimaryFlux

Reduced GSF proton/neutron model (globalsplinefit.ReducedGSF).

The nucleon flux of :class:GlobalSplineFit with its uncertainty represented by relative-flux deviations theta at pivot energies; theta[i] = reduction.sigma[i] is a one-sigma variation. Nucleon fluxes only: nucleus_ids is empty and nucleus_flux raises.

Use for reference: Fedynitch et al., arXiv:2609.32649 (2026).

Parameters

version : str globalsplinefit model version; default "2026.1". theta : array-like, optional Reduced-model parameter vector; the default is the central flux. reduction : globalsplinefit.ReducedGSF, optional Reuse an existing total p/n reduction across many parameter variations (time interval, cutoff and basis are then taken from it). time_interval : tuple or str, optional Solar-modulation interval, or "LIS". None uses Solar Cycle 24. geomagnetic_cutoff : float, optional Rigidity cutoff in GV, applied by globalsplinefit. basis : str Interpolation basis ("spline" or "hat") when building a reduction.

Source code in crflux/models.py 
 921
 922
 923
 924
 925
 926
 927
 928
 929
 930
 931
 932
 933
 934
 935
 936
 937
 938
 939
 940
 941
 942
 943
 944
 945
 946
 947
 948
 949
 950
 951
 952
 953
 954
 955
 956
 957
 958
 959
 960
 961
 962
 963
 964
 965
 966
 967
 968
 969
 970
 971
 972
 973
 974
 975
 976
 977
 978
 979
 980
 981
 982
 983
 984
 985
 986
 987
 988
 989
 990
 991
 992
 993
 994
 995
 996
 997
 998
 999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
class GlobalSplineFitReduced(GlobalSplineFit, PrimaryFlux):
    """Reduced GSF proton/neutron model (``globalsplinefit.ReducedGSF``).

    The nucleon flux of :class:`GlobalSplineFit` with its uncertainty
    represented by relative-flux deviations ``theta`` at pivot energies;
    ``theta[i] = reduction.sigma[i]`` is a one-sigma variation. Nucleon fluxes
    only: ``nucleus_ids`` is empty and ``nucleus_flux`` raises.

    Use for reference: Fedynitch et al., arXiv:2609.32649 (2026).

    Parameters
    ----------
    version : str
        globalsplinefit model version; default "2026.1".
    theta : array-like, optional
        Reduced-model parameter vector; the default is the central flux.
    reduction : globalsplinefit.ReducedGSF, optional
        Reuse an existing total p/n reduction across many parameter variations
        (time interval, cutoff and basis are then taken from it).
    time_interval : tuple or str, optional
        Solar-modulation interval, or "LIS". None uses Solar Cycle 24.
    geomagnetic_cutoff : float, optional
        Rigidity cutoff in GV, applied by globalsplinefit.
    basis : str
        Interpolation basis ("spline" or "hat") when building a reduction.
    """

    def __init__(
        self,
        version="2026.1",
        theta=None,
        *,
        reduction=None,
        time_interval=None,
        geomagnetic_cutoff=None,
        basis="spline",
    ):
        PrimaryFlux.__init__(self, geomagnetic_cutoff=geomagnetic_cutoff)
        from globalsplinefit import (
            GSFEnergyPerNucleon,
            GSFKineticEnergyPerNucleon,
            ReducedGSF,
        )

        if reduction is None:
            reduction = ReducedGSF(
                GSFEnergyPerNucleon(
                    version=version,
                    default_time_interval=time_interval,
                    default_rigidity_cutoff=geomagnetic_cutoff,
                ),
                basis=basis,
            )
        elif (
            time_interval is not None
            or geomagnetic_cutoff is not None
            or basis != "spline"
        ):
            raise ValueError(
                "Configure time_interval/geomagnetic_cutoff/basis on the supplied reduction"
            )
        if not isinstance(reduction, ReducedGSF) or reduction.per_group:
            raise ValueError("A total proton/neutron ReducedGSF is required")
        if not isinstance(reduction.model, GSFEnergyPerNucleon) or isinstance(
            reduction.model, GSFKineticEnergyPerNucleon
        ):
            raise ValueError(
                "MCEq requires total energy per nucleon, not kinetic energy"
            )
        self.reduction = reduction
        self.gsf = reduction.model
        self.nucleus_ids = []
        self.theta = (
            np.zeros(reduction.n_params)
            if theta is None
            else np.array(theta, dtype=float, copy=True)
        )
        if (
            self.theta.shape != (reduction.n_params,)
            or not np.isfinite(self.theta).all()
        ):
            raise ValueError(
                f"theta must contain {reduction.n_params} finite components"
            )
        self.theta.setflags(write=False)
        self.name = f"Global Spline Fit {reduction.model.version} reduced"
        self.sname = f"GSF{reduction.model.version}"

    def _nucleus_flux(self, corsika_id, E):
        raise NotImplementedError(
            "ReducedGSF provides only proton/neutron nucleon fluxes"
        )

    def p_and_n_flux(self, E):
        energy = np.atleast_1d(np.asarray(E, dtype=float))
        proton, neutron = self.reduction.flux(energy, self.theta)
        total = proton + neutron
        fraction = np.divide(proton, total, out=np.zeros_like(total), where=total > 0)
        return fraction, proton, neutron

total_flux(E)

Returns total flux of nuclei, the "all-particle-flux".

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of particles in GeV

 required

Returns: (numpy.ndarray): particle flux in :math:\Phi_{particles} in :math:(\text{m}^2 \text{s sr GeV})^{-1}

Source code in crflux/models.py 
107
108
109
110
111
112
113
114
115
116
def total_flux(self, E):
    """Returns total flux of nuclei, the "all-particle-flux".

    Args:
      E (float or numpy.ndarray): laboratory energy of particles in GeV
    Returns:
      (numpy.ndarray): particle flux in :math:`\\Phi_{particles}` in
      :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
    """
    return sum(self.nucleus_flux(corsika_id, E) for corsika_id in self.nucleus_ids)

nucleon_gamma(E, rel_delta=0.01)

Returns differential spectral index of all-nucleon-flux obtained from a numerical derivative.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required
rel_delta float

range of derivative relative to log10(E)

 0.01

Returns: (numpy.ndarray): spectral index :math:\gamma of nucleons

Source code in crflux/models.py 
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
def nucleon_gamma(self, E, rel_delta=0.01):
    """Returns differential spectral index of all-nucleon-flux
    obtained from a numerical derivative.

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
      rel_delta (float): range of derivative relative to log10(E)
    Returns:
      (numpy.ndarray): spectral index :math:`\\gamma` of nucleons
    """
    delta = rel_delta * E
    fl = self.tot_nucleon_flux
    return np.log10(fl(E + delta) / fl(E - delta)) / np.log10(
        (E + delta) / (E - delta)
    )

nucleus_gamma(E, corsika_id, rel_delta=0.01)

Returns differential spectral index of nuclei obtained from a numerical derivative.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nuclei in GeV

 required
corsika_id int

corsika id of nucleus/mass group

 required
rel_delta float

range of derivative relative to log10(E)

 0.01

Returns: (numpy.ndarray): spectral index :math:\gamma of nuclei

Source code in crflux/models.py 
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
def nucleus_gamma(self, E, corsika_id, rel_delta=0.01):
    """Returns differential spectral index of nuclei
    obtained from a numerical derivative.

    Args:
      E (float or numpy.ndarray): laboratory energy of nuclei in GeV
      corsika_id (int): corsika id of nucleus/mass group
      rel_delta (float): range of derivative relative to log10(E)
    Returns:
      (numpy.ndarray): spectral index :math:`\\gamma` of nuclei
    """
    delta = rel_delta * E
    return np.log10(
        self.nucleus_flux(corsika_id, E + delta)
        / self.nucleus_flux(corsika_id, E - delta)
    ) / np.log10((E + delta) / (E - delta))

delta_0(E)

Returns proton excess.

The proton excess is defined as :math:\delta_0 = \frac{\Phi_p - \Phi_n}{\Phi_p + \Phi_n}.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required

Returns: (numpy.ndarray): proton excess :math:\delta_0

Source code in crflux/models.py 
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
def delta_0(self, E):
    """Returns proton excess.

    The proton excess is defined as
    :math:`\\delta_0 = \\frac{\\Phi_p - \\Phi_n}{\\Phi_p + \\Phi_n}`.

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
    Returns:
      (numpy.ndarray): proton excess :math:`\\delta_0`
    """
    p_0 = 0.0
    n_0 = 0.0

    p_0 += self.nucleus_flux(14, E)

    p_0 += 2.0**2 * self.nucleus_flux(402, E * 4.0)
    n_0 += 2.0**2 * self.nucleus_flux(402, E * 4.0)

    p_0 += 6.0**2 * self.nucleus_flux(1206, E * 12.0)
    n_0 += 6.0**2 * self.nucleus_flux(1206, E * 12.0)

    p_0 += 14.0**2 * self.nucleus_flux(2814, E * 28.0)
    n_0 += 14.0**2 * self.nucleus_flux(2814, E * 28.0)

    a_fe = self.Z_A(5426)[1]
    p_0 += 26.0**2 * self.nucleus_flux(5426, E * a_fe)
    n_0 += 26.0**2 * self.nucleus_flux(5426, E * a_fe)

    return (p_0 - n_0) / (p_0 + n_0)

lnA(E)

Returns mean logarithmic mass .

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of particles in GeV

 required

Returns: (numpy.ndarray): mean (natural) logarithmic mass

Source code in crflux/models.py 
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
def lnA(self, E):
    """Returns mean logarithmic mass <ln A>.

    Args:
      E (float or numpy.ndarray): laboratory energy of particles in GeV
    Returns:
      (numpy.ndarray): mean (natural) logarithmic mass
    """
    sum_weight = 0.0

    for cid in self.nucleus_ids:
        if cid == 14:
            continue  # p has lnA = 0
        sum_weight += np.log(self.Z_A(cid)[1]) * self.nucleus_flux(cid, E)

    return sum_weight / self.total_flux(E)

Z_A(corsika_id)

Returns charge :math:Z and mass number :math:A corresponding to corsika_id

Parameters:

Name Type Description Default
corsika_id int

corsika id of nucleus/mass group

 required

Returns: (int,int): (Z, A) tuple

Source code in crflux/models.py 
274
275
276
277
278
279
280
281
282
283
284
285
286
287
def Z_A(self, corsika_id):
    """Returns charge :math:`Z` and mass number :math:`A` corresponding
    to ``corsika_id``

    Args:
      corsika_id (int): corsika id of nucleus/mass group
    Returns:
      (int,int): (Z, A) tuple
    """
    Z, A = 1, 1
    if corsika_id > 14:
        Z = corsika_id % 100
        A = (corsika_id - Z) // 100
    return Z, A

GlobalSplineFitBeta

Bases: PrimaryFlux

Data-driven fit of direct and indirect measurements of the cosmic ray flux and composition. Tracks the mass composition using four leading elements (p, He, O, Fe), whose flux is modeled by shaped spline functions. Assumes fixed flux ratios in rigidity for the flux of subleading elements. Covers the whole rigidity range from 10 GV to 10^11 GeV.

.. deprecated:: Tabulated nucleon flux of GSF 2017 at Earth, modulated with a fixed force-field potential of 554 MV (Usoskin, July 1998). It equals globalsplinefit GSFEnergyPerNucleon(version="2017") with time_interval=(199807, 199808) (protons to 0.2%; neutrons differ by a few percent from the nucleon-count convention). Use GlobalSplineFit(version="2017") (or "2019", or a current version) instead. The class picks the most recent spline file in the package directory.

Use for reference: Dembinski et al., PoS ICRC2017 533 https://inspirehep.net/literature/1639832

Source code in crflux/models.py 
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
class GlobalSplineFitBeta(PrimaryFlux):
    """Data-driven fit of direct and indirect measurements of the
    cosmic ray flux and composition. Tracks the mass composition using
    four leading elements (p, He, O, Fe), whose flux is modeled by shaped
    spline functions. Assumes fixed flux ratios in rigidity for the flux
    of subleading elements. Covers the whole rigidity range from
    10 GV to 10^11 GeV.

    .. deprecated::
        Tabulated nucleon flux of GSF 2017 at Earth, modulated with a fixed
        force-field potential of 554 MV (Usoskin, July 1998). It equals
        globalsplinefit ``GSFEnergyPerNucleon(version="2017")`` with
        ``time_interval=(199807, 199808)`` (protons to 0.2%; neutrons differ by
        a few percent from the nucleon-count convention). Use
        ``GlobalSplineFit(version="2017")`` (or "2019", or a current version)
        instead. The class picks the most recent spline file in the package
        directory.

    Use for reference: Dembinski et al., PoS ICRC2017 533
    https://inspirehep.net/literature/1639832
    """

    name = "Global Spline Fit"
    sname = "GSF"

    def __init__(self, spl_fname=None, **kwargs):
        warnings.warn(
            "GlobalSplineFitBeta is deprecated: it tabulates the GSF 2017 "
            "nucleon flux at phi = 554 MV. Use globalsplinefit "
            "crflux.models.GlobalSplineFit(version='2017') instead (or '2019', "
            "or a current version).",
            DeprecationWarning,
            stacklevel=2,
        )
        super().__init__(**kwargs)
        import bz2
        import os
        import os.path as path
        import pickle

        base_path = path.dirname(path.abspath(__file__))

        if spl_fname is None:
            # Look for all files starting with GSF_spline
            gsf_files = [
                fname
                for fname in os.listdir(base_path)
                if fname.startswith("GSF_spline_")
            ]

            if not gsf_files:
                raise Exception(
                    self.__class__.__name__ + ": No spline files found in " + base_path
                )

            spl_fname = gsf_files[0]

            # Find out file datetag
            def fdate(fn):
                return int(os.path.splitext(os.path.splitext(fn)[0])[0].split("_")[-1])

            # Pick the latest
            for fn in gsf_files:
                if fdate(fn) >= fdate(spl_fname):
                    spl_fname = fn
            spl_fname = os.path.join(base_path, spl_fname)

        self.p_frac_spl, self.p_flux_spl, self.n_flux_spl = pickle.load(
            bz2.BZ2File(spl_fname), encoding="latin1"
        )

        self.nucleus_ids = []

    def p_and_n_flux(self, E):
        """Returns tuple with proton fraction, proton flux and neutron flux.

        The proton fraction is defined as
        :math:`\\frac{\\Phi_p}{\\Phi_p + \\Phi_n}`.

        Args:
          E (float or numpy.ndarray): laboratory energy of nucleons in GeV
        Returns:
          (float,float,float): proton fraction, proton flux, neutron flux
        """
        E = np.atleast_1d(np.asarray(E, dtype=float))
        p_frac = self.p_frac_spl(np.log(E))
        p_flux = np.exp(self.p_flux_spl(np.log(E)))
        n_flux = np.exp(self.n_flux_spl(np.log(E)))

        if self.geomagnetic_cutoff is not None:
            # For nucleon fluxes, apply proton cutoff (Z=1, A=1)
            E_cut = self.geomagnetic_cutoff
            mask = E >= E_cut
            p_flux = np.where(mask, p_flux, 0.0)
            n_flux = np.where(mask, n_flux, 0.0)

        return p_frac, p_flux, n_flux

    def tot_nucleon_flux(self, E):
        """Returns total flux of nucleons, the "all-nucleon-flux".

        Args:
          E (float or numpy.ndarray): laboratory energy of nucleons in GeV
        Returns:
          (numpy.ndarray): nucleon flux :math:`\\Phi_{nucleons}` in
          :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
        """
        return np.sum(self.p_and_n_flux(E)[1:], axis=0)

    def _nucleus_flux(self, corsika_id, E):
        """Dummy function, since particle fluxes are not supported
        in the spline interface version."""
        return np.zeros_like(E)

p_and_n_flux(E)

Returns tuple with proton fraction, proton flux and neutron flux.

The proton fraction is defined as :math:\frac{\Phi_p}{\Phi_p + \Phi_n}.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required

Returns: (float,float,float): proton fraction, proton flux, neutron flux

Source code in crflux/models.py 
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
def p_and_n_flux(self, E):
    """Returns tuple with proton fraction, proton flux and neutron flux.

    The proton fraction is defined as
    :math:`\\frac{\\Phi_p}{\\Phi_p + \\Phi_n}`.

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
    Returns:
      (float,float,float): proton fraction, proton flux, neutron flux
    """
    E = np.atleast_1d(np.asarray(E, dtype=float))
    p_frac = self.p_frac_spl(np.log(E))
    p_flux = np.exp(self.p_flux_spl(np.log(E)))
    n_flux = np.exp(self.n_flux_spl(np.log(E)))

    if self.geomagnetic_cutoff is not None:
        # For nucleon fluxes, apply proton cutoff (Z=1, A=1)
        E_cut = self.geomagnetic_cutoff
        mask = E >= E_cut
        p_flux = np.where(mask, p_flux, 0.0)
        n_flux = np.where(mask, n_flux, 0.0)

    return p_frac, p_flux, n_flux

tot_nucleon_flux(E)

Returns total flux of nucleons, the "all-nucleon-flux".

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required

Returns: (numpy.ndarray): nucleon flux :math:\Phi_{nucleons} in :math:(\text{m}^2 \text{s sr GeV})^{-1}

Source code in crflux/models.py 
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
def tot_nucleon_flux(self, E):
    """Returns total flux of nucleons, the "all-nucleon-flux".

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
    Returns:
      (numpy.ndarray): nucleon flux :math:`\\Phi_{nucleons}` in
      :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
    """
    return np.sum(self.p_and_n_flux(E)[1:], axis=0)

nucleus_flux(corsika_id, E)

Returns the flux of nuclei corresponding to the corsika_id at energy E, with optional geomagnetic cutoff.

Parameters:

Name Type Description Default
corsika_id int

see :mod:crflux for description.

 required
E float or ndarray

laboratory energy of nucleus in GeV

 required

Returns: (numpy.ndarray): flux of single nucleus type :math:\Phi_{nucleus} in :math:(\text{m}^2 \text{s sr GeV})^{-1}

Source code in crflux/models.py 
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
def nucleus_flux(self, corsika_id, E):
    """Returns the flux of nuclei corresponding to
    the ``corsika_id`` at energy ``E``, with optional geomagnetic cutoff.

    Args:
      corsika_id (int): see :mod:`crflux` for description.
      E (float or numpy.ndarray): laboratory energy of nucleus in GeV
    Returns:
      (numpy.ndarray): flux of single nucleus type :math:`\\Phi_{nucleus}`
      in :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
    """
    E = np.atleast_1d(np.asarray(E, dtype=float))
    flux = self._nucleus_flux(corsika_id, E)
    if self.geomagnetic_cutoff is not None:
        Z, A = self.Z_A(corsika_id)
        E_cut = Z * self.geomagnetic_cutoff
        flux = np.where(E >= E_cut, flux, 0.0)
    return flux

total_flux(E)

Returns total flux of nuclei, the "all-particle-flux".

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of particles in GeV

 required

Returns: (numpy.ndarray): particle flux in :math:\Phi_{particles} in :math:(\text{m}^2 \text{s sr GeV})^{-1}

Source code in crflux/models.py 
107
108
109
110
111
112
113
114
115
116
def total_flux(self, E):
    """Returns total flux of nuclei, the "all-particle-flux".

    Args:
      E (float or numpy.ndarray): laboratory energy of particles in GeV
    Returns:
      (numpy.ndarray): particle flux in :math:`\\Phi_{particles}` in
      :math:`(\\text{m}^2 \\text{s sr GeV})^{-1}`
    """
    return sum(self.nucleus_flux(corsika_id, E) for corsika_id in self.nucleus_ids)

nucleon_gamma(E, rel_delta=0.01)

Returns differential spectral index of all-nucleon-flux obtained from a numerical derivative.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required
rel_delta float

range of derivative relative to log10(E)

 0.01

Returns: (numpy.ndarray): spectral index :math:\gamma of nucleons

Source code in crflux/models.py 
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
def nucleon_gamma(self, E, rel_delta=0.01):
    """Returns differential spectral index of all-nucleon-flux
    obtained from a numerical derivative.

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
      rel_delta (float): range of derivative relative to log10(E)
    Returns:
      (numpy.ndarray): spectral index :math:`\\gamma` of nucleons
    """
    delta = rel_delta * E
    fl = self.tot_nucleon_flux
    return np.log10(fl(E + delta) / fl(E - delta)) / np.log10(
        (E + delta) / (E - delta)
    )

nucleus_gamma(E, corsika_id, rel_delta=0.01)

Returns differential spectral index of nuclei obtained from a numerical derivative.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nuclei in GeV

 required
corsika_id int

corsika id of nucleus/mass group

 required
rel_delta float

range of derivative relative to log10(E)

 0.01

Returns: (numpy.ndarray): spectral index :math:\gamma of nuclei

Source code in crflux/models.py 
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
def nucleus_gamma(self, E, corsika_id, rel_delta=0.01):
    """Returns differential spectral index of nuclei
    obtained from a numerical derivative.

    Args:
      E (float or numpy.ndarray): laboratory energy of nuclei in GeV
      corsika_id (int): corsika id of nucleus/mass group
      rel_delta (float): range of derivative relative to log10(E)
    Returns:
      (numpy.ndarray): spectral index :math:`\\gamma` of nuclei
    """
    delta = rel_delta * E
    return np.log10(
        self.nucleus_flux(corsika_id, E + delta)
        / self.nucleus_flux(corsika_id, E - delta)
    ) / np.log10((E + delta) / (E - delta))

delta_0(E)

Returns proton excess.

The proton excess is defined as :math:\delta_0 = \frac{\Phi_p - \Phi_n}{\Phi_p + \Phi_n}.

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of nucleons in GeV

 required

Returns: (numpy.ndarray): proton excess :math:\delta_0

Source code in crflux/models.py 
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
def delta_0(self, E):
    """Returns proton excess.

    The proton excess is defined as
    :math:`\\delta_0 = \\frac{\\Phi_p - \\Phi_n}{\\Phi_p + \\Phi_n}`.

    Args:
      E (float or numpy.ndarray): laboratory energy of nucleons in GeV
    Returns:
      (numpy.ndarray): proton excess :math:`\\delta_0`
    """
    p_0 = 0.0
    n_0 = 0.0

    p_0 += self.nucleus_flux(14, E)

    p_0 += 2.0**2 * self.nucleus_flux(402, E * 4.0)
    n_0 += 2.0**2 * self.nucleus_flux(402, E * 4.0)

    p_0 += 6.0**2 * self.nucleus_flux(1206, E * 12.0)
    n_0 += 6.0**2 * self.nucleus_flux(1206, E * 12.0)

    p_0 += 14.0**2 * self.nucleus_flux(2814, E * 28.0)
    n_0 += 14.0**2 * self.nucleus_flux(2814, E * 28.0)

    a_fe = self.Z_A(5426)[1]
    p_0 += 26.0**2 * self.nucleus_flux(5426, E * a_fe)
    n_0 += 26.0**2 * self.nucleus_flux(5426, E * a_fe)

    return (p_0 - n_0) / (p_0 + n_0)

lnA(E)

Returns mean logarithmic mass .

Parameters:

Name Type Description Default
E float or ndarray

laboratory energy of particles in GeV

 required

Returns: (numpy.ndarray): mean (natural) logarithmic mass

Source code in crflux/models.py 
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
def lnA(self, E):
    """Returns mean logarithmic mass <ln A>.

    Args:
      E (float or numpy.ndarray): laboratory energy of particles in GeV
    Returns:
      (numpy.ndarray): mean (natural) logarithmic mass
    """
    sum_weight = 0.0

    for cid in self.nucleus_ids:
        if cid == 14:
            continue  # p has lnA = 0
        sum_weight += np.log(self.Z_A(cid)[1]) * self.nucleus_flux(cid, E)

    return sum_weight / self.total_flux(E)

Z_A(corsika_id)

Returns charge :math:Z and mass number :math:A corresponding to corsika_id

Parameters:

Name Type Description Default
corsika_id int

corsika id of nucleus/mass group

 required

Returns: (int,int): (Z, A) tuple

Source code in crflux/models.py 
274
275
276
277
278
279
280
281
282
283
284
285
286
287
def Z_A(self, corsika_id):
    """Returns charge :math:`Z` and mass number :math:`A` corresponding
    to ``corsika_id``

    Args:
      corsika_id (int): corsika id of nucleus/mass group
    Returns:
      (int,int): (Z, A) tuple
    """
    Z, A = 1, 1
    if corsika_id > 14:
        Z = corsika_id % 100
        A = (corsika_id - Z) // 100
    return Z, A