Examples

Fitting a minimal radiative model

#!/usr/bin/env python
import numpy as np
import naima
import astropy.units as u
from astropy.io import ascii
from naima.models import InverseCompton, ExponentialCutoffPowerLaw

## Read data

data = ascii.read('RXJ1713_HESS_2007.dat')


def ElectronIC(pars, data):
    """
    Define particle distribution model, radiative model, and return model flux
    at data energy values
    """

    ECPL = ExponentialCutoffPowerLaw(pars[0] / u.eV, 10. * u.TeV, pars[1],
                                     10**pars[2] * u.TeV)
    IC = InverseCompton(ECPL, seed_photon_fields=['CMB'])

    return IC.flux(data, distance=1.0 * u.kpc)


def lnprior(pars):
    # Limit amplitude to positive domain
    logprob = naima.uniform_prior(pars[0], 0., np.inf)
    return logprob


if __name__ == '__main__':

    ## Set initial parameters and labels
    p0 = np.array((1e30, 3.0, np.log10(30),))
    labels = ['norm', 'index', 'log10(cutoff)']

    ## Run sampler
    sampler, pos = naima.run_sampler(data_table=data,
                                     p0=p0,
                                     labels=labels,
                                     model=ElectronIC,
                                     prior=lnprior,
                                     nwalkers=32,
                                     nburn=100,
                                     nrun=20,
                                     threads=4,
                                     prefit=True,
                                     interactive=False)
    ## Save run results
    out_root = 'RXJ1713_IC_minimal'
    naima.save_run(out_root, sampler)

    ## Save diagnostic plots and results table
    naima.save_diagnostic_plots(out_root, sampler, sed=False)
    naima.save_results_table(out_root, sampler)

Fitting IC emission from an electron distribution

#!/usr/bin/env python
import numpy as np
import naima
import astropy.units as u
from astropy.io import ascii

## Read data

data = ascii.read('RXJ1713_HESS_2007.dat')

## Model definition

from naima.models import InverseCompton, ExponentialCutoffPowerLaw


def ElectronIC(pars, data):

    # Match parameters to ECPL properties, and give them the appropriate units
    amplitude = pars[0] / u.eV
    alpha = pars[1]
    e_cutoff = (10**pars[2]) * u.TeV

    # Initialize instances of the particle distribution and radiative model
    ECPL = ExponentialCutoffPowerLaw(amplitude, 10. * u.TeV, alpha, e_cutoff)
    # Compute IC on CMB and on a FIR component with values from GALPROP for the
    # position of RXJ1713
    IC = InverseCompton(
        ECPL,
        seed_photon_fields=['CMB', ['FIR', 26.5 * u.K, 0.415 * u.eV / u.cm**3]],
        Eemin=100 * u.GeV)

    # compute flux at the energies given in data['energy'], and convert to units
    # of flux data
    model = IC.flux(data, distance=1.0 * u.kpc).to(data['flux'].unit)

    # Save this realization of the particle distribution function
    elec_energy = np.logspace(11, 15, 100) * u.eV
    nelec = ECPL(elec_energy)

    # Compute and save total energy in electrons above 1 TeV
    We = IC.compute_We(Eemin=1 * u.TeV)

    # The first array returned will be compared to the observed spectrum for
    # fitting. All subsequent objects will be stores in the sampler metadata
    # blobs.
    return model, (elec_energy, nelec), We

## Prior definition


def lnprior(pars):
    """
    Return probability of parameter values according to prior knowledge.
    Parameter limits should be done here through uniform prior ditributions
    """

    logprob = naima.uniform_prior(pars[0], 0., np.inf) \
                + naima.uniform_prior(pars[1], -1, 5)

    return logprob


if __name__ == '__main__':

    ## Set initial parameters and labels
    p0 = np.array((1e30, 3.0, np.log10(30),))
    labels = ['norm', 'index', 'log10(cutoff)']

    ## Run sampler
    sampler, pos = naima.run_sampler(data_table=data,
                                     p0=p0,
                                     labels=labels,
                                     model=ElectronIC,
                                     prior=lnprior,
                                     nwalkers=32,
                                     nburn=100,
                                     nrun=20,
                                     threads=4,
                                     prefit=True)

    ## Save run results to HDF5 file (can be read later with naima.read_run)
    naima.save_run('RXJ1713_IC_run.hdf5', sampler)

    ## Diagnostic plots with labels for the metadata blobs
    naima.save_diagnostic_plots(
        'RXJ1713_IC',
        sampler,
        sed=True,
        last_step=False,
        blob_labels=['Spectrum', 'Electron energy distribution',
                     '$W_e (E_e>1\, \mathrm{TeV})$'])
    naima.save_results_table('RXJ1713_IC', sampler)

Fitting Synchrotron and IC emission from an electron distribution

#!/usr/bin/env python
import numpy as np
import astropy.units as u
from astropy.constants import m_e, c
from astropy.io import ascii

import naima

## Read data

# We only consider every fifth X-ray spectral point to speed-up calculations for this example
# DO NOT do this for a final analysis!
soft_xray = ascii.read('RXJ1713_Suzaku-XIS.dat')[::5]
vhe = ascii.read('RXJ1713_HESS_2007.dat')

## Model definition

from naima.models import InverseCompton, Synchrotron, ExponentialCutoffPowerLaw


def ElectronSynIC(pars, data):

    # Match parameters to ECPL properties, and give them the appropriate units
    amplitude = 10**pars[0] / u.eV
    alpha = pars[1]
    e_cutoff = (10**pars[2]) * u.TeV
    B = pars[3] * u.uG

    # Initialize instances of the particle distribution and radiative models
    ECPL = ExponentialCutoffPowerLaw(amplitude, 10. * u.TeV, alpha, e_cutoff)
    # Compute IC on CMB and on a FIR component with values from GALPROP for the
    # position of RXJ1713
    IC = InverseCompton(
        ECPL,
        seed_photon_fields=['CMB', ['FIR', 26.5 * u.K, 0.415 * u.eV / u.cm**3]],
        Eemin=100 * u.GeV)
    SYN = Synchrotron(ECPL, B=B)

    # compute flux at the energies given in data['energy']
    model = (IC.flux(data,
                     distance=1.0 * u.kpc) + SYN.flux(data,
                                                      distance=1.0 * u.kpc))

    # The first array returned will be compared to the observed spectrum for
    # fitting. All subsequent objects will be stored in the sampler metadata
    # blobs.
    return model, IC.compute_We(Eemin=1 * u.TeV)

## Prior definition


def lnprior(pars):
    """
    Return probability of parameter values according to prior knowledge.
    Parameter limits should be done here through uniform prior ditributions
    """
    # Limit norm and B to be positive
    logprob = naima.uniform_prior(pars[0], 0., np.inf) \
                + naima.uniform_prior(pars[1], -1, 5) \
                + naima.uniform_prior(pars[3], 0, np.inf)

    return logprob


if __name__ == '__main__':

    ## Set initial parameters and labels
    # Estimate initial magnetic field and get value in uG
    B0 = 2 * naima.estimate_B(soft_xray, vhe).to('uG').value

    p0 = np.array((33, 2.5, np.log10(48.0), B0))
    labels = ['log10(norm)', 'index', 'log10(cutoff)', 'B']

    ## Run sampler
    sampler, pos = naima.run_sampler(data_table=[soft_xray, vhe],
                                     p0=p0,
                                     labels=labels,
                                     model=ElectronSynIC,
                                     prior=lnprior,
                                     nwalkers=32,
                                     nburn=100,
                                     nrun=20,
                                     threads=4,
                                     prefit=True,
                                     interactive=False)

    ## Save run results to HDF5 file (can be read later with naima.read_run)
    naima.save_run('RXJ1713_SynIC', sampler)

    ## Diagnostic plots
    naima.save_diagnostic_plots('RXJ1713_SynIC',
                                sampler,
                                sed=True,
                                blob_labels=['Spectrum', '$W_e$($E_e>1$ TeV)'])
    naima.save_results_table('RXJ1713_SynIC', sampler)

Additional model expressions

import numpy as np
import astropy.units as u
import naima

################################################################################
#
# This file shows a few example model functions (with associated priors, labels
# and p0 vector), that can be used as input for naima.run_sampler
#
################################################################################

#
# RADIATIVE MODELS
#
# Pion decay
# ==========

PionDecay_ECPL_p0 = np.array((46, 2.34, np.log10(80.)))
PionDecay_ECPL_labels = ['log10(norm)', 'index', 'log10(cutoff)']

# Prepare an energy array for saving the particle distribution
proton_energy = np.logspace(-3, 2, 50) * u.TeV


def PionDecay_ECPL(pars, data):
    amplitude = 10**pars[0] / u.TeV
    alpha = pars[1]
    e_cutoff = 10**pars[2] * u.TeV

    ECPL = naima.models.ExponentialCutoffPowerLaw(amplitude, 30 * u.TeV, alpha,
                                                  e_cutoff)
    PP = naima.models.PionDecay(ECPL, nh=1.0 * u.cm** -3)

    model = PP.flux(data, distance=1.0 * u.kpc)
    # Save a realization of the particle distribution to the metadata blob
    proton_dist = PP.particle_distribution(proton_energy)
    # Compute the total energy in protons above 1 TeV for this realization
    Wp = PP.compute_Wp(Epmin=1 * u.TeV)

    # Return the model, proton distribution and energy in protons to be stored
    # in metadata blobs
    return model, (proton_energy, proton_dist), Wp


def PionDecay_ECPL_lnprior(pars):
    logprob = naima.uniform_prior(pars[1], -1, 5)
    return logprob

# Inverse Compton with the energy in electrons as the normalization parameter
# ===========================================================================

IC_We_p0 = np.array((40, 3.0, np.log10(30)))
IC_We_labels = ['log10(We)', 'index', 'log10(cutoff)']


def IC_We(pars, data):
    # Example of a model that is normalized though the total energy in electrons

    # Match parameters to ECPL properties, and give them the appropriate units
    We = 10**pars[0] * u.erg
    alpha = pars[1]
    e_cutoff = 10**pars[2] * u.TeV

    # Initialize instances of the particle distribution and radiative model
    # set a bogus normalization that will be changed in third line
    ECPL = ExponentialCutoffPowerLaw(1 / u.eV, 10. * u.TeV, alpha, e_cutoff)
    IC = InverseCompton(ECPL, seed_photon_fields=['CMB'])
    IC.set_We(We, Eemin=1 * u.TeV)

    # compute flux at the energies given in data['energy']
    model = IC.flux(data, distance=1.0 * u.kpc)

    # Save this realization of the particle distribution function
    elec_energy = np.logspace(11, 15, 100) * u.eV
    nelec = ECPL(elec_energy)

    return model, (elec_energy, nelec)


def IC_We_lnprior(pars):
    logprob = naima.uniform_prior(pars[1], -1, 5)
    return logprob

#
# FUNCTIONAL MODELS
#
# Exponential cutoff powerlaw
# ===========================

ECPL_p0 = np.array((1e-12, 2.4, np.log10(15.0)))
ECPL_labels = ['norm', 'index', 'log10(cutoff)']


def ECPL(pars, data):
    # Get the units of the flux data and match them in the model amplitude
    amplitude = pars[0] * data['flux'].unit
    alpha = pars[1]
    e_cutoff = (10**pars[2]) * u.TeV
    ECPL = naima.models.ExponentialCutoffPowerLaw(amplitude, 1 * u.TeV, alpha,
                                                  e_cutoff)

    return ECPL(data)


def ECPL_lnprior(pars):
    logprob = naima.uniform_prior(pars[0], 0., np.inf) \
            + naima.uniform_prior(pars[1], -1, 5)
    return logprob

# Log-Parabola or Curved Powerlaw
# ===============================

LP_p0 = np.array((1.5e-12, 2.7, 0.12,))
LP_labels = ['norm', 'alpha', 'beta']


def LP(pars, data):
    amplitude = pars[0] * data['flux'].unit
    alpha = pars[1]
    beta = pars[2]
    LP = naima.models.LogParabola(amplitude, 1 * u.TeV, alpha, beta)
    return LP(data)


def LP_lnprior(pars):
    logprob = naima.uniform_prior(pars[0], 0., np.inf) \
                + naima.uniform_prior(pars[1], -1, 5)
    return logprob

Crab Nebula SSC model

import numpy as np
from astropy.io import ascii
from astropy.constants import c
import astropy.units as u
import matplotlib.pyplot as plt
import naima
from naima.models import *

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('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
ax.loglog(energy, IC.sed(energy, 2 * u.kpc) + SYN.sed(energy, 2 * u.kpc),
          lw=3, c='k', label='Total')
for i, seed, ls in zip(
        range(4), ['CMB', 'FIR', 'NIR', 'SSC'], ['--', '-.', ':', '-']):
    ax.loglog(energy, IC.sed(energy, 2 * u.kpc, seed=seed),
              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()
figure.savefig('CrabNebula_SynSSC.png')