Source code for gammapy.spectrum.cosmic_ray

# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""Simple models for cosmic ray spectra at Earth.

For measurements, the "Database of Charged Cosmic Rays (CRDB)" is a great resource:
http://lpsc.in2p3.fr/cosmic-rays-db/
"""
from __future__ import absolute_import, division, print_function, unicode_literals
from collections import OrderedDict
import numpy as np
from astropy.units import Quantity

__all__ = ['cosmic_ray_flux']


def _power_law(E, N, k):
    E = Quantity(E, 'TeV')
    E0 = Quantity(1, 'TeV')
    N = Quantity(N, 'm^-2 s^-1 TeV^-1 sr^-1')
    flux = N * (E / E0) ** (-k)
    return flux


def _log_normal(E, L, E_p, w):
    E = Quantity(E, 'TeV')
    E_p = Quantity(E_p, 'TeV')
    L = Quantity(L, 'm^-2 s^-1 sr^-1')
    term1 = L / (E * w * np.sqrt(2 * np.pi))
    term2 = np.exp(-np.log(E / E_p) ** 2 / (2 * w ** 2))
    return term1 * term2


def _electron_spectrum(E, N, k, L, E_p, w):
    flux = _power_law(E, N, k)
    flux += _log_normal(E, L, E_p, w)
    return flux


[docs]def cosmic_ray_flux(energy, particle='proton'): """Cosmic ray flux at Earth. These are the spectra assumed in this CTA study: Table 3 in http://adsabs.harvard.edu/abs/2013APh....43..171B The hadronic spectra are simple power-laws, the electron spectrum is the sum of a power law and a log-normal component to model the "Fermi shoulder". Parameters ---------- energy : `~astropy.units.Quantity` Particle energy particle : {'electron', 'proton', 'He', 'N', 'Si', 'Fe'} Particle type Returns ------- flux : `~astropy.units.Quantity` Cosmic ray flux in unit ``m^-2 s^-1 TeV^-1 sr^-1`` """ pars = OrderedDict() pars['electron'] = OrderedDict(N=6.85e-5, k=3.21, L=3.19e-3, E_p=0.107, w=0.776) pars['proton'] = OrderedDict(N=0.096, k=2.70) pars['N'] = OrderedDict(N=0.0719, k=2.64) pars['Si'] = OrderedDict(N=0.0284, k=2.66) pars['Fe'] = OrderedDict(N=0.0134, k=2.63) if particle == 'electron': return _electron_spectrum(energy, **pars['electron']) elif particle in ['proton', 'He', 'N', 'Si', 'Fe']: return _power_law(energy, **pars[particle]) else: raise ValueError('Invalid argument for particle: {0}'.format(particle))