magnetic_dist.py

# -*- coding: utf-8 -*-

from naima.models import Synchrotron
from naima.utils import trapz_loglog
from naima.extern.validator import validate_scalar
from naima.radiative import _validate_ene

##use naima.radiative logger
#import naima
#log = naima.radiative.log
#log.setLevel(1) # set to DEBUG

class log(object):
    def debug(msg):
        #print('DEBUG: '+msg)
        pass

import numpy as np
import astropy.units as u
from astropy.constants import c, G, m_e, h, hbar, k_B, R_sun, sigma_sb, e, m_p, M_sun
e = e.gauss

mec2 = (m_e * c ** 2).cgs
mec2_unit = u.Unit(mec2)

ar = (4 * sigma_sb / c).to('erg/(cm3 K4)')

heaviside = lambda x: (np.sign(x) + 1) / 2.


class PowerLawB(object):
    """Magnetic field powerlaw distribution

    Parameters
    ----------
    index : float
        Index of the powerlaw
    Bmin : `~astropy.units.Quantity`
        Minimum magnetic field strength
    Bmax : `~astropy.units.Quantity`
        Maximum magnetic field strength
    """
    def __init__(self, index, Bmin=0.1*u.uG, Bmax=300*u.uG):
        self.index = validate_scalar('index',index)
        self.Bmin = validate_scalar('Bmin',Bmin,physical_type='magnetic flux density').to('G')
        self.Bmax = validate_scalar('Bmax',Bmax,physical_type='magnetic flux density').to('G')

    def pdf(self,B):
        g = self.index
        x = B.to('G').value
        a = self.Bmin.to('G').value
        b = self.Bmax.to('G').value
        norm=((1 - g) / (b ** (1 - g) - a ** (1 - g)))
        pdf = norm * np.where((x >= a) * (x <= b), x ** -g, np.zeros_like(x))
        return pdf

    def rms(self):
        g = self.index
        a = self.Bmin.to('G').value
        b = self.Bmax.to('G').value
        rms = np.sqrt((1 - g) / (3 - g) * (b ** (3 - g) - a ** (3 - g)) /
                      (b ** (1 - g) - a ** (1 - g)))
        return rms * u.G


class SynchrotronBdist(Synchrotron):
    """Synchrotron emission from an electron population in a turbulent magnetic field.

    Parameters
    ----------
    particle_distribution : function
        Particle distribution function, taking electron energies as a
        `~astropy.units.Quantity` array or float, and returning the particle
        energy density in units of number of electrons per unit energy as a
        `~astropy.units.Quantity` array or float.

    B_distribution : class instance
        Distribution of the magnetic field strength.
        Must have ``Bmin``, ``Bmax``, ``rms``, and ``pdf`` attributes.  ``pdf``
        will take magnetic field strength values and return the PDF for that
        value.

    Other parameters
    ----------------
    log10gmin : float
        Base 10 logarithm of the minimum Lorentz factor for the electron
        distribution. Default is 4 (:math:`E_e ≈ 5` GeV).

    log10gmax : float
        Base 10 logarithm of the maximum Lorentz factor for the electron
        distribution. Default is 9 (:math:`E_e ≈ 510` TeV).

    ngamd : scalar
        Number of points per decade in energy for the electron energy and
        distribution arrays. Default is 100.
    """
    def __init__(self, particle_distribution, B_distribution, Bmin=0.1*u.uG,
                 Bmax=300*u.uG, **kwargs):
        self._memoize=True
        self._cache = {}
        self._queue = []
        self.particle_distribution = particle_distribution
        # check that the particle distribution returns particles per unit energy
        P = self.particle_distribution(1*u.TeV)
        validate_scalar('particle distribution', P, physical_type='differential energy')
        validate_scalar('Bmin', Bmin, physical_type='magnetic flux density')
        validate_scalar('Bmax', Bmax, physical_type='magnetic flux density')
        self.B_distribution = B_distribution
        self.Eemin = 1 * u.GeV
        self.Eemax = 1e9 * mec2
        self.nEed = 100
        self.__dict__.update(**kwargs)

    def _spectrum(self, photon_energy):
        """Compute intrinsic synchrotron differential spectrum for energies in ``photon_energy``

        Compute synchrotron for random magnetic field according to approximation
        of Aharonian, Kelner, and Prosekin 2010, PhysRev D 82, 3002
        (`arXiv:1006.1045 <http://arxiv.org/abs/1006.1045>`_).

        Parameters
        ----------
        photon_energy : :class:`~astropy.units.Quantity` instance
            Photon energy array.
        """

        outspecene = _validate_ene(photon_energy)

        from scipy.special import cbrt

        if not hasattr(self,'nB'):
            self.nB = 11

        def Gtilde(x):
            """
            AKP10 Eq. D7

            Factor ~2 performance gain in using cbrt(x)**n vs x**(n/3.)
            """
            gt1 = 1.808 * cbrt(x) / np.sqrt(1 + 3.4 * cbrt(x) ** 2.)
            gt2 = 1 + 2.210 * cbrt(x) ** 2. + 0.347 * cbrt(x) ** 4.
            gt3 = 1 + 1.353 * cbrt(x) ** 2. + 0.217 * cbrt(x) ** 4.
            return gt1 * (gt2 / gt3) * np.exp(-x)

        log.debug('calc_sy: Starting synchrotron computation with AKB2010...')

        def _calc_sy(B):
            # strip units, ensuring correct conversion
            # astropy units do not convert correctly for gyroradius calculation when using
            # cgs (SI is fine, see https://github.com/astropy/astropy/issues/1687)
            CS1_0 = np.sqrt(3) * e.value ** 3 * B.to('G').value
            CS1_1 = (2 * np.pi * m_e.cgs.value * c.cgs.value ** 2 *
                     hbar.cgs.value * outspecene.to('erg').value)
            CS1 = CS1_0/CS1_1
            Ec = 3 * e.value * hbar.cgs.value * B.to('G').value * self._gam ** 2
            Ec /= 2 * (m_e * c).cgs.value
            EgEc = outspecene.to('erg').value / np.vstack(Ec)
            dNdE = CS1 * Gtilde(EgEc)
            # return units
            specB = trapz_loglog(np.vstack(self._nelec) * dNdE, self._gam, axis=0) / u.s / u.erg

            return specB.to('1/(s eV)')

        Bdist = self.B_distribution

        if Bdist.Bmin == Bdist.Bmax:
            total_spec = _calc_sy(Bdist.Bmin)
            self.Brms = Bdist.Bmin
        else:
            self.Brms = Bdist.rms()

            Bs=np.logspace(np.log10(Bdist.Bmin.to('G').value),
                           np.log10(Bdist.Bmax.to('G').value), self.nB) * u.G
            #Bfunc=lambda x: plpdf(x,1.5,Bs.min(),Bs.max()) # alpha=3/2 --> MHD turbulence spectrum: Kraichnan 1965

            sampleBrms = np.sqrt(trapz_loglog(Bs**2 * Bdist.pdf(Bs), Bs) /
                                 trapz_loglog(Bdist.pdf(Bs), Bs))
            log.debug('Relative difference between ideal B RMS and sample B RMS: {0}'.format(
                np.abs(self.Brms - sampleBrms) / self.Brms))
            log.debug('Brms={0}, sampleBrms={1}'.format(self.Brms,sampleBrms))

            spec = []
            for B in Bs:
                spec.append(Bdist.pdf(B) * _calc_sy(B))
            total_spec = trapz_loglog(u.Quantity(spec), Bs.to('G').value, axis=0)

            # not really needed for trapz_loglog, do it just in case
            #samplingcorr = 1/trapz_loglog(Bdist.pdf(Bs), Bs.to('G').value)
            samplingcorr = self.Brms**2 / sampleBrms**2
            log.debug('Sampling correction: {0}'.format(samplingcorr))
            total_spec *= samplingcorr

            flux = trapz_loglog(outspecene*total_spec,outspecene).to('erg/s')
            log.debug('Flux: {0}'.format(flux))

        return total_spec


if __name__ == "__main__":
# test all above
    from pylab import figure

    from naima.models import ExponentialCutoffPowerLaw as ECPL

    pdist = ECPL(4.3e30/u.eV,10*u.TeV,0.835,30*u.TeV)
    bdist = PowerLawB(1.5)
    sync = SynchrotronBdist(pdist,bdist)

    print('We = {0}'.format(sync.We))

    f = figure()
    ax = f.add_subplot(111)

    ene = np.logspace(-3,5,1000)*u.keV

    ax.loglog(ene,sync.sed(ene),label='Brms = {0:.2e}'.format(sync.B_distribution.rms()))

    sync.B_distribution.Bmax /= 10

    ax.loglog(ene,sync.sed(ene),label='Brms = {0:.2e}'.format(sync.B_distribution.rms()))

    Brms = sync.B_distribution.rms()
    sync.B_distribution.Bmax = Brms
    sync.B_distribution.Bmin = Brms

    ax.loglog(ene,sync.sed(ene),label='B = {0:.2e}'.format(Brms))

    ax.legend(loc='lower left')
    ax.set_xlabel('Energy [{0}]'.format(ene.unit))
    ax.set_ylabel('$E^2 dN/dE$')

    f.savefig('bdist_test_plots/magnetic_dist_test.png')


    f = figure()
    ax = f.add_subplot(111)

    pdist = ECPL(4.3e30/u.eV,10*u.TeV,0.835,30*u.TeV)
    bdist = PowerLawB(1.5)
    sync = SynchrotronBdist(pdist,bdist)
    ene = np.logspace(-3,5,1000)*u.keV

    for nB in [5,11,21,51,101,201]:
        sync.nB = nB
        sed = sync.sed(ene)
        ax.loglog(ene,sed,label='nB = {0}'.format(nB))
    ax.legend(loc='lower left')
    ax.set_xlabel('Energy [{0}]'.format(ene.unit))
    ax.set_ylabel('$E^2 dN/dE$')
    ax.set_ylim(bottom=np.max(sed.value)/1e4)

    f.savefig('bdist_test_plots/magnetic_dist_nB.png')