Fix Disk_on_sphere.evaluate() call with units

[israelmcmc]

@PreisTo pointed out in the #238 that the following code using dev (d756924) returns an unexpected result:

from astromodels.sources import ExtendedSource
from astromodels.functions import Disk_on_sphere, Powerlaw
import astropy.units as u

pl = Powerlaw()
disk = Disk_on_sphere()
es = ExtendedSource("test",spectral_shape = pl, spatial_shape = disk)

disk.radius = 1*u.deg
pl.K = 1*u.Unit("keV-1 cm-2 s-1")
pl.piv = 1*u.keV

ra = [0]*u.deg
dec = [0]*u.deg
E = [1,10,100]*u.keV

print(es(ra,dec,E)) # [1.04494972e+03 1.02116376e+01 9.97919231e-02]  1 / (keV s cm2 deg2)

It should be either
[1.04494972e+03 1.02116376e+01 9.97919231e-02] 1 / (keV s cm2 sr)
or
[3.18309886e-01 3.11064269e-03 3.03983581e-05] 1 / (keV s cm2 deg2)

The core of the bug is that the Disk_on_sphere definition::

$$ f(\vec{x}) = \left(\frac{180}{\pi}\right)^2 \frac{1}{\pi~({\rm radius})^2} $$

implies a conversion from a radius in degrees to rad, such that the result is in 1/sr. Because the extra factor (180/pi)^2 doesn't have units --i.e. deg2/sr-- astropy cannot get the right result by itself.

I modified the Disk_on_sphere.evaluate() to recognize if the inputs have unit and perform the necessary conversion.

Now the results are:

print(es(ra,dec,E)) # [1.04494972e+03 1.02116376e+01 9.97919231e-02]  1 / (keV s cm2 deg2)
print(es([0],[0],[1,10,100])) # [1.04494972e+03 1.02116376e+01 9.97919231e-02]  (Implicit  1 / (keV s cm2 sr) )

This fixes the bug when called with units, while keeping it backward compatible when called without units -- I think the (180/pi)^2 shouldn't even be there, but I wanted to keep thing backwards compatible. Maybe this could be fixed in a future major release.

I also did various tricks to limit the overhead from allowing input with and without units, as well as the overhead when actually calling with with unit. This is the benchmark:

And this the new version calling it without units (<10% overhead) and using units (almost 2x, mostly from Quantity.__new__()). It's a small overhead on an already cheap operation, and it will be negligible for a call with multiple values.

I think this shows that Function.__call__() shouldn't be the one handling the units for performance reasons, but instead leaving it to the child classes to do it appropriately. I think it's less error prone and it's just as fast.

Lastly, while this PR fixes the individual issue that @PreisTo brought up, the original issue in #238 remains: ExtendedSource is returning units of 1/sr by default, while the documentation and default units say it should be 1/deg2.

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 82.76%. Comparing base (6ebdf11) to head (86db75f).
⚠️ Report is 250 commits behind head on dev.

Additional details and impacted files

@@            Coverage Diff             @@
##              dev     #242      +/-   ##
==========================================
+ Coverage   82.09%   82.76%   +0.67%     
==========================================
  Files          66       59       -7     
  Lines        8204     8049     -155     
==========================================
- Hits         6735     6662      -73     
+ Misses       1469     1387      -82     

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.