1
|
|
|
# Licensed under a 3-clause BSD style license - see LICENSE.rst |
2
|
|
|
import pytest |
3
|
|
|
from numpy.testing import assert_allclose |
4
|
|
|
import numpy as np |
5
|
|
|
import astropy.units as u |
6
|
|
|
from ....maps import Map, WcsGeom |
7
|
|
|
from ....utils.testing import requires_data |
8
|
|
|
from ..core import ( |
9
|
|
|
SkyPointSource, |
10
|
|
|
SkyGaussian, |
11
|
|
|
SkyDisk, |
12
|
|
|
SkyEllipse, |
13
|
|
|
SkyShell, |
14
|
|
|
SkyDiffuseConstant, |
15
|
|
|
SkyDiffuseMap, |
16
|
|
|
) |
17
|
|
|
|
18
|
|
|
|
19
|
|
|
def test_sky_point_source(): |
20
|
|
|
model = SkyPointSource(lon_0="2.5 deg", lat_0="2.5 deg") |
21
|
|
|
lat, lon = np.mgrid[0:6, 0:6] * u.deg |
22
|
|
|
val = model(lon, lat) |
23
|
|
|
assert val.unit == "deg-2" |
24
|
|
|
assert_allclose(val.sum().value, 1) |
25
|
|
|
radius = model.evaluation_radius |
26
|
|
|
assert radius.unit == "deg" |
27
|
|
|
assert_allclose(radius.value, 0) |
28
|
|
|
assert model.frame == "galactic" |
29
|
|
|
|
30
|
|
|
assert_allclose(model.position.l.deg, 2.5) |
31
|
|
|
assert_allclose(model.position.b.deg, 2.5) |
32
|
|
|
|
33
|
|
|
|
34
|
|
|
def test_sky_gaussian(): |
35
|
|
|
sigma = 1 * u.deg |
36
|
|
|
model = SkyGaussian(lon_0="5 deg", lat_0="15 deg", sigma=sigma) |
37
|
|
|
assert model.parameters["sigma"].min == 0 |
38
|
|
|
val_0 = model(5 * u.deg, 15 * u.deg) |
39
|
|
|
val_sigma = model(5 * u.deg, 16 * u.deg) |
40
|
|
|
assert val_0.unit == "sr-1" |
41
|
|
|
ratio = val_0 / val_sigma |
42
|
|
|
assert_allclose(ratio, np.exp(0.5)) |
43
|
|
|
radius = model.evaluation_radius |
44
|
|
|
assert radius.unit == "deg" |
45
|
|
|
assert_allclose(radius.value, 5 * sigma.value) |
46
|
|
|
|
47
|
|
|
|
48
|
|
|
def test_sky_disk(): |
49
|
|
|
r_0 = 2 * u.deg |
50
|
|
|
model = SkyDisk(lon_0="1 deg", lat_0="45 deg", r_0=r_0) |
51
|
|
|
lon = [1, 5, 359] * u.deg |
52
|
|
|
lat = 46 * u.deg |
53
|
|
|
val = model(lon, lat) |
54
|
|
|
assert val.unit == "sr-1" |
55
|
|
|
desired = [261.263956, 0, 261.263956] |
56
|
|
|
assert_allclose(val.value, desired) |
57
|
|
|
radius = model.evaluation_radius |
58
|
|
|
assert radius.unit == "deg" |
59
|
|
|
assert_allclose(radius.value, r_0.value) |
60
|
|
|
|
61
|
|
|
def test_sky_disk_edge(): |
62
|
|
|
r_0 = 2 * u.deg |
63
|
|
|
model = SkyDisk(lon_0="0 deg", lat_0="0 deg", r_0=r_0) |
64
|
|
|
value_center = model(0 * u.deg, 0 * u.deg) |
65
|
|
|
value_edge = model(r_0, 0 * u.deg) |
66
|
|
|
assert_allclose((value_edge / value_center).to_value(""), 0.5) |
67
|
|
|
|
68
|
|
|
|
69
|
|
|
def test_sky_ellipse(): |
70
|
|
|
pytest.importorskip("astropy", minversion="3.1.1") |
71
|
|
|
# test the normalization for an elongated ellipse near the Galactic Plane |
72
|
|
|
m_geom_1 = WcsGeom.create( |
73
|
|
|
binsz=0.015, width=(20, 20), skydir=(2, 2), coordsys="GAL", proj="AIT" |
74
|
|
|
) |
75
|
|
|
coords = m_geom_1.get_coord() |
76
|
|
|
lon = coords.lon * u.deg |
77
|
|
|
lat = coords.lat * u.deg |
78
|
|
|
semi_major = 10 * u.deg |
79
|
|
|
model_1 = SkyEllipse(2 * u.deg, 2 * u.deg, semi_major, 0.4, 30 * u.deg) |
80
|
|
|
vals_1 = model_1(lon, lat) |
81
|
|
|
assert vals_1.unit == "sr-1" |
82
|
|
|
mymap_1 = Map.from_geom(m_geom_1, data=vals_1.value) |
83
|
|
|
assert_allclose( |
84
|
|
|
np.sum(mymap_1.quantity * u.sr ** -1 * m_geom_1.solid_angle()), 1, rtol=1.0e-3 |
85
|
|
|
) |
86
|
|
|
|
87
|
|
|
radius = model_1.evaluation_radius |
88
|
|
|
assert radius.unit == "deg" |
89
|
|
|
assert_allclose(radius.value, semi_major.value) |
90
|
|
|
|
91
|
|
|
# test the normalization for a disk (ellipse with e=0) at the Galactic Pole, |
92
|
|
|
# both analytically and comparing with the SkyDisk model |
93
|
|
|
m_geom_2 = WcsGeom.create( |
94
|
|
|
binsz=0.1, width=(6, 6), skydir=(0, 90), coordsys="GAL", proj="AIT" |
95
|
|
|
) |
96
|
|
|
coords = m_geom_2.get_coord() |
97
|
|
|
lon = coords.lon * u.deg |
98
|
|
|
lat = coords.lat * u.deg |
99
|
|
|
|
100
|
|
|
semi_major = 5 * u.deg |
101
|
|
|
model_2 = SkyEllipse(0 * u.deg, 90 * u.deg, semi_major, 0.0, 0.0 * u.deg) |
102
|
|
|
vals_2 = model_2(lon, lat) |
103
|
|
|
mymap_2 = Map.from_geom(m_geom_2, data=vals_2.value) |
104
|
|
|
|
105
|
|
|
disk = SkyDisk(lon_0="0 deg", lat_0="90 deg", r_0="5 deg") |
106
|
|
|
vals_disk = disk(lon, lat) |
107
|
|
|
mymap_disk = Map.from_geom(m_geom_2, data=vals_disk.value) |
108
|
|
|
|
109
|
|
|
solid_angle = 2 * np.pi * (1 - np.cos(5 * u.deg)) |
110
|
|
|
assert_allclose(np.max(vals_2).value * solid_angle, 1) |
111
|
|
|
|
112
|
|
|
assert_allclose( |
113
|
|
|
np.sum(mymap_2.quantity * u.sr ** -1 * m_geom_2.solid_angle()), |
114
|
|
|
np.sum(mymap_disk.quantity * u.sr ** -1 * m_geom_2.solid_angle()), |
115
|
|
|
) |
116
|
|
|
|
117
|
|
|
|
118
|
|
|
def test_sky_ellipse_edge(): |
119
|
|
|
pytest.importorskip("astropy", minversion="3.1.1") |
120
|
|
|
r_0 = 2 * u.deg |
121
|
|
|
model = SkyEllipse(lon_0="0 deg", lat_0="0 deg", semi_major=r_0, e=0.5, theta="0 deg") |
122
|
|
|
value_center = model(0 * u.deg, 0 * u.deg) |
123
|
|
|
value_edge = model(r_0, 0 * u.deg) |
124
|
|
|
assert_allclose((value_edge / value_center).to_value(""), 0.5) |
125
|
|
|
|
126
|
|
|
|
127
|
|
|
def test_sky_shell(): |
128
|
|
|
width = 2 * u.deg |
129
|
|
|
rad = 2 * u.deg |
130
|
|
|
model = SkyShell(lon_0="1 deg", lat_0="45 deg", radius=rad, width=width) |
131
|
|
|
lon = [1, 2, 4] * u.deg |
132
|
|
|
lat = 45 * u.deg |
133
|
|
|
val = model(lon, lat) |
134
|
|
|
assert val.unit == "deg-2" |
135
|
|
|
desired = [55.979449, 57.831651, 94.919895] |
136
|
|
|
assert_allclose(val.to_value("sr-1"), desired) |
137
|
|
|
radius = model.evaluation_radius |
138
|
|
|
assert radius.unit == "deg" |
139
|
|
|
assert_allclose(radius.value, rad.value + width.value) |
140
|
|
|
|
141
|
|
|
|
142
|
|
|
def test_sky_diffuse_constant(): |
143
|
|
|
model = SkyDiffuseConstant(value="42 sr-1") |
144
|
|
|
lon = [1, 2] * u.deg |
145
|
|
|
lat = 45 * u.deg |
146
|
|
|
val = model(lon, lat) |
147
|
|
|
assert val.unit == "sr-1" |
148
|
|
|
assert_allclose(val.value, 42) |
149
|
|
|
radius = model.evaluation_radius |
150
|
|
|
assert radius is None |
151
|
|
|
|
152
|
|
|
|
153
|
|
|
@requires_data("gammapy-data") |
154
|
|
|
def test_sky_diffuse_map(): |
155
|
|
|
filename = "$GAMMAPY_DATA/catalogs/fermi/Extended_archive_v18/Templates/RXJ1713_2016_250GeV.fits" |
156
|
|
|
model = SkyDiffuseMap.read(filename, normalize=False) |
157
|
|
|
lon = [258.5, 0] * u.deg |
158
|
|
|
lat = -39.8 * u.deg |
159
|
|
|
val = model(lon, lat) |
160
|
|
|
assert val.unit == "sr-1" |
161
|
|
|
desired = [3269.178107, 0] |
162
|
|
|
assert_allclose(val.value, desired) |
163
|
|
|
radius = model.evaluation_radius |
164
|
|
|
assert radius.unit == "deg" |
165
|
|
|
assert_allclose(radius.value, 0.64, rtol=1.0e-2) |
166
|
|
|
assert model.frame == "fk5" |
167
|
|
|
|
168
|
|
|
|
169
|
|
|
@requires_data("gammapy-data") |
170
|
|
|
def test_sky_diffuse_map_normalize(): |
171
|
|
|
# define model map with a constant value of 1 |
172
|
|
|
model_map = Map.create(map_type="wcs", width=(10, 5), binsz=0.5) |
173
|
|
|
model_map.data += 1.0 |
174
|
|
|
model = SkyDiffuseMap(model_map) |
175
|
|
|
|
176
|
|
|
# define data map with a different spatial binning |
177
|
|
|
data_map = Map.create(map_type="wcs", width=(10, 5), binsz=1) |
178
|
|
|
coords = data_map.geom.get_coord() |
179
|
|
|
solid_angle = data_map.geom.solid_angle() |
180
|
|
|
vals = model(coords.lon * u.deg, coords.lat * u.deg) * solid_angle |
181
|
|
|
|
182
|
|
|
assert vals.unit == "" |
183
|
|
|
integral = vals.sum() |
184
|
|
|
assert_allclose(integral.value, 1, rtol=1e-4) |
185
|
|
|
|