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 
 58
 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
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 
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
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 
106
107
108
109
110
111
112
113
114
115
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 
117
118
119
120
121
122
123
124
125
126
127
128
129
130
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 
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
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 
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
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 
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
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 
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
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 
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
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 
273
274
275
276
277
278
279
280
281
282
283
284
285
286
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 
289
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
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 
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
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 
106
107
108
109
110
111
112
113
114
115
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 
117
118
119
120
121
122
123
124
125
126
127
128
129
130
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 
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
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 
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
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 
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
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 
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
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 
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
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 
273
274
275
276
277
278
279
280
281
282
283
284
285
286
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 
381
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
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 
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
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 
106
107
108
109
110
111
112
113
114
115
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 
117
118
119
120
121
122
123
124
125
126
127
128
129
130
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 
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
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 
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
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 
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
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 
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
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 
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
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 
273
274
275
276
277
278
279
280
281
282
283
284
285
286
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 
450
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
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 
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
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 
106
107
108
109
110
111
112
113
114
115
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 
117
118
119
120
121
122
123
124
125
126
127
128
129
130
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 
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
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 
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
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 
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
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 
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
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 
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
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 
273
274
275
276
277
278
279
280
281
282
283
284
285
286
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 
484
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
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 
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
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 
106
107
108
109
110
111
112
113
114
115
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 
117
118
119
120
121
122
123
124
125
126
127
128
129
130
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 
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
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 
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
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 
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
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 
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
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 
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
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 
273
274
275
276
277
278
279
280
281
282
283
284
285
286
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 
590
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
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 
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
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 
106
107
108
109
110
111
112
113
114
115
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 
117
118
119
120
121
122
123
124
125
126
127
128
129
130
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 
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
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 
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
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 
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
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 
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
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 
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
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 
273
274
275
276
277
278
279
280
281
282
283
284
285
286
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 
638
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
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 
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
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 
106
107
108
109
110
111
112
113
114
115
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 
117
118
119
120
121
122
123
124
125
126
127
128
129
130
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 
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
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 
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
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 
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
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 
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
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 
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
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 
273
274
275
276
277
278
279
280
281
282
283
284
285
286
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 
756
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
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 
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
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 
106
107
108
109
110
111
112
113
114
115
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 
117
118
119
120
121
122
123
124
125
126
127
128
129
130
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 
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
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 
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
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 
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
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 
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
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 
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
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 
273
274
275
276
277
278
279
280
281
282
283
284
285
286
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 
789
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
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 
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
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 
106
107
108
109
110
111
112
113
114
115
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 
117
118
119
120
121
122
123
124
125
126
127
128
129
130
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 
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
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 
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
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 
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
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 
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
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 
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
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 
273
274
275
276
277
278
279
280
281
282
283
284
285
286
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 
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
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 
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
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 
106
107
108
109
110
111
112
113
114
115
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 
117
118
119
120
121
122
123
124
125
126
127
128
129
130
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 
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
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 
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
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 
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
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 
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
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 
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
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 
273
274
275
276
277
278
279
280
281
282
283
284
285
286
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.

Tabulated ICRC 2017 version. The full interface will become available, when the model is published. The class picks the most recent spline file in the current directory.

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

Source code in crflux/models.py 
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
917
918
919
920
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
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.

    Tabulated ICRC 2017 version. The full interface will become
    available, when the model is published. The class picks the most recent
    spline file in the current 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):
        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 
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
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 
935
936
937
938
939
940
941
942
943
944
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 
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
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 
106
107
108
109
110
111
112
113
114
115
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 
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
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 
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
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 
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
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 
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
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 
273
274
275
276
277
278
279
280
281
282
283
284
285
286
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