|
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
|
|
|
|
gammapy /
gammapy
