Skip to content

SpatialTemplate_2D, ExtendedSource and default units #238

@israelmcmc

Description

@israelmcmc

The default solid angles unit for ExtendedSource are 1/deg^2. This is mentioned in the documentation, and in the code they are set here, based on the default units defined here.

SpatialTemplate_2D template ,however, expects a template normalized based on 1/sr:

# test if the map is normalized as expected
area = wcs.utils.proj_plane_pixel_area(self._wcs)
dOmega = (area * u.deg * u.deg).to(u.sr).value
total = self._map.sum() * dOmega
if not np.isclose(total, 1, rtol=1e-2):
log.warning(
"2D template read from {} is normalized to {} (expected: 1)".format(
self._fitsfile, total
)
)

Later on in the code, you can see that SpatialTemplate_2D is only sampling from self._map, which hasn't been converted to 1/deg2.

out[iz] = self._map[Ypix[iz], Xpix[iz]]
return np.multiply(K, out)

The K at this point, the K value has the right units, but the values do not match. e.g.

from astromodels.functions import SpatialTemplate_2D
from astromodels import ExtendedSource, Powerlaw

lon = 275.9326171875
lat = -12.941631985423356
energy = 500.

# File from https://astropy.stsci.edu/data/photometry/spitzer_example_image.fits
shape = SpatialTemplate_2D(fits_file='spitzer_example_image.fits')
spectrum = Powerlaw()

print(shape.K.unit) # = None
print(spectrum(energy) * shape(lon, lat)) # = 0.0001474072331793441

source = ExtendedSource('simple_source', spectral_shape=spectrum, spatial_shape = shape )

print(shape.K.unit) # = 1 / deg2
print(source(lon, lat, energy)) # = 0.0001474072331793441

This causes that a ExtendedModel using SpatialTemplate_2D returns 1/sr by default.

Is this a bug or are we missing something? At the very least, it looks error-prone.

For context, this came up here: cositools/cosipy#438

Metadata

Metadata

Assignees

Labels

No labels
No labels

Type

No type

Projects

No projects

Milestone

No milestone

Relationships

None yet

Development

No branches or pull requests

Issue actions