Source code for cripticpy.sr

"""
Special relativistic conversion functions
"""

import numpy as np
try:
    import astropy.units as u
    import astropy.constants as const
    from astropy.constants import m_p, m_e, c
except:
    u = None

# Particle types used in criptic; note masses are in units of m_p,
# charges are in units of e, as used internally in criptic
parttypes = [
    { 'name' : 'proton',
      'symbol' : r'$p$',
      'm' : 1.0,
      'Z' : 1 },
    { 'name' : 'electron',
      'symbol' : r'$e^-$',
      'm' : 1.0 / 1836.15267343,  # CODATA 2018
      'Z' : -1 },
    { 'name' : 'positron',
      'symbol' : r'$e^+$',
      'm' : 1.0 / 1836.15267343,  # CODATA 2018
      'Z' : 1 },
    { 'name': 'boron10',
      'symbol': r'${^{10}\mathrm{B}}$',
      'm': 10.0129370 / 1.00727647,  # Nuclear Data Center at KAERI, CODATA
      'Z': 5 },   
    { 'name': 'boron11',
      'symbol': r'${^{11}\mathrm{B}}$',
      'm': 11.0093055 / 1.00727647,  # Nuclear Data Center at KAERI, CODATA
      'Z': 5 },
    { 'name': 'carbon12',
      'symbol': r'${^{12}\mathrm{C}}$',
      'm': 12.0000000 / 1.00727647,  # Nuclear Data Center at KAERI, CODATA
      'Z': 6 },
    { 'name': 'nitrogen14',
      'symbol': r'${^{14}\mathrm{N}}$',
      'm': 14.0030740 / 1.00727647,  # Nuclear Data Center at KAERI, CODATA
      'Z': 7 }, 
    { 'name': 'oxygen16',
      'symbol': r'${^{16}\mathrm{O}}$',
      'm': 15.9949146 / 1.00727647,  # Nuclear Data Center at KAERI, CODATA
      'Z': 8 }
    ]

def _Tp(p, ptype):
    """
    Compute a kinetic energy from a momentum for a given particle
    type

    Parameters
       p : array
          momentum in units where m_p = 1 and c = 1
       ptype : array
          array of particle types

    Returns
       T : array
          kinetic energies in units where m_p = 1 and c = 1
    """
    m = np.array([ parttypes[i]['m'] for i in ptype ])
    pm = p / m
    T = m * (np.sqrt(1+pm*pm) - 1)
    idx = pm < 1e-6
    T[idx] = p[idx]**2 / (2.0*m[idx]) * (1 - pm[idx]**2/4.0)
    return T

def _pT(T, ptype):
    """
    Compute a momentum from a kinetic energy for a given particle
    type

    Parameters
       T : array
          kinetic energies in units where m_p = 1 and c = 1
       ptype : array
          array of particle types

    Returns
       p : array
          momentum in units where m_p = 1 and c = 1
    """
    m = np.array([ parttypes[i]['m'] for i in ptype ])
    Tm = T / m
    p = m * np.sqrt(Tm * (2 + Tm));
    idx = Tm < 1e-6
    p[idx] = np.sqrt(Tm[idx]) + m * Tm[idx]/2 * np.sqrt(Tm[idx]/2)
    return p

def _vp(p, ptype):
    """
    Compute a velocity from a kinetic energy for a given particle type

    Parameters
       p : array
          momentum in units where m_p = 1 and c = 1
       ptype : array
          array of particle types

    Returns
       v : array
          velocity in units where c = 1
    """
    m = np.array([ parttypes[i]['m'] for i in ptype ])
    pm = p / m
    v = 1.0 / np.sqrt(1 + 1/pm**2)
    idx = pm < 1e-6
    v[idx] = pm[idx] - 0.5 / pm[idx]**3
    return v

def _dTdp(p, ptype):
    """
    Compute derivative of kinetic energy with respect to momentum for
    a given particle type

    Parameters:
       p : array
          momentum in units where m_p = 1 and c = 1
       ptype : array
          array of particle types

    Returns
       dTdp : array
          dT/dp in units where m_p = 1 and c = 1
    """
    m = np.array([ parttypes[i]['m'] for i in ptype ])
    pm = p / m
    dTdp = pm / np.sqrt(1 + pm**2)
    idx = pm < 1e-6
    dTdp[idx] = pm[idx] * (1 - pm[idx]**2/2)
    return dTdp


def _ptype_convert(ptype):
    """
    Convert a particle type in human-readable input to the internal
    numbering scheme used in criptic

    Parameters
       ptype : int, string, or array of int or string
          particle type(s); can be specified by number as used
          internally in criptic, or by name; known names are 'proton',
          'electron', 'positron', 'boron10', 'boron11', 'carbon12', 'nitrogen14', 'oxygen16'

    Returns
       ptarr : array
          array of integers giving particle type in criptic
          representation
    """
    if type(ptype) is str:
        if ptype == 'proton':
            ptype_ = np.array([0])
        elif ptype == 'electron':
            ptype_ = np.array([1])
        elif ptype == 'positron':
            ptype_ = np.array([2])
        elif ptype == 'boron10':
            ptype_ = np.array([3])
        elif ptype == 'boron11':
            ptype_ = np.array([4])
        elif ptype == 'carbon12':   
            ptype_ = np.array([5])
        elif ptype == 'nitrogen14':   
            ptype_ = np.array([6])
        elif ptype == 'oxygen16':
            ptype_ = np.array([7])
        else:
            raise ValueError("unknown particle type" + ptype)
    elif hasattr(ptype, '__iter__'):
        ptype_ = []
        for p in ptype:
            if type(p) is str:
                if p == 'proton':
                    ptype_.append(0)
                elif p == 'electron':
                    ptype_.append(1)
                elif p == 'positron':
                    ptype_.append(2)
                elif p == 'boron10':
                    ptype_.append(3)
                elif p == 'boron11':
                    ptype_.append(4)
                elif p == 'carbon12':
                    ptype_.append(5)
                elif p == 'nitrogen14':
                    ptype_.append(6)
                elif p == 'oxygen16':
                    ptype_.append(7)
                else:
                    raise ValueError("unknown particle type" + p)
            else:
                ptype_.append(p)
        ptype_ = np.array(ptype)
    else:
        ptype_ = np.array([ptype])
    return ptype_

        
[docs] def Tp(p, ptype): """ Compute a kinetic energy from a momentum for a given particle type Parameters p : astropy.Quantity momentum ptype : int, string, or array of int or string particle type(s); can be specified by number as used internally in criptic, or by name; known names are 'proton', 'electron', 'positron', 'boron10', 'boron11', 'carbon12', 'nitrogen14', 'oxygen16' Returns T : astropy.Quantity kinetic energies """ pdim = np.atleast_1d((p / (m_p * c)).to('')) pt = _ptype_convert(ptype) Tdim = _Tp(pdim, pt) T = Tdim * m_p * c**2 try: p[0] except: T = T[0] return T
[docs] def pT(T, ptype): """ Compute a momentum from a kinetic energy for a given particle type Parameters T : astropy.Quantity kinetic energy ptype : int, string, or array of int or string particle type(s); can be specified by number as used internally in criptic, or by name; known names are 'proton', 'electron', 'positron', 'boron10', 'boron11', 'carbon12', 'nitrogen14', 'oxygen16' Returns p : astropy.Quantity momentum """ Tdim = np.atleast_1d((T / (m_p * c**2)).to('')) pt = _ptype_convert(ptype) pdim = _pT(Tdim, pt) p = pdim * m_p * c try: T[0] except: p = p[0] return p
[docs] def vp(p, ptype): """ Compute a velocity from a momentum for a given particle type Parameters p : astropy.Quantity momentum ptype : int, string, or array of int or string particle type(s); can be specified by number as used internally in criptic, or by name; known names are 'proton', 'electron', 'positron', 'boron10', 'boron11', 'carbon12', 'nitrogen14', 'oxygen16' Returns v : astropy.Quantity velocities """ pdim = np.atleast_1d((p / (m_p * c)).to('')) pt = _ptype_convert(ptype) vdim = _vp(pdim, pt) v = vdim * c try: p[0] except: v = v[0] return v
def dTdp(p, ptype): """ Compute derivative of kinetic energy with respect to momentum for a given particle type Parameters p : astropy.Quantity momentum ptype : int, string, or array of int or string particle type(s); can be specified by number as used internally in criptic, or by name; known names are 'proton', 'electron', 'positron', 'boron10', 'boron11', 'carbon12', 'nitrogen14', 'oxygen16' Returns v : astropy.Quantity derivative of kinetic energy with respect to momentum """ pdim = np.atleast_1d((p / (m_p * c)).to('')) pt = _ptype_convert(ptype) dTdpdim = _dTdp(pdim, pt) dTdp = dTdpdim * c try: p[0] except: dTdp = dTdp[0] return dTdp