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())
21:22:56 WARNING The naima package is not available. Models that depend on it will not be functions.py:47 available
WARNING The GSL library or the pygsl wrapper cannot be loaded. Models that depend on it functions.py:68 will not be available.
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
- ra:
- 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:
- position:
[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)
234.32058209778674 11.365144109043943
123.19999334752814 -13.199993514950487
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'>
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()
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()
“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]
[ ]: