Point sources

In astromodels, a point source is described by its position in the sky and its spectral features.

Creating a point source

A simple source with a power law spectrum can be created like this, using J2000 R.A. and Dec (ICRS), which is the default coordinate system:

[1]:
from astromodels import *
simple_source_icrs = PointSource('simple_source', ra=123.2, dec=-13.2, spectral_shape=Powerlaw())
22:59:27 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)

We can also use Galactic coordinates:

[2]:
simple_source_gal = PointSource('simple_source', l=234.320573, b=11.365142, spectral_shape=Powerlaw())

As spectral shape we can use any function or any composite function (see “Creating and modifying functions”)

Getting info about a point source

Info about a point source can be obtained with the .display() method (which will use the richest representation available), or by printing it which will display a text-only representation:

[3]:
simple_source_icrs.display()
  • simple_source (point source):
    • position:
      • ra:
        • value: 123.2
        • desc: Right Ascension
        • min_value: 0.0
        • max_value: 360.0
        • unit: deg
        • is_normalization: False
      • dec:
        • value: -13.2
        • desc: Declination
        • min_value: -90.0
        • max_value: 90.0
        • unit: deg
        • is_normalization: False
      • equinox: J2000
    • 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
[4]:
print(simple_source_icrs)
  * simple_source (point source):
    * position:
      * ra:
        * value: 123.2
        * desc: Right Ascension
        * min_value: 0.0
        * max_value: 360.0
        * unit: deg
        * is_normalization: false
      * dec:
        * value: -13.2
        * desc: Declination
        * min_value: -90.0
        * max_value: 90.0
        * unit: deg
        * is_normalization: false
      * equinox: J2000
    * spectrum:
      * main:
        * Powerlaw:
          * K:
            * value: 1.0
            * desc: Normalization (differential flux at the pivot value)
            * min_value: 1.0e-30
            * max_value: 1000.0
            * unit: keV-1 s-1 cm-2
            * is_normalization: true
          * piv:
            * value: 1.0
            * desc: Pivot value
            * min_value: null
            * max_value: null
            * unit: keV
            * is_normalization: false
          * index:
            * value: -2.01
            * desc: Photon index
            * min_value: -10.0
            * max_value: 10.0
            * unit: ''
            * is_normalization: false
        * polarization: {}

As you can see we have created a point source with one component automatically named “main”, with a power law spectrum, at the specified position.

Converting between coordinate systems

By default the coordinates of the point source are displayed in the same system used during creation. However, you can always obtain R.A, Dec or L,B like this:

[5]:
simple_source_icrs.position.display()
l = simple_source_icrs.position.get_l()
b = simple_source_icrs.position.get_b()
print(l,b)

simple_source_gal.position.display()
ra = simple_source_gal.position.get_ra()
dec = simple_source_gal.position.get_dec()
print(ra,dec)
Sky direction (R.A., Dec.) = (123.20000, -13.20000) (J2000)
234.32058209778668 11.36514410904394
Sky direction (l, b) = (234.32057, 11.36514) (J2000)
123.19999334752814 -13.199993514950481

For more control on the output and many more options, such as transform to local frames or other equinoxes, you can obtain an instance of astropy.coordinates.SkyCoord by using the sky_coord property of the position object:

[6]:
sky_coord_instance = simple_source_icrs.position.sky_coord
ra = sky_coord_instance.transform_to('icrs').ra
dec = sky_coord_instance.transform_to('icrs').dec
print(ra.deg)
123.2

Gotcha while accessing coordinates

Please note that using get_ra() and .ra (or the equivalent methods for the other coordinates) is not the same. While get_ra() will always return a single float value corresponding to the R.A. of the source, the .ra property will exist only if the source has been created using R.A, Dec as input coordinates and will return a Parameter instance:

[7]:
parameter_ra = simple_source_icrs.position.ra
parameter_dec = simple_source_icrs.position.dec

print( type(parameter_ra) )
parameter_ra.display()
parameter_dec.display()
<class 'astromodels.core.parameter.Parameter'>
Parameter ra = 123.2 [deg] (min_value = 0.0, max_value = 360.0, delta = 12.32, free = False)
Parameter dec = -13.2 [deg] (min_value = -90.0, max_value = 90.0, delta = 1.32, free = False)

The following would instead throw AttributeError, since simple_source_icrs was instanced using R.A. and Dec. and hence does not have the l, b parameters:

[8]:
try:
    print( simple_source_icrs.position.l)
except Exception as e:
    print(e)
Accessing an element l of the node that does not exist

In all cases, independently on how the source was instanced, you can obtain the values of coordinates in degrees as floating point numbers using get_ra(), get_dec(), get_l(), get_b(). However, you can only directly assign coordinates in the same system that the source direction was originally created, e.g.:

[9]:
simple_source_icrs.position.display()
simple_source_icrs.position.ra =   simple_source_icrs.position.ra.value + 1.0
simple_source_icrs.position.dec = simple_source_icrs.position.dec.value - 1.0
simple_source_icrs.position.display()
Sky direction (R.A., Dec.) = (123.20000, -13.20000) (J2000)
Sky direction (R.A., Dec.) = (124.20000, -14.20000) (J2000)

Fitting the source position

Source coordinates, like any parameters, can be set to be free or fixed during the fit. By default, coordinates are set to be fixed. If you would like to fit them as free parameters during the likelihood fit, they can be freed as any other parameter. Note that param.free = True and param.fix = False are equivalent.

[10]:
print("Free parameters (before freeing position):", simple_source_icrs.free_parameters.keys())
simple_source_icrs.position.ra.free = True
simple_source_icrs.position.dec.fix = False
print("Free parameters (after freeing position):", simple_source_icrs.free_parameters.keys())

Free parameters (before freeing position): odict_keys(['simple_source.spectrum.main.Powerlaw.K', 'simple_source.spectrum.main.Powerlaw.index'])
Free parameters (after freeing position): odict_keys(['simple_source.spectrum.main.Powerlaw.K', 'simple_source.spectrum.main.Powerlaw.index', 'simple_source.position.ra', 'simple_source.position.dec'])

For a source created in Galactic coordinates, instead use the following:

[11]:
print("Free parameters (before freeing position):", simple_source_gal.free_parameters.keys())
simple_source_gal.position.l.free = True
simple_source_gal.position.b.fix = False
print("Free parameters (after freeing position):", simple_source_gal.free_parameters.keys())

Free parameters (before freeing position): odict_keys(['simple_source.spectrum.main.Powerlaw.K', 'simple_source.spectrum.main.Powerlaw.index'])
Free parameters (after freeing position): odict_keys(['simple_source.spectrum.main.Powerlaw.K', 'simple_source.spectrum.main.Powerlaw.index', 'simple_source.position.l', 'simple_source.position.b'])

By default, the allowed range for the Right Ascension is from 0˚ to 360˚ and allowed declination values range from -90˚ to 90˚. If fitting the source position, it is strongly recommended to restrict the coordinates to be inside 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. For example:

[12]:
simple_source_icrs.position.ra.bounds = ( simple_source_icrs.position.ra.value - 5.0, simple_source_icrs.position.ra.value + 5.0 )
simple_source_icrs.position.dec.bounds = ( simple_source_icrs.position.dec.value - 5.0, simple_source_icrs.position.dec.value + 5.0 )
simple_source_icrs.position.ra.display()
simple_source_icrs.position.dec.display()

Parameter ra = 124.2 [deg] (min_value = 119.2, max_value = 129.2, delta = 12.32, free = True)
Parameter dec = -14.2 [deg] (min_value = -19.2, max_value = -9.2, delta = 1.32, free = True)

“Calling” a point source

Both the point source object itself as well as the compontents are callable functions who take as argument(s) an array of energies and return the differential flux dN/dE at those energies. Energies can be provided with or without units. If no units are provided, energies are assumed to be in keV and fluxes are returned in 1/(cm2 keV s).

[13]:
from astropy import units as u
E = [1, 10, 100]*u.keV

print("Energy in keV:")
print( "With units:", simple_source_icrs(E) )
print( "Without units:", simple_source_icrs(E.value) )

print( "With units:", simple_source_icrs.spectrum.main.shape(E) )
print( "Without units:", simple_source_icrs.spectrum.main.shape(E.value) )

print("")
print("Energy in TeV:")
E_TeV = E.to(u.TeV)
print( "With units:", simple_source_icrs(E_TeV) )
print( "With units:", simple_source_icrs(E_TeV).to(1/u.cm**2/u.TeV/u.s) )
print( "Without units:", simple_source_icrs(E_TeV.value) )



Energy in keV:
With units: [1.00000000e+00 9.77237221e-03 9.54992586e-05] 1 / (keV s cm2)
Without units: [1.00000000e+00 9.77237221e-03 9.54992586e-05]
With units: [1.00000000e+00 9.77237221e-03 9.54992586e-05] 1 / (keV s cm2)
Without units: [1.00000000e+00 9.77237221e-03 9.54992586e-05]

Energy in TeV:
With units: [1.00000000e+00 9.77237221e-03 9.54992586e-05] 1 / (keV s cm2)
With units: [1.00000000e+09 9.77237221e+06 9.54992586e+04] 1 / (TeV s cm2)
Without units: [1.23026877e+18 1.20226443e+16 1.17489755e+14]
[ ]: