# 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)]')
```

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

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.

```
[ ]:
```

```
```