# Extended source tutorial

## Introduction

Extended sources with constant (not energy-dependent) morphology can be instanced by providing a spectral shape (1D function) and a spatial shape (2D function):

```
[1]:
```

```
from astromodels import *
disk_source = ExtendedSource('simple_source', spectral_shape=Powerlaw(), spatial_shape = Disk_on_sphere() )
disk_source.display()
```

22:58:53 WARNING The naima package is not available. Models that depend on it will not be functions.py:48 available

WARNING The GSL library or the pygsl wrapper cannot be loaded. Models that depend on it functions.py:69 will not be available.

```
/Users/runner/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages/numba/core/decorators.py:262: NumbaDeprecationWarning: numba.generated_jit is deprecated. Please see the documentation at: https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-generated-jit for more information and advice on a suitable replacement.
warnings.warn(msg, NumbaDeprecationWarning)
```

```
/Users/runner/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages/numba/core/decorators.py:262: NumbaDeprecationWarning: numba.generated_jit is deprecated. Please see the documentation at: https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-generated-jit for more information and advice on a suitable replacement.
warnings.warn(msg, NumbaDeprecationWarning)
```

- simple_source (extended source):
- shape:
- lon0:
- value: 0.0
- desc: Longitude of the center of the source
- min_value: 0.0
- max_value: 360.0
- unit: deg
- is_normalization: False

- lat0:
- value: 0.0
- desc: Latitude of the center of the source
- min_value: -90.0
- max_value: 90.0
- unit: deg
- is_normalization: False

- radius:
- value: 15.0
- desc: Radius of the disk
- min_value: 0.0
- max_value: 20.0
- unit: deg
- is_normalization: False

- lon0:
- spectrum:
- main:
- Powerlaw:
- K:
- value: 1.0
- desc: Normalization (differential flux at the pivot value)
- min_value: 1e-30
- max_value: 1000.0
- unit: keV-1 s-1 cm-2
- is_normalization: True

- piv:
- value: 1.0
- desc: Pivot value
- min_value: None
- max_value: None
- unit: keV
- is_normalization: False

- index:
- value: -2.01
- desc: Photon index
- min_value: -10.0
- max_value: 10.0
- unit:
- is_normalization: False

- K:

- Powerlaw:

- main:

- shape:

Instead of the `position`

property, extended sources have a `spatial_shape`

property, which contains the function that was specified as the `shape`

argument earlier. Information can be printed either with the regular `print`

function or using `.display()`

:

```
[2]:
```

```
print("Name", disk_source.spatial_shape.name)
print("Class", type(disk_source.spatial_shape) )
disk_source.spatial_shape.display()
```

```
Name Disk_on_sphere
Class <class 'astromodels.functions.functions_2D.Disk_on_sphere'>
```

- description: A bidimensional disk/tophat function on a sphere (in spherical coordinates)
- formula: $$ f(\vec{x}) = \left(\frac{180}{\pi}\right)^2 \frac{1}{\pi~({\rm radius})^2} ~\left\{\begin{matrix} 1 & {\rm if}& {\rm | \vec{x} - \vec{x}_0| \le {\rm radius}} \\ 0 & {\rm if}& {\rm | \vec{x} - \vec{x}_0| > {\rm radius}} \end{matrix}\right. $$
- parameters:
- lon0:
- value: 0.0
- desc: Longitude of the center of the source
- min_value: 0.0
- max_value: 360.0
- unit: deg
- is_normalization: False
- delta: 0.1
- free: True

- lat0:
- value: 0.0
- desc: Latitude of the center of the source
- min_value: -90.0
- max_value: 90.0
- unit: deg
- is_normalization: False
- delta: 0.1
- free: True

- radius:
- value: 15.0
- desc: Radius of the disk
- min_value: 0.0
- max_value: 20.0
- unit: deg
- is_normalization: False
- delta: 1.5
- free: True

- lon0:

Spatial shape parameters can be listed, accessed, changed as any function parameters, for example:

```
[3]:
```

```
import astropy.units as u
print( disk_source.spatial_shape.parameters.keys() )
print( disk_source.parameters.keys() )
disk_source.spatial_shape.radius.display()
disk_source.spatial_shape.radius = 0.3 * u.deg
disk_source.spatial_shape.radius.bounds = (0.1, 1.0)*u.deg
disk_source.spatial_shape.radius.display()
```

```
odict_keys(['lon0', 'lat0', 'radius'])
odict_keys(['simple_source.spectrum.main.Powerlaw.K', 'simple_source.spectrum.main.Powerlaw.piv', 'simple_source.spectrum.main.Powerlaw.index', 'simple_source.Disk_on_sphere.lon0', 'simple_source.Disk_on_sphere.lat0', 'simple_source.Disk_on_sphere.radius'])
```

## Available shapes

Here’s a list of available 2D functions:

```
[4]:
```

```
print([c.__name__ for c in functions.function.Function2D.__subclasses__()])
```

```
['Latitude_galactic_diffuse', 'Gaussian_on_sphere', 'Asymm_Gaussian_on_sphere', 'Disk_on_sphere', 'Ellipse_on_sphere', 'SpatialTemplate_2D', 'Power_law_on_sphere']
```

Check the API or use the `display()`

function to find out more about the shapes and their free parameters.

### “Moveable” Functions

`Gaussian_on_sphere`

, `Asymm_Gaussian_on_sphere`

, `Disk_on_sphere`

, `Ellipse_on_sphere`

, and `Power_law_on_sphere`

each have parameters called `lat0`

and `lon0`

, describing the longitude (R.A.) and latitude (Declination) of the function’s “center” (center of mass/center of symmetry). These parameters can be set free to allow the source position to vary during the fit, or fixed to a known position. **If fitting the source position, it is strongly recommended to restrict these
parameters so that source overlaps the region of interest (ROI) at all times.** The source moving far enough from the ROI during the fit can lead to issues such as the minimizer getting “stuck” due to the likelihood surface being flat.

Each of these functions also has other parameter(s) related to the source’s extend and (if applicable) its orientation.

All of these functions are normalized so that their integral is 1. The 1D `spectrum`

function of an extended source can thus be interpreted as the spectrum that would be observed by a non-imaging instrument with a field of view much larger than the source in question.

### Fixed-position Functions

`Latitude_galactic_diffuse`

and `SpatialTemplate_2D`

do not have parameters named `lat0`

and `lon0`

. Their positions are always fixed throughout the fit.

`Latitude_galactic_diffuse`

is a function that is constant in Galactic longitude between `l_min`

and `l_max`

and follows a Gaussian shape in Galactic latitude, with a “width” parameter (`sigma_b`

).

```
[5]:
```

```
diff = Latitude_galactic_diffuse()
diff.display()
```

- description: A Gaussian distribution in Galactic latitude around the Galactic plane
- formula: $ K \exp{\left( \frac{-b^2}{2 \sigma_b^2} \right)} $
- parameters:
- K:
- value: 1.0
- desc: normalization
- min_value: None
- max_value: None
- unit:
- is_normalization: False
- delta: 0.1
- free: True

- sigma_b:
- value: 1.0
- desc: Sigma for
- min_value: None
- max_value: None
- unit:
- is_normalization: False
- delta: 0.1
- free: True

- l_min:
- value: 10.0
- desc: min Longtitude
- min_value: None
- max_value: None
- unit:
- is_normalization: False
- delta: 1.0
- free: True

- l_max:
- value: 30.0
- desc: max Longtitude
- min_value: None
- max_value: None
- unit:
- is_normalization: False
- delta: 3.0
- free: True

- K:

`SpatialTemplate_2D`

is designed to read in a user-provided fits file with an image (in WCS coordinates) using the `load_file()`

function. The function value will be 0 outside the WCS and equal to the value of the pixel containing the given coordinates inside the WCS. The provided template should be in units of 1/deg2 and normalized so that its integral is 1. If that is not the case, the normalization parameter `K`

may be set accordingly so that the overall function is normalized as
expected.

```
[6]:
```

temp = SpatialTemplate_2D(fits_file='spitzer_example_image.fits')
temp.display()

22:58:54 WARNING 2D template read from spitzer_example_image.fits is normalized to nan functions_2D.py:747 (expected: 1)

- description: User input Spatial Template. Expected to be normalized to 1/sr
- formula: $ hi $
- parameters:
- K:
- value: 1.0
- desc: normalization
- min_value: None
- max_value: None
- unit:
- is_normalization: False
- delta: 0.1
- free: False

- hash:
- value: 7.428580553431191e+66
- desc: hash of model map [needed for memoization]
- min_value: None
- max_value: None
- unit:
- is_normalization: False
- delta: 0.1
- free: False

- ihdu:
- value: 0.0
- desc: header unit index of fits file
- min_value: 0.0
- max_value: None
- unit:
- is_normalization: False
- delta: 0.1
- free: False

- fits_file:
- value: spitzer_example_image.fits
- desc: fits file to load
- allowed values: None
- defer: True
- function: _load_file

- frame:
- value: icrs
- desc: coordinate frame
- allowed values: ['icrs', 'galactic', 'fk5', 'fk4', 'fk4_no_e']
- defer: False
- function: None

- K:

## “Calling” extended sources

Extended sources can be called as functions. They take 3 arguments: A list or `np.array`

of right ascensions, a list of declinations, and a list of energies. The first two lists must have the same dimensions. The result contains the value(s) of the double differential flux dN/dE/dOmega at all combinations of coordinates and energies. See the example below. There, we picked two positions (one at the center of the disk source, one outside of it) and three energy values. **If no units are
provided, energies are assumed to be in keV, coordinates in degrees, and fluxes are returned in 1/(cm2 keV s deg2).**

```
[7]:
```

```
import numpy as np
ra_center, dec_center = disk_source.spatial_shape.lon0.value, disk_source.spatial_shape.lat0.value
ra = [ra_center, ra_center+2]*u.deg
dec = [dec_center, dec_center+2]*u.deg
E = [1, 10, 100]*u.keV
print( "Whole source:" )
print( disk_source(ra, dec, E) )
print ("")
print ("Spatial only:" )
spatial = disk_source.spatial_shape(ra, dec)
print( spatial )
print ("")
print ("Spectral only:")
spectral = disk_source.spatial_shape.get_total_spatial_integral( E ) * disk_source.spectrum.main.shape( E )
print ( spectral )
print ("")
print ( "Spatial transposed times spectral")
print( np.atleast_2d(spatial).T * np.atleast_2d(spectral) )
```

```
Whole source:
```

```
[[1.16105524e+04 1.13462640e+02 1.10879915e+00]
[0.00000000e+00 0.00000000e+00 0.00000000e+00]] 1 / (keV s cm2 deg2)
Spatial only:
[11610.55239595 0. ] 1 / deg2
Spectral only:
[1.00000000e+00 9.77237221e-03 9.54992586e-05] 1 / (keV s cm2)
Spatial transposed times spectral
[[1.16105524e+04 1.13462640e+02 1.10879915e+00]
[0.00000000e+00 0.00000000e+00 0.00000000e+00]] 1 / (keV s cm2 deg2)
```

## Caveats and gotchas

### Large extended sources

The Gaussian shapes provided here (`Gaussian_on_sphere`

, `Asymm_Gaussian_on_sphere`

) are only valid as long as their width (`sigma`

/`a`

parameter) is much smaller than 360˚. For very large extended sources, we would need to implement the Kent distribution, the equivalent of a Gaussian on the sphere.

### Small extended sources

Most threeML plugins supporting imaging instruments perform a convolution of the `astromodels`

source(s) with the instrument’s point spread function (PSF). These convolutions are typically performed on a grid, and features in the source morphology smaller than the grid size are lost. One should take care that the extent of the source(s) in the model is constrained to be no smaller than the grid spacing used by the plugin.

### Source support

Each 2D function class has a member function called `get_boundaries`

, which returns the (minimal) range in RA and Dec over which the convolution with the PSF has to be computed (see example below). The exact calculation depends on the function but the range typically depends on the centroid and the **maximum value of the width/radius parameter**. That means that even if the radius or gaussian width parameter is fixed during the fit, it is helpful to set a maximum that is close to the fixed
value. This reduces computation time during the likelihood fit.

```
[8]:
```

```
disk_source.spatial_shape.lon0 = 80*u.deg
disk_source.spatial_shape.lat0 = 20*u.deg
disk_source.spatial_shape.radius.bounds = (0.1, 2)*u.deg
bounds = (disk_source.spatial_shape.get_boundaries() )
print ("RA range:", bounds[0])
print ("Dec range:", bounds[1])
```

```
RA range: (75.62145488597581, 84.37854511402419)
Dec range: (16.0, 24.0)
```

## Energy-Dependent Morphology

`astromodels`

functions support energy-dependent morphology in two different ways.

Extended sources can be instanced with a 3D function spatial shape and a separate energy spectrum. In that case, the 3D function is expected to be a function of RA, Dec, Energy (in that order) and have units of 1/deg2. The double-differential flux dN/dE/dOmega is again given by the product of the spatial and spectral parts. See for example

``Continuous_injection_diffusion`

<https:link>`__.Extended sources can be instanced with just a 3D function and no energy spectrum. In that case, the 3D function is interpreted as the double-differential flux dN/dE/dOmega.

```
[ ]:
```

```
```