From ccafdba8dda1d859274313f1d0181d5b33f16f6d Mon Sep 17 00:00:00 2001 From: luca-giunti Date: Fri, 23 Aug 2019 19:37:22 +0200 Subject: [PATCH] Applied suggested changes --- gammapy/image/models/core.py | 69 +++++++++++++------------ gammapy/image/models/tests/test_core.py | 27 +++++----- 2 files changed, 47 insertions(+), 49 deletions(-) diff --git a/gammapy/image/models/core.py b/gammapy/image/models/core.py index ee406336e0..31b83fbf39 100644 --- a/gammapy/image/models/core.py +++ b/gammapy/image/models/core.py @@ -112,9 +112,9 @@ class SkyGaussian(SkySpatialModel): .. math:: \phi(\text{lon}, \text{lat}) = N \times \text{exp}\left\{-\frac{1}{2} - \frac{1-\text{cos}(r(\text{lon}, \text{lat}))}{1-\text{cos}\sigma}\right\}\,, + \frac{1-\text{cos}(\theta(\text{lon}, \text{lat}))}{1-\text{cos}\sigma}\right\}\,, - where :math:`r(\text{lon}, \text{lat})` is the angular separation between the center of the Gaussian + where :math:`\theta(\text{lon}, \text{lat})` is the angular separation between the center of the Gaussian and the evaluation point. This angle is calculated on the celestial sphere using the function `angular.separation` defined in `astropy.coordinates.angle_utilities`. The Gaussian is normalized to 1 on @@ -125,12 +125,12 @@ class SkyGaussian(SkySpatialModel): a = 1-\text{cos}\sigma\,. The normalization factor is in units of :math:`\text{sr}^{-1}`. - In the limit of small :math:`r(\text{lon}, \text{lat})` and :math:`\sigma`, this definition + In the limit of small :math:`\theta` and :math:`\sigma`, this definition reduces to the usual form: .. math:: \phi(\text{lon}, \text{lat}) = \frac{1}{2\pi\sigma^2} \exp{\left(-\frac{1}{2} - \frac{r(\text{lon}, \text{lat})^2}{\sigma^2}\right)} + \frac{\theta(\text{lon}, \text{lat})^2}{\sigma^2}\right)} Parameters ---------- @@ -179,20 +179,20 @@ class SkyGaussianElongated(SkySpatialModel): .. math:: \phi(\text{lon}, \text{lat}) = N \times \text{exp}\left\{-\frac{1}{2} - \frac{1-\text{cos}(r(\text{lon}, \text{lat}))}{1-\text{cos}\sigma_{eff}}\right\}\,, + \frac{1-\text{cos}(\theta(\text{lon}, \text{lat}))}{1-\text{cos}\sigma_{eff}}\right\}\,, - where :math:`r(\text{lon}, \text{lat})` is the angular separation between the center of + where :math:`\theta(\text{lon}, \text{lat})` is the angular separation between the center of the Gaussian and the evaluation point. The effective radius of the Gaussian, used for the evaluation of the model, is: .. math:: \sigma_{eff}(\text{lon}, \text{lat}) = \sqrt{ - (\sigma_M \text{sin}(\Delta \theta))^2 + - (\sigma_m \text{cos}(\Delta \theta))^2 + (\sigma_M \text{sin}(\Delta \phi))^2 + + (\sigma_m \text{cos}(\Delta \phi))^2 }, where :math:`\sigma_M` (:math:`\sigma_m`) is the major (minor) semiaxis of the Gaussian, and - :math:`\Delta \theta` is the difference between `theta`, the position angle of the Gaussian, and the + :math:`\Delta \phi` is the difference between `phi`, the position angle of the Gaussian, and the position angle of the evaluation point. **Caveat:** The model is normalized to 1 on the plane, i.e. in small angle approximation: @@ -211,8 +211,8 @@ class SkyGaussianElongated(SkySpatialModel): Length of the major semiaxis of the Gaussian, in angular units. e : `float` Eccentricity of the Gaussian (:math:`0< e< 1`). - theta : `~astropy.coordinates.Angle` - :math:`\theta`: + phi : `~astropy.coordinates.Angle` + :math:`\phi`: Rotation angle of the major semiaxis. The rotation angle increases counter-clockwise from the North direction. frame : {"galactic", "icrs"} @@ -232,8 +232,8 @@ class SkyGaussianElongated(SkySpatialModel): from gammapy.maps import Map, WcsGeom m_geom = WcsGeom.create(binsz=0.01, width=(5, 5), skydir=(2, 2), coordsys="GAL", proj="AIT") - theta = Angle("30 deg") - model = SkyGaussianElongated("2 deg", "2 deg", "1 deg", 0.7, theta, frame='galactic') + phi = Angle("30 deg") + model = SkyGaussianElongated("2 deg", "2 deg", "1 deg", 0.7, phi, frame='galactic') coords = m_geom.get_coord() @@ -247,17 +247,17 @@ class SkyGaussianElongated(SkySpatialModel): transform = ax.get_transform('galactic') ax.scatter(2, 2, transform=transform, s=20, edgecolor='red', facecolor='red') ax.text(1.5, 1.85, r"$(l_0, b_0)$", transform=transform, ha="center") - ax.plot([2, 2 + np.sin(theta)], [2, 2 + np.cos(theta)], color="r", transform=transform) + ax.plot([2, 2 + np.sin(phi)], [2, 2 + np.cos(phi)], color="r", transform=transform) ax.vlines(x=2, color='r', linestyle='--', transform=transform, ymin=-5, ymax=5) - ax.text(2.25, 2.45, r"$\theta$", transform=transform); + ax.text(2.25, 2.45, r"$\phi$", transform=transform); ax.contour(skymap.data, cmap='coolwarm', levels=10, alpha=0.6) plt.show() """ - __slots__ = ["frame", "lon_0", "lat_0", "sigma_semi_major", "e", "theta"] + __slots__ = ["frame", "lon_0", "lat_0", "sigma_semi_major", "e", "phi"] - def __init__(self, lon_0, lat_0, sigma_semi_major, e, theta, frame="galactic"): + def __init__(self, lon_0, lat_0, sigma_semi_major, e, phi, frame="galactic"): self.frame = frame self.lon_0 = Parameter( "lon_0", Longitude(lon_0).wrap_at("180d"), min=-180, max=180 @@ -265,10 +265,10 @@ def __init__(self, lon_0, lat_0, sigma_semi_major, e, theta, frame="galactic"): self.lat_0 = Parameter("lat_0", Latitude(lat_0), min=-90, max=90) self.sigma_semi_major = Parameter("sigma_semi_major", Angle(sigma_semi_major)) self.e = Parameter("e", e, min=0, max=1) - self.theta = Parameter("theta", Angle(theta)) + self.phi = Parameter("phi", Angle(phi)) super().__init__( - [self.lon_0, self.lat_0, self.sigma_semi_major, self.e, self.theta] + [self.lon_0, self.lat_0, self.sigma_semi_major, self.e, self.phi] ) @property @@ -279,17 +279,18 @@ def evaluation_radius(self): """ return 5 * self.parameters["sigma_semi_major"].quantity - def evaluate(self, lon, lat, lon_0, lat_0, sigma_semi_major, e, theta): + @staticmethod + def evaluate(lon, lat, lon_0, lat_0, sigma_semi_major, e, phi): """Evaluate the model (static function).""" sep = angular_separation(lon, lat, lon_0, lat_0) - theta_0 = position_angle(lon_0, lat_0, lon, lat) - d_theta = theta - theta_0 + phi_0 = position_angle(lon_0, lat_0, lon, lat) + d_phi = phi - phi_0 sigma_semi_minor = Angle(sigma_semi_major * np.sqrt(1 - e ** 2)) # Effective radius, used for model evaluation as in the symmetric case - a2 = (sigma_semi_major * np.sin(d_theta)) ** 2 - b2 = (sigma_semi_minor * np.cos(d_theta)) ** 2 + a2 = (sigma_semi_major * np.sin(d_phi)) ** 2 + b2 = (sigma_semi_minor * np.cos(d_phi)) ** 2 denominator = np.sqrt(a2 + b2) sigma_eff = sigma_semi_major * sigma_semi_minor / denominator @@ -431,8 +432,8 @@ class SkyEllipse(SkySpatialModel): :math:`a`: length of the major semiaxis, in angular units. e : `float` Eccentricity of the ellipse (:math:`0< e< 1`). - theta : `~astropy.coordinates.Angle` - :math:`\theta`: + phi : `~astropy.coordinates.Angle` + :math:`\phi`: Rotation angle of the major semiaxis. The rotation angle increases counter-clockwise from the North direction. edge : `~astropy.coordinates.Angle` @@ -469,15 +470,15 @@ class SkyEllipse(SkySpatialModel): ax.text(1.7, 1.85, r"$(l_0, b_0)$", transform=transform, ha="center") ax.plot([2, 2 + np.sin(np.pi / 6)], [2, 2 + np.cos(np.pi / 6)], color="r", transform=transform) ax.vlines(x=2, color='r', linestyle='--', transform=transform, ymin=0, ymax=5) - ax.text(2.15, 2.3, r"$\theta$", transform=transform); + ax.text(2.15, 2.3, r"$\phi$", transform=transform); plt.show() """ - __slots__ = ["frame", "lon_0", "lat_0", "semi_major", "e", "theta", "_offset_by"] + __slots__ = ["frame", "lon_0", "lat_0", "semi_major", "e", "phi", "_offset_by"] def __init__( - self, lon_0, lat_0, semi_major, e, theta, edge="0.01 deg", frame="galactic" + self, lon_0, lat_0, semi_major, e, phi, edge="0.01 deg", frame="galactic" ): try: from astropy.coordinates.angle_utilities import offset_by @@ -493,11 +494,11 @@ def __init__( self.lat_0 = Parameter("lat_0", Latitude(lat_0), min=-90, max=90) 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.phi = Parameter("phi", Angle(phi)) 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, self.edge] + [self.lon_0, self.lat_0, self.semi_major, self.e, self.phi, self.edge] ) @property @@ -525,12 +526,12 @@ 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, edge): + def evaluate(self, lon, lat, lon_0, lat_0, semi_major, e, phi, edge): """Evaluate the model (static function).""" # find the foci of the ellipse c = semi_major * e - lon_1, lat_1 = self._offset_by(lon_0, lat_0, theta, c) - lon_2, lat_2 = self._offset_by(lon_0, lat_0, 180 * u.deg + theta, c) + lon_1, lat_1 = self._offset_by(lon_0, lat_0, phi, c) + lon_2, lat_2 = self._offset_by(lon_0, lat_0, 180 * u.deg + phi, c) sep_1 = angular_separation(lon, lat, lon_1, lat_1) sep_2 = angular_separation(lon, lat, lon_2, lat_2) diff --git a/gammapy/image/models/tests/test_core.py b/gammapy/image/models/tests/test_core.py index 4a031e82e3..d076a524e4 100644 --- a/gammapy/image/models/tests/test_core.py +++ b/gammapy/image/models/tests/test_core.py @@ -52,15 +52,15 @@ def test_sky_gaussian_elongated(): binsz=0.05, width=(20, 20), skydir=(2, 2), coordsys="GAL", proj="AIT" ) coords = m_geom_1.get_coord() + angles = m_geom_1.solid_angle() lon = coords.lon * u.deg lat = coords.lat * u.deg semi_major = 3 * u.deg model_1 = SkyGaussianElongated(2 * u.deg, 2 * u.deg, semi_major, 0.8, 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.0e-3 + np.sum(vals_1 * angles), 1, rtol=1.0e-3 ) radius = model_1.evaluation_radius @@ -94,7 +94,7 @@ def test_sky_gaussian_elongated(): binsz=0.05, width=(25, 25), skydir=(0, 90), coordsys="GAL", proj="AIT" ) coords = m_geom_4.get_coord() - angles = m_geom_4.solid_angle().value + angles = m_geom_4.solid_angle() lon = coords.lon * u.deg lat = coords.lat * u.deg @@ -107,11 +107,8 @@ def test_sky_gaussian_elongated(): vals_4_el = model_4_el(lon, lat) vals_4_sym = model_4_sym(lon, lat) - mymap_4_el = Map.from_geom(m_geom_4, data=vals_4_el.value) - mymap_4_sym = Map.from_geom(m_geom_4, data=vals_4_sym.value) - - int_elongated = np.sum(mymap_4_el.quantity * angles) - int_symmetric = np.sum(mymap_4_sym.quantity * angles) + int_elongated = np.sum(vals_4_el * angles) + int_symmetric = np.sum(vals_4_sym * angles) assert_allclose(int_symmetric, int_elongated, rtol=1e-3) @@ -152,15 +149,15 @@ def test_sky_ellipse(): binsz=0.015, width=(20, 20), skydir=(2, 2), coordsys="GAL", proj="AIT" ) coords = m_geom_1.get_coord() + angles = m_geom_1.solid_angle() lon = coords.lon * u.deg lat = coords.lat * 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.0e-3 + np.sum(vals_1 * angles), 1, rtol=1.0e-3 ) radius = model_1.evaluation_radius @@ -181,24 +178,24 @@ def test_sky_ellipse(): binsz=0.1, width=(6, 6), skydir=(0, 90), coordsys="GAL", proj="AIT" ) coords = m_geom_2.get_coord() + angles = m_geom_2.solid_angle() + lon = coords.lon * u.deg lat = coords.lat * 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) - 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(vals_2 * angles), + np.sum(vals_disk * angles), ) @@ -206,7 +203,7 @@ 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" + lon_0="0 deg", lat_0="0 deg", semi_major=r_0, e=0.5, phi="0 deg" ) value_center = model(0 * u.deg, 0 * u.deg) value_edge = model(0 * u.deg, r_0)