Merge pull request #2141 from adonath/smooth_disk_ellipse_edge

Implement SkyDisk and SkyEllipse edge parameter

gammapy/image/models/core.py

@@ -8,6 +8,7 @@
 from ...maps import Map
 from scipy.integrate import quad
 
+
 __all__ = [
     "SkySpatialModel",
     "SkyPointSource",
@@ -23,6 +24,11 @@
 log = logging.getLogger(__name__)
 
 
+def smooth_edge(x, width):
+    value = (x / width).to_value("")
+    return 0.5 - np.clip(value, -0.5, 0.5)
+
+
 class SkySpatialModel(Model):
     """Sky spatial model base class."""
 
@@ -187,7 +193,7 @@ class SkyDisk(SkySpatialModel):
                     0 & \text{for } \theta > r_0
                 \end{cases}
 
-    where :math:`\theta` is the sky separation
+    where :math:`\theta` is the sky separation.
 
     Parameters
     ----------
@@ -197,21 +203,24 @@ class SkyDisk(SkySpatialModel):
         :math:`lat_0`
     r_0 : `~astropy.coordinates.Angle`
         :math:`r_0`
+    edge : `~astropy.coordinates.Angle`
+        Width of the edge.
     frame : {"galactic", "icrs"}
         Coordinate frame of `lon_0` and `lat_0`.
     """
 
     __slots__ = ["frame", "lon_0", "lat_0", "r_0"]
 
-    def __init__(self, lon_0, lat_0, r_0, frame="galactic"):
+    def __init__(self, lon_0, lat_0, r_0, edge="0.01 deg", frame="galactic"):
         self.frame = frame
         self.lon_0 = Parameter(
             "lon_0", Longitude(lon_0).wrap_at("180d"), min=-180, max=180
         )
         self.lat_0 = Parameter("lat_0", Latitude(lat_0), min=-90, max=90)
         self.r_0 = Parameter("r_0", Angle(r_0))
+        self.edge = Parameter("edge", Angle(edge), min=0.01, frozen=True)
 
-        super().__init__([self.lon_0, self.lat_0, self.r_0])
+        super().__init__([self.lon_0, self.lat_0, self.r_0, self.edge])
 
     @property
     def evaluation_radius(self):
@@ -226,13 +235,15 @@ def evaluation_radius(self):
         return self.parameters["r_0"].quantity
 
     @staticmethod
-    def evaluate(lon, lat, lon_0, lat_0, r_0):
+    def evaluate(lon, lat, lon_0, lat_0, r_0, edge):
         """Evaluate the model (static function)."""
         sep = angular_separation(lon, lat, lon_0, lat_0)
 
         # Surface area of a spherical cap, see https://en.wikipedia.org/wiki/Spherical_cap
         norm = 1.0 / (2 * np.pi * (1 - np.cos(r_0)))
-        return u.Quantity(norm.value * (sep <= r_0), "sr-1", copy=False)
+
+        in_disk = smooth_edge(sep - r_0, edge)
+        return u.Quantity(norm.value * in_disk, "sr-1", copy=False)
 
 
 class SkyEllipse(SkySpatialModel):
@@ -270,6 +281,8 @@ class SkyEllipse(SkySpatialModel):
         :math:`\theta`:
         Rotation angle of the major semiaxis.  The rotation angle increases clockwise
         (i.e., East of North) from the positive `lon` axis.
+    edge : `~astropy.coordinates.Angle`
+        Width of the edge.
     frame : {"galactic", "icrs"}
         Coordinate frame of `lon_0` and `lat_0`.
 
@@ -306,7 +319,7 @@ class SkyEllipse(SkySpatialModel):
 
     __slots__ = ["frame", "lon_0", "lat_0", "semi_major", "e", "theta", "_offset_by"]
 
-    def __init__(self, lon_0, lat_0, semi_major, e, theta, frame="galactic"):
+    def __init__(self, lon_0, lat_0, semi_major, e, theta, edge="0.01 deg", frame="galactic"):
         try:
             from astropy.coordinates.angle_utilities import offset_by
 
@@ -322,8 +335,9 @@ def __init__(self, lon_0, lat_0, semi_major, e, theta, frame="galactic"):
         self.semi_major = Parameter("semi_major", Angle(semi_major))
         self.e = Parameter("e", e, min=0, max=1)
         self.theta = Parameter("theta", Angle(theta))
+        self.edge = Parameter("edge", Angle(edge), frozen=True, min=0.01)
 
-        super().__init__([self.lon_0, self.lat_0, self.semi_major, self.e, self.theta])
+        super().__init__([self.lon_0, self.lat_0, self.semi_major, self.e, self.theta, self.edge])
 
     @property
     def evaluation_radius(self):
@@ -356,7 +370,7 @@ def integral_fcn(x, a, b):
             2 * quad(lambda x: integral_fcn(x, semi_major, semi_minor), 0, np.pi)[0]
         ) ** -1
 
-    def evaluate(self, lon, lat, lon_0, lat_0, semi_major, e, theta):
+    def evaluate(self, lon, lat, lon_0, lat_0, semi_major, e, theta, edge):
         """Evaluate the model (static function)."""
         # find the foci of the ellipse
         c = semi_major * e
@@ -365,7 +379,8 @@ def evaluate(self, lon, lat, lon_0, lat_0, semi_major, e, theta):
 
         sep_1 = angular_separation(lon, lat, lon_1, lat_1)
         sep_2 = angular_separation(lon, lat, lon_2, lat_2)
-        in_ellipse = sep_1 + sep_2 <= 2 * semi_major
+
+        in_ellipse = smooth_edge(sep_1 + sep_2 - 2 * semi_major, edge)
 
         norm = SkyEllipse.compute_norm(semi_major, e)
         return u.Quantity(norm * in_ellipse, "sr-1", copy=False)

gammapy/image/models/tests/test_core.py

@@ -58,6 +58,13 @@ def test_sky_disk():
     assert radius.unit == "deg"
     assert_allclose(radius.value, r_0.value)
 
+def test_sky_disk_edge():
+    r_0 = 2 * u.deg
+    model = SkyDisk(lon_0="0 deg", lat_0="0 deg", r_0=r_0)
+    value_center = model(0 * u.deg, 0 * u.deg)
+    value_edge = model(r_0, 0 * u.deg)
+    assert_allclose((value_edge / value_center).to_value(""), 0.5)
+
 
 def test_sky_ellipse():
     pytest.importorskip("astropy", minversion="3.1.1")
@@ -108,6 +115,15 @@ def test_sky_ellipse():
     )
 
 
+def test_sky_ellipse_edge():
+    pytest.importorskip("astropy", minversion="3.1.1")
+    r_0 = 2 * u.deg
+    model = SkyEllipse(lon_0="0 deg", lat_0="0 deg", semi_major=r_0, e=0.5, theta="0 deg")
+    value_center = model(0 * u.deg, 0 * u.deg)
+    value_edge = model(r_0, 0 * u.deg)
+    assert_allclose((value_edge / value_center).to_value(""), 0.5)
+
+
 def test_sky_shell():
     width = 2 * u.deg
     rad = 2 * u.deg