Examples

Fitting a minimal radiative model

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

import naima
from naima.models import ExponentialCutoffPowerLaw, InverseCompton

## 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.0 * 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.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 + ".hdf5", 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 astropy.units as u
import numpy as np
from astropy.io import ascii

import naima
from naima.models import ExponentialCutoffPowerLaw, InverseCompton

## Read data

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

## Model definition


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.0 * 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.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 astropy.units as u
import numpy as np
from astropy.io import ascii

import naima
from naima.models import ExponentialCutoffPowerLaw, InverseCompton, Synchrotron

## 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


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.0 * 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.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.hdf5", 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 astropy.units as u
import numpy as np

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.0)))
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 = naima.models.ExponentialCutoffPowerLaw(
        1 / u.eV, 10.0 * u.TeV, alpha, e_cutoff
    )
    IC = naima.models.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.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.0, np.inf) + naima.uniform_prior(
        pars[1], -1, 5
    )
    return logprob

Crab Nebula SSC model

import astropy.units as u
import numpy as np
from astropy.constants import c
from astropy.io import ascii

import naima
from naima.models import (
    ExponentialCutoffBrokenPowerLaw,
    InverseCompton,
    Synchrotron,
)

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.0,
)

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")