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