Merge pull request #2065 from luca-giunti/Define-evaluation_radius-fo…

…r-SkyEllipse

Added evaluation radius for the SkyEllipse model, and a test

gammapy/image/models/core.py

@@ -287,6 +287,7 @@ class SkyEllipse(SkySpatialModel):
     def __init__(self, lon_0, lat_0, semi_major, e, theta):
         try:
             from astropy.coordinates.angle_utilities import offset_by
+
             self._offset_by = offset_by
         except ImportError:
             raise ImportError("The SkyEllipse model requires astropy>=3.1")
@@ -301,6 +302,20 @@ def __init__(self, lon_0, lat_0, semi_major, e, theta):
             ]
         )
 
+    @property
+    def evaluation_radius(self):
+        r"""Returns the effective radius of the sky region where the model evaluates to non-zero.
+        For an elliptical source, we fix it to the length of the semi-major axis.
+
+        Returns
+        -------
+        radius : `~astropy.coordinates.Angle`
+            Radius in angular units
+
+        """
+        radius = self.parameters["semi_major"].quantity
+        return radius
+
     def compute_norm(semi_major, e):
         """Compute the normalization factor."""
         semi_minor = semi_major * np.sqrt(1 - e ** 2)
@@ -385,10 +400,7 @@ def evaluation_radius(self):
             Radius in angular units
 
         """
-        radius = (
-            self.parameters["radius"].quantity
-            + self.parameters["width"].quantity
-        )
+        radius = self.parameters["radius"].quantity + self.parameters["width"].quantity
         return radius
 
     @staticmethod

gammapy/image/models/tests/test_core.py

@@ -63,41 +63,46 @@ def test_sky_ellipse():
     coords = m_geom_1.get_coord()
     lon = coords.lon * u.deg
     lat = coords.lat * u.deg
-    model_1 = SkyEllipse(2 * u.deg, 2 * u.deg, 10 * u.deg, 0.4, 30 * u.deg)
+    semi_major = 10 * u.deg
+    model_1 = SkyEllipse(2 * u.deg, 2 * u.deg, semi_major, 0.4, 30 * u.deg)
     vals_1 = model_1(lon, lat)
     assert vals_1.unit == "sr-1"
     mymap_1 = Map.from_geom(m_geom_1, data=vals_1.value)
     assert_allclose(
-        np.sum(mymap_1.quantity * u.sr ** -1 * m_geom_1.solid_angle()), 1, rtol=1.e-3
+        np.sum(mymap_1.quantity * u.sr ** -1 * m_geom_1.solid_angle()), 1, rtol=1.0e-3
     )
 
+    radius = model_1.evaluation_radius
+    assert radius.unit == "deg"
+    assert_allclose(radius.value, semi_major.value)
+
     # test the normalization for a disk (ellipse with e=0) at the Galactic Pole,
     # both analytically and comparing with the SkyDisk model
     m_geom_2 = WcsGeom.create(
-        binsz=0.1, width=(6,6), skydir=(0, 90), coordsys="GAL", proj="AIT"
+        binsz=0.1, width=(6, 6), skydir=(0, 90), coordsys="GAL", proj="AIT"
     )
     coords = m_geom_2.get_coord()
     lon = coords.lon * u.deg
     lat = coords.lat * u.deg
 
-    model_2 = SkyEllipse(0 * u.deg, 90 * u.deg, 5 * u.deg, 0.0, 0.0 * u.deg)
+    semi_major = 5 * u.deg
+    model_2 = SkyEllipse(0 * u.deg, 90 * u.deg, semi_major, 0.0, 0.0 * u.deg)
     vals_2 = model_2(lon, lat)
     mymap_2 = Map.from_geom(m_geom_2, data=vals_2.value)
 
     disk = SkyDisk(lon_0="0 deg", lat_0="90 deg", r_0="5 deg")
-    vals_disk = disk(lon,lat)
+    vals_disk = disk(lon, lat)
     mymap_disk = Map.from_geom(m_geom_2, data=vals_disk.value)
 
     solid_angle = 2 * np.pi * (1 - np.cos(5 * u.deg))
     assert_allclose(np.max(vals_2).value * solid_angle, 1)
 
     assert_allclose(
         np.sum(mymap_2.quantity * u.sr ** -1 * m_geom_2.solid_angle()),
-        np.sum(mymap_disk.quantity * u.sr ** -1 * m_geom_2.solid_angle())
+        np.sum(mymap_disk.quantity * u.sr ** -1 * m_geom_2.solid_angle()),
     )
 
 
-
 def test_sky_shell():
     width = 2 * u.deg
     rad = 2 * u.deg