Skip to content

Commit

Permalink
Added a model and a test
Browse files Browse the repository at this point in the history
  • Loading branch information
luca-giunti committed Aug 23, 2019
1 parent 9dac07b commit 0807e34
Show file tree
Hide file tree
Showing 2 changed files with 157 additions and 42 deletions.
128 changes: 86 additions & 42 deletions gammapy/image/models/core.py
View file
Expand Up @@ -2,7 +2,7 @@
import logging
import numpy as np
import astropy.units as u
from astropy.coordinates.angle_utilities import angular_separation
from astropy.coordinates.angle_utilities import angular_separation, position_angle
from astropy.coordinates import Angle, Longitude, Latitude, SkyCoord
from ...utils.fitting import Parameter, Model
from ...maps import Map
Expand Down Expand Up @@ -112,9 +112,10 @@ class SkyGaussian(SkySpatialModel):
.. math::
\phi(\text{lon}, \text{lat}) = N \times \text{exp}\left\{-\frac{1}{2}
\frac{1-\text{cos}\theta}{1-\text{cos}\sigma}\right\}\,,
\frac{1-\text{cos}(r(\text{lon}, \text{lat}))}{1-\text{cos}\sigma}\right\}\,,
where :math:`\theta` is the angular separation between the center of the Gaussian and the evaluation point.
where :math:`r(\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
the sphere:
Expand All @@ -124,11 +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:`\theta` and :math:`\sigma`, this definition reduces to the usual form:
In the limit of small :math:`r(\text{lon}, \text{lat})` 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{\theta^2}{\sigma^2}\right)}
\frac{r(\text{lon}, \text{lat})^2}{\sigma^2}\right)}
Parameters
----------
Expand Down Expand Up @@ -171,53 +173,91 @@ def evaluate(lon, lat, lon_0, lat_0, sigma):
exponent = -0.5 * ((1 - np.cos(sep)) / a)
return u.Quantity(norm.value * np.exp(exponent).value, "sr-1", copy=False)


class SkyGaussianElongated(SkySpatialModel):
r"""Two-dimensional elongated Gaussian model
.. math::
\phi(\text{lon}, \text{lat}) = N \times \text{exp}\left\{-\frac{1}{2}
\frac{1-\text{cos}\theta}{1-\text{cos}\sigma}\right\}\,,
\frac{1-\text{cos}(r(\text{lon}, \text{lat}))}{1-\text{cos}\sigma_{eff}}\right\}\,,
where :math:`\theta` 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
the sphere:
where :math:`r(\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::
N = \frac{1}{4\pi a\left[1-\text{exp}(-1/a)\right]}\,,\,\,\,\,
a = 1-\text{cos}\sigma\,.
\sigma_{eff}(\text{lon}, \text{lat}) = \sqrt{
(\sigma_M \text{sin}(\Delta \theta))^2 +
(\sigma_m \text{cos}(\Delta \theta))^2
},
The normalization factor is in units of :math:`\text{sr}^{-1}`.
In the limit of small :math:`\theta` and :math:`\sigma`, this definition reduces to the usual form:
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
position angle of the evaluation point.
.. math::
\phi(\text{lon}, \text{lat}) = \frac{1}{2\pi\sigma^2} \exp{\left(-\frac{1}{2}
\frac{\theta^2}{\sigma^2}\right)}
**Caveat:** The model is normalized to 1 on the plane, i.e. in small angle approximation:
:math:`N = 1/(2 \pi \sigma_M \sigma_m)`. This means that for huge elongated Gaussians on the sky
this model is not correctly normalized. However, this approximation is perfectly acceptable for the more
common case of models with modest dimensions: indeed, the error introduced by normalizing on the plane
rather than on the sphere is below 0.1\% for Gaussians with radii smaller than ~ 5 deg.
Parameters
----------
lon_0 : `~astropy.coordinates.Longitude`
:math:`\text{lon}_0`
:math:`\text{lon}_0`: `lon` coordinate for the center of the Gaussian.
lat_0 : `~astropy.coordinates.Latitude`
:math:`\text{lat}_0`
sigma : `~astropy.coordinates.Angle`
:math:`\sigma`
:math:`\text{lat}_0`: `lat` coordinate for the center of the Gaussian.
sigma_semi_major : `~astropy.coordinates.Angle`
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`:
Rotation angle of the major semiaxis. The rotation angle increases counter-clockwise
from the North direction.
frame : {"galactic", "icrs"}
Coordinate frame of `lon_0` and `lat_0`.
"""
__slots__ = ["frame", "lon_0", "lat_0", "sigma_semi_major", "e", "theta", "_offset_by"]
def __init__(
self, lon_0, lat_0, sigma_semi_major, e, theta, frame="galactic"
):
try:
from astropy.coordinates.angle_utilities import offset_by
Examples
--------
.. plot::
:include-source:
self._offset_by = offset_by
except ImportError:
raise ImportError("The SkyGaussianElongated model requires astropy>=3.1")
import numpy as np
import matplotlib.pyplot as plt
import astropy.units as u
from astropy.coordinates import Angle
from gammapy.image.models.core import SkyEllipse
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')
coords = m_geom.get_coord()
lon = coords.lon * u.deg
lat = coords.lat * u.deg
vals = model(lon, lat)
skymap = Map.from_geom(m_geom, data=vals.value)
_, ax, _ = skymap.smooth("0.05 deg").plot()
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.vlines(x=2, color='r', linestyle='--', transform=transform, ymin=-5, ymax=5)
ax.text(2.25, 2.45, r"$\theta$", 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"]

def __init__(self, lon_0, lat_0, sigma_semi_major, e, theta, frame="galactic"):
self.frame = frame
self.lon_0 = Parameter(
"lon_0", Longitude(lon_0).wrap_at("180d"), min=-180, max=180
Expand All @@ -230,6 +270,7 @@ def __init__(
super().__init__(
[self.lon_0, self.lat_0, self.sigma_semi_major, self.e, self.theta]
)

@property
def evaluation_radius(self):
r"""Evaluation radius (`~astropy.coordinates.Angle`).
Expand All @@ -240,22 +281,25 @@ def evaluation_radius(self):

def evaluate(self, lon, lat, lon_0, lat_0, sigma_semi_major, e, theta):
"""Evaluate the model (static function)."""
# find the foci of the ellipse corresponding to the 1 sigma isocontour of the Gaussian
c = sigma_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)
sep = angular_separation(lon, lat, lon_0, lat_0)

sep_1 = angular_separation(lon, lat, lon_1, lat_1)
sep_2 = angular_separation(lon, lat, lon_2, lat_2)
theta_0 = position_angle(lon_0, lat_0, lon, lat)
d_theta = theta - theta_0
sigma_semi_minor = Angle(sigma_semi_major * np.sqrt(1 - e ** 2))

# effective distance for model evaluation similar to the symmetric Gaussian
sep = 0.5 * (sep_1 + sep_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
denominator = np.sqrt(a2 + b2)
sigma_eff = sigma_semi_major * sigma_semi_minor / denominator

sigma_semi_minor = sigma_semi_major * np.sqrt(1 - e ** 2)
norm = 1 / (2 * np.pi * sigma_semi_major * sigma_semi_minor)
a = 1.0 - np.cos(sigma_semi_major)

a = 1.0 - np.cos(sigma_eff)
exponent = -0.5 * ((1 - np.cos(sep)) / a)
return u.Quantity(norm.to_value('sr-1') * np.exp(exponent).value, "sr-1", copy=False)
return u.Quantity(
norm.to_value("sr-1") * np.exp(exponent).value, "sr-1", copy=False
)


class SkyDisk(SkySpatialModel):
Expand Down
71 changes: 71 additions & 0 deletions gammapy/image/models/tests/test_core.py
View file
Expand Up @@ -8,6 +8,7 @@
from ..core import (
SkyPointSource,
SkyGaussian,
SkyGaussianElongated,
SkyDisk,
SkyEllipse,
SkyShell,
Expand Down Expand Up @@ -45,6 +46,76 @@ def test_sky_gaussian():
assert_allclose(radius.value, 5 * sigma.value)


def test_sky_gaussian_elongated():
# test the normalization for an elongated Gaussian near the Galactic Plane
m_geom_1 = WcsGeom.create(
binsz=0.05, width=(20, 20), skydir=(2, 2), coordsys="GAL", proj="AIT"
)
coords = m_geom_1.get_coord()
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
)

radius = model_1.evaluation_radius
assert radius.unit == "deg"
assert_allclose(radius.value, 5 * semi_major.value)

# check the ratio between the value at the peak and on the 1-sigma isocontour
semi_major = 4 * u.deg
semi_minor = 2 * u.deg
e = np.sqrt(1 - (semi_minor / semi_major) ** 2)
model_2 = SkyGaussianElongated(0 * u.deg, 0 * u.deg, semi_major, e, 0 * u.deg)
val_0 = model_2(0 * u.deg, 0 * u.deg)
val_major = model_2(0 * u.deg, 4 * u.deg)
val_minor = model_2(2 * u.deg, 0 * u.deg)
assert val_0.unit == "sr-1"
ratio_major = val_0 / val_major
ratio_minor = val_0 / val_minor

assert_allclose(ratio_major, np.exp(0.5))
assert_allclose(ratio_minor, np.exp(0.5))

# check the rotation
model_3 = SkyGaussianElongated(0 * u.deg, 0 * u.deg, semi_major, e, 90 * u.deg)
val_minor_rotated = model_3(0 * u.deg, 2 * u.deg)
ratio_minor_rotated = val_0 / val_minor_rotated
assert_allclose(ratio_minor_rotated, np.exp(0.5))

# compare the normalization of a symmetric Gaussian (ellipse with e=0) and an
# elongated Gaussian with null eccentricity, both defined at the Galactic Pole
m_geom_4 = WcsGeom.create(
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
lon = coords.lon * u.deg
lat = coords.lat * u.deg

semi_major = 5 * u.deg
model_4_el = SkyGaussianElongated(
0 * u.deg, 90 * u.deg, semi_major, 0.0, 0.0 * u.deg
)
model_4_sym = SkyGaussian(0 * u.deg, 90 * u.deg, semi_major)

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)

assert_allclose(int_symmetric, int_elongated, rtol=1e-3)


def test_sky_disk():
r_0 = 2 * u.deg
model = SkyDisk(lon_0="1 deg", lat_0="45 deg", r_0=r_0)
Expand Down

0 comments on commit 0807e34

Please sign in to comment.