Source model interface

icecube_tools has a simple source modelling interface built in that we demonstrate here.

[1]:
import numpy as np
from matplotlib import pyplot as plt

from icecube_tools.source.flux_model import PowerLawFlux, BrokenPowerLawFlux
from icecube_tools.source.source_model import PointSource, DiffuseSource

Spectral shape

We start by defining a spectral shape, such as a power law or broken power law. Let’s start with the definition of a simple power law flux.

[2]:
# Parameters of power law flux
flux_norm = 1e-18  # Flux normalisation in units of GeV^-1 cm^-2 s^-1 (sr^-1)
norm_energy = 1e5  # Energy of normalisation in units of GeV
spectral_index = 2.0  # Assumed negative slope
min_energy = 1e4  # GeV
max_energy = 1e8  # GeV

# Instantiate
power_law = PowerLawFlux(flux_norm, norm_energy, spectral_index, min_energy, max_energy)
[3]:
energies = np.geomspace(min_energy, max_energy)

fig, ax = plt.subplots()
ax.plot(energies, [power_law.spectrum(e) for e in energies])
ax.axhline(flux_norm, color="k", linestyle=":")
ax.axvline(norm_energy, color="k", linestyle=":")
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_xlabel("E [GeV]")
ax.set_ylabel("F [GeV^-1 cm^-2 s^-1 (sr^-1)]")
[3]:
Text(0, 0.5, 'F [GeV^-1 cm^-2 s^-1 (sr^-1)]')
../_images/notebooks_source_model_5_1.png

We can also use the PowerLawFlux class to perform some simple calculations, such as integration of the flux.

[4]:
total_flux = power_law.integrated_spectrum(min_energy, max_energy)  # cm^-2 s^-1 (sr^-1)
total_flux
[4]:
9.999000000000002e-13
[5]:
total_energy_flux = power_law.total_flux_density()  # GeV cm^-2 s^-1 (sr^-1)
total_energy_flux
[5]:
9.210340371976183e-08

Sampling from the power law shape is also possible:

[6]:
samples = power_law.sample(1000)

fig, ax = plt.subplots()
ax.hist(samples, bins=energies, density=True, label="Samples")
ax.plot(energies, [power_law.spectrum(e) / total_flux for e in energies], label="Model")
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_xlabel("E [GeV]")
ax.legend()
[6]:
<matplotlib.legend.Legend at 0x7fac9e7f26d0>
../_images/notebooks_source_model_10_1.png

The BrokenPowerLaw class is also available and behaves in a very similar way.

Diffuse and point sources

Once the spectral shape is defined, we can specify either a DiffuseSource or a PointSource. It is assumed that diffuse sources are isotropic and the flux model describes the per-steradian flux over the entire \(4\pi\) sky. We also specify a redshift of the source such that adiabatic neutrino energy losses can be accounted for. Naturally, PointSource objects also have a direction specified in (ra, dec) coordinates.

[7]:
diffuse_source = DiffuseSource(power_law, z=0.0)

ra = np.deg2rad(50)
dec = np.deg2rad(-10)
point_source = PointSource(power_law, z=0.5, coord=(ra, dec))

The original flux model can now be accessed from within the source along with its other properties:

[8]:
diffuse_source.flux_model
[8]:
<icecube_tools.source.flux_model.PowerLawFlux at 0x7facd82dcc10>

Sources and lists of sources can be used as input to simulations, as demonstrated in the simulation notebook.

[ ]: