In [1]:

    

import numpy as np
from astropy.io import ascii
from astropy.constants import c
import astropy.units as u
import naima
from naima.models import (ExponentialCutoffBrokenPowerLaw, Synchrotron,
                          InverseCompton)

%matplotlib inline

ECBPL = ExponentialCutoffBrokenPowerLaw(amplitude=3.699e36 / u.eV,
                                        e_0=1 * u.TeV,
                                        e_break=0.265 * u.TeV,
                                        alpha_1=1.5,
                                        alpha_2=3.233,
                                        e_cutoff=1863 * u.TeV,
                                        beta=2.)

eopts = {'Eemax': 50 * u.PeV, 'Eemin': 0.1 * u.GeV}

SYN = Synchrotron(ECBPL, B=125 * u.uG, Eemax=50 * u.PeV, Eemin=0.1 * u.GeV)

# Compute photon density spectrum from synchrotron emission assuming R=2.1 pc
Rpwn = 2.1 * u.pc
Esy = np.logspace(-7, 9, 100) * u.eV
Lsy = SYN.flux(Esy, distance=0 * u.cm)  # use distance 0 to get luminosity
phn_sy = Lsy / (4 * np.pi * Rpwn**2 * c) * 2.24

IC = InverseCompton(ECBPL,
                    seed_photon_fields=['CMB',
                                        ['FIR', 70 * u.K, 0.5 * u.eV / u.cm**3],
                                        ['NIR', 5000 * u.K, 1 * u.eV / u.cm**3],
                                        ['SSC', Esy, phn_sy]],
                    Eemax=50 * u.PeV, Eemin=0.1 * u.GeV)

# Use plot_data from naima to plot the observed spectra
data = ascii.read('/home/giacomov/software/naima/examples/CrabNebula_spectrum.ecsv')
figure = naima.plot_data(data, e_unit=u.eV)
ax = figure.axes[0]

# Plot the computed model emission
energy = np.logspace(-7, 15, 100) * u.eV

y = IC.sed(energy, 2 * u.kpc) + SYN.sed(energy, 2 * u.kpc)

ax.loglog(energy.value, y.value,
          lw=3, c='k', label='Total')

ax.loglog(energy.value, SYN.sed(energy, 2 * u.kpc).value, lw=2,c='blue',label='Synch')

for i, seed, ls in zip(
        range(4), ['CMB', 'FIR', 'NIR', 'SSC'], ['--', '-.', ':', '-']):
    ax.loglog(energy.value, IC.sed(energy, 2 * u.kpc, seed=seed).value,
              lw=2, c=naima.plot.color_cycle[i + 1], label=seed, ls=ls)


ax.set_ylim(1e-12, 1e-7)
ax.legend(loc='upper right', frameon=False)
figure.tight_layout()


    

In [14]:

    

import matplotlib.pyplot as plt

energies = np.logspace(0,17,100)

plt.loglog(energies,ECBPL(energies * u.eV).value)

plt.ylim([1e-12,1e56])


    

    
Out[14]:





(1e-12, 1e+56)






    

In [16]:

    

ECBPL(energies * u.eV)


    

    
Out[16]:





$[3.699 \times 10^{54},~2.0441303 \times 10^{54},~1.1296212 \times 10^{54},~6.2424785 \times 10^{53},~3.4496998 \times 10^{53},~1.9063628 \times 10^{53},~1.0534885 \times 10^{53},~5.8217564 \times 10^{52},~3.2172016 \times 10^{52},~1.7778803 \times 10^{52},~9.8248687 \times 10^{51},~5.4293895 \times 10^{51},~3.0003729 \times 10^{51},~1.6580571 \times 10^{51},~9.1627054 \times 10^{50},~5.0634668 \times 10^{50},~2.7981578 \times 10^{50},~1.5463096 \times 10^{50},~8.5451698 \times 10^{49},~4.7222061 \times 10^{49},~2.6095714 \times 10^{49},~1.4420935 \times 10^{49},~7.9692539 \times 10^{48},~4.4039452 \times 10^{48},~2.433695 \times 10^{48},~1.3449012 \times 10^{48},~7.4321529 \times 10^{47},~4.107134 \times 10^{47},~2.2696721 \times 10^{47},~1.2542594 \times 10^{47},~6.9312506 \times 10^{46},~3.830327 \times 10^{46},~2.1167038 \times 10^{46},~1.1697265 \times 10^{46},~6.4641076 \times 10^{45},~3.5721757 \times 10^{45},~1.974045 \times 10^{45},~1.0908909 \times 10^{45},~6.0284483 \times 10^{44},~3.331423 \times 10^{44},~1.841001 \times 10^{44},~1.0173685 \times 10^{44},~5.6221511 \times 10^{43},~3.1068963 \times 10^{43},~1.7169237 \times 10^{43},~9.4880123 \times 10^{42},~5.2432369 \times 10^{42},~2.8975019 \times 10^{42},~1.6012088 \times 10^{42},~8.8485521 \times 10^{41},~4.8898603 \times 10^{41},~2.7022199 \times 10^{41},~1.4932927 \times 10^{41},~8.2521893 \times 10^{40},~4.5603001 \times 10^{40},~2.5200994 \times 10^{40},~1.3926498 \times 10^{40},~7.6960193 \times 10^{39},~4.2529511 \times 10^{39},~2.3502531 \times 10^{39},~1.2987898 \times 10^{39},~7.1773334 \times 10^{38},~3.9663165 \times 10^{38},~2.1918539 \times 10^{38},~1.2112557 \times 10^{38},~6.6936051 \times 10^{37},~3.699 \times 10^{37},~1.4748272 \times 10^{37},~4.1075183 \times 10^{36},~1.1439785 \times 10^{36},~3.1860762 \times 10^{35},~8.8734883 \times 10^{34},~2.4713394 \times 10^{34},~6.8828765 \times 10^{33},~1.9169315 \times 10^{33},~5.3387681 \times 10^{32},~1.4868625 \times 10^{32},~4.1408564 \times 10^{31},~1.1531521 \times 10^{31},~3.2109419 \times 10^{30},~8.9385443 \times 10^{29},~2.4868824 \times 10^{29},~6.9103731 \times 10^{28},~1.9149263 \times 10^{28},~5.2743149 \times 10^{27},~1.4333954 \times 10^{27},~3.7822067 \times 10^{26},~9.3508952 \times 10^{25},~2.002715 \times 10^{25},~3.1254528 \times 10^{24},~2.42696 \times 10^{23},~4.0433362 \times 10^{21},~2.261247 \times 10^{18},~7.1029884 \times 10^{11},~0.015160849,~5.0587853 \times 10^{-32},~2.3613558 \times 10^{-96},~1.5282419 \times 10^{-237},~0,~0] \; \mathrm{\frac{1}{eV}}$

In [ ]: