In [1]:
from astromodels import *
In [2]:
# Using J2000 R.A. and Dec (ICRS), which is the default coordinate system:
simple_source_icrs = PointSource('simple_source', ra=123.2, dec=-13.2, spectral_shape=powerlaw())
We can also use Galactic coordinates:
In [3]:
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")
In [4]:
simple_source_icrs.display()
# or print(simple_source_icrs) for a text-only representation
As you can see we have created a point source with one component (see below) automatically named "main", with a power law spectrum, at the specified position.
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:
In [5]:
l = simple_source_icrs.position.get_l()
b = simple_source_icrs.position.get_b()
ra = simple_source_gal.position.get_ra()
dec = simple_source_gal.position.get_dec()
type(ra)
Out[5]:
The get_ra, get_dec, get_l and get_b return either a Latitude or Longitude object of astropy.coordinates, from which you can obtain all formats for the coordinates, like:
In [6]:
# Decimal R.A.
print("Decimal R.A. is %s" % ra.deg)
print("Sexadecimal R.A. is %.0f:%.0f:%s" % (ra.dms.d, ra.dms.m, ra.dms.s))
For more control on the output and many more options, such as transform to local frames or other equinoxes, compute distances between points in the sky, and so on, you can obtain an instance of astropy.coordinates.SkyCoord by using the sky_coord property of the position object:
In [7]:
# Refer to the documentation of the astropy.coordinates.SkyCoord class:
# http://docs.astropy.org/en/stable/coordinates/
# for all available options.
sky_coord_instance = simple_source_icrs.position.sky_coord
# Now you can transform to another reference use transform_to.
# Here for example we compute the altitude of our source for HAWC at 2 am on 2013-07-01
from astropy.coordinates import SkyCoord, EarthLocation, AltAz
from astropy.time import Time
hawc_site = EarthLocation(lat=19*u.deg, lon=97.3*u.deg, height=4100*u.m)
utcoffset = -5*u.hour # Hour at HAWC is CDT, which is UTC - 5 hours
time = Time('2013-7-01 02:00:00') - utcoffset
src_altaz = sky_coord_instance.transform_to(AltAz(obstime=time,location=hawc_site))
print("Source Altitude at HAWC : {0.alt:.5}".format(src_altaz))
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:
In [8]:
# These will return two Parameter instances corresponding to the parameters ra and dec
# NOT the corresponding floating point numbers:
parameter_ra = simple_source_icrs.position.ra
parameter_dec = simple_source_icrs.position.dec
# This would instead throw AttributeError, since simple_source_icrs was instanced using
# R.A. and Dec. and hence does not have the l,b parameters:
# error = simple_source_icrs.position.l
# error = simple_source_icrs.position.b
# Similarly this will throw AttributeError, because simple_source_gal was instanced using
# Galactic coordinates:
# error = simple_source_gal.position.ra
# error = simple_source_gal.position.dec
# In all cases, independently on how the source was instanced, you can obtain the coordinates
# as normal floating point numbers using:
ra1 = simple_source_icrs.position.get_ra().value
dec1 = simple_source_icrs.position.get_dec().value
l1 = simple_source_icrs.position.get_l().value
b1 = simple_source_icrs.position.get_b().value
ra2 = simple_source_gal.position.get_ra().value
dec2 = simple_source_gal.position.get_dec().value
l2 = simple_source_gal.position.get_l().value
b2 = simple_source_gal.position.get_b().value
A multi-component source is a point source which has different spectral components. For example, in a Gamma-Ray Burst you can have a Synchrotron component and a Inverse Compton component, which come from different zones and are described by different spectra. Depending on the needs of your analysis, you might model this situation using a single component constituted by the sum of the two spectra, or you might want to model them independently. The latter choice allows you to measure for instance the fluxes from the two components independently. Also, each components has its own polarization, which can be useful when studying polarized sources (to be implemented). Representing a source with more than one component is easy in astromodels:
In [9]:
# Create the two different components
#(of course the shape can be any function, or any composite function)
component1 = SpectralComponent('synchrotron',shape=powerlaw())
component2 = SpectralComponent('IC',shape=powerlaw())
# Create a multi-component source
multicomp_source = PointSource('multicomp_source', ra=123.2, dec=-13.2, components=[component1,component2])
multicomp_source.display()
Starting from the source instance you can modify any of its components, or its position, in a straightforward way:
In [10]:
# Change position
multicomp_source.position.ra = 124.5
multicomp_source.position.dec = -11.5
# Change values for the parameters
multicomp_source.spectrum.synchrotron.powerlaw.logK = -1.2
multicomp_source.spectrum.IC.powerlaw.index = -1.0
# To avoid having to write that much, you can create a "shortcut" for a function
po = multicomp_source.spectrum.synchrotron.powerlaw
# Now you can modify its parameters more easily
# (see "Creating and modifying functions" for more info on what you can to with a parameter)
po.K = 1e-5
# Change the minimum using explicit units
po.K.min_value = 1e-6 * 1 / (u.MeV * u.cm**2 * u.s)
# GOTCHA
# Creating a shortcut directly to the parameter will not work:
# p1 = multicomp_source.spectrum.synchrotron.powerlaw.logK
# p1 = -1.3 # this does NOT change the value of logK, but instead assign -1.3 to p1 (i.e., destroy the shortcut)
# However you can change the value of p1 like this:
# p1.value = -1.3 # This will work
multicomp_source.display()
In [ ]: