We are currently solving the following issue within our repository. Here is the issue text:
--- BEGIN ISSUE ---
Incorrect results for SkyPointSource
The SkyPointSource spatial model gives a wrong npred.

eg:

`spatial_model = SkyGaussian(
    lon_0='0.2 deg',
    lat_0='0.1 deg',
    sigma='0.01 deg',
)
spectral_model = PowerLaw(
    index=3,
    amplitude='1e-11 cm-2 s-1 TeV-1',
    reference='1 TeV',
)
sky_model = SkyModel(
    spatial_model=spatial_model,
    spectral_model=spectral_model,
)'
'evaluator = MapEvaluator(sky_model, exposure=exposure_map)'
'npred = evaluator.compute_npred()'
'print(npred.sum())`

gives  `3846` predicted events in total.

`spatial_model1 = SkyPointSource(
    lon_0='0.2 deg',
    lat_0='0.1 deg',
)
spectral_model = PowerLaw(
    index=3,
    amplitude='1e-11 cm-2 s-1 TeV-1',
    reference='1 TeV',
)
sky_model1 = SkyModel(
    spatial_model=spatial_model1,
    spectral_model=spectral_model,
)'
`evaluator = MapEvaluator(sky_model1, exposure=exposure_map)`
`npred = evaluator.compute_npred()`
`print(npred.sum())`

gives `0.0004824` predicted events in total



(for the detailed code, see https://gist.github.com/AtreyeeS/d5161ebe30f29c1a19b34757a1916301)





--- END ISSUE ---

Below are some code segments, each from a relevant file. One or more of these files may contain bugs.
--- BEGIN FILES ---
Path: `gammapy/image/models/new.py`
Content:
```
1 # Licensed under a 3-clause BSD style license - see LICENSE.rst
2 """
3 Spatial models for astrophysical gamma-ray sources.
4 """
5 from __future__ import absolute_import, division, print_function, unicode_literals
6 import copy
7 import abc
8 import numpy as np
9 import astropy.units as u
10 from astropy.coordinates.angle_utilities import angular_separation
11 from astropy.coordinates import Angle, Longitude, Latitude
12 from ...extern import six
13 from ...utils.modeling import Parameter, ParameterList
14 from ...maps import Map
15 
16 __all__ = [
17     'SkySpatialModel',
18     'SkyPointSource',
19     'SkyGaussian',
20     'SkyDisk',
21     'SkyShell',
22     'SkyDiffuseConstant',
23     'SkyDiffuseMap',
24 ]
25 
26 
27 @six.add_metaclass(abc.ABCMeta)
28 class SkySpatialModel(object):
29     """SkySpatial model base class.
30     """
31 
32     def __str__(self):
33         ss = self.__class__.__name__
34         ss += '\n\nParameters: \n\n\t'
35 
36         table = self.parameters.to_table()
37         ss += '\n\t'.join(table.pformat())
38 
39         if self.parameters.covariance is not None:
40             ss += '\n\nCovariance: \n\n\t'
41             covar = self.parameters.covariance_to_table()
42             ss += '\n\t'.join(covar.pformat())
43         return ss
44 
45     def __call__(self, lon, lat):
46         """Call evaluate method"""
47         kwargs = dict()
48         for par in self.parameters.parameters:
49             kwargs[par.name] = par.quantity
50 
51         return self.evaluate(lon, lat, **kwargs)
52 
53     def copy(self):
54         """A deep copy."""
55         return copy.deepcopy(self)
56 
57 
58 class SkyPointSource(SkySpatialModel):
59     r"""Point Source.
60 
61     .. math::
62 
63         \phi(lon, lat) = \delta{(lon - lon_0, lat - lat_0)}
64 
65     A tolerance of 1 arcsecond is accepted for numerical stability
66 
67     Parameters
68     ----------
69     lon_0 : `~astropy.coordinates.Longitude`
70         :math:`lon_0`
71     lat_0 : `~astropy.coordinates.Latitude`
72         :math:`lat_0`
73     """
74 
75     def __init__(self, lon_0, lat_0):
76         self.parameters = ParameterList([
77             Parameter('lon_0', Longitude(lon_0)),
78             Parameter('lat_0', Latitude(lat_0))
79         ])
80 
81     @staticmethod
82     def evaluate(lon, lat, lon_0, lat_0):
83         """Evaluate the model (static function)."""
84 
85         wrapval=lon_0 + 180*u.deg
86         lon = Angle(lon).wrap_at(wrapval)
87 
88         _, grad_lon = np.gradient(lon)
89         grad_lat, _ = np.gradient(lat)
90         lon_diff = np.abs((lon - lon_0) / grad_lon)
91         lat_diff = np.abs((lat - lat_0) / grad_lat)
92 
93         lon_val = np.select([lon_diff < 1], [1 - lon_diff], 0)
94         lat_val = np.select([lat_diff < 1], [1 - lat_diff], 0)
95         val = lon_val * lat_val
96         return val * u.Unit('sr-1')
97 
98 
99 class SkyGaussian(SkySpatialModel):
100     r"""Two-dimensional symmetric Gaussian model.
101 
102     .. math::
103 
104         \phi(lon, lat) = \frac{1}{2\pi\sigma^2} \exp{\left(-\frac{1}{2}
105             \frac{\theta^2}{\sigma^2}\right)}
106 
107     where :math:`\theta` is the sky separation
108 
109     Parameters
110     ----------
111     lon_0 : `~astropy.coordinates.Longitude`
112         :math:`lon_0`
113     lat_0 : `~astropy.coordinates.Latitude`
114         :math:`lat_0`
115     sigma : `~astropy.coordinates.Angle`
116         :math:`\sigma`
117     """
118 
119     def __init__(self, lon_0, lat_0, sigma):
120         self.parameters = ParameterList([
121             Parameter('lon_0', Longitude(lon_0)),
122             Parameter('lat_0', Latitude(lat_0)),
123             Parameter('sigma', Angle(sigma))
124         ])
125 
126     @staticmethod
127     def evaluate(lon, lat, lon_0, lat_0, sigma):
128         """Evaluate the model (static function)."""
129         sep = angular_separation(lon, lat, lon_0, lat_0)
130         sep = sep.to('rad').value
131         sigma = sigma.to('rad').value
132 
133         norm = 1 / (2 * np.pi * sigma ** 2)
134         exponent = -0.5 * (sep / sigma) ** 2
135         val = norm * np.exp(exponent)
136 
137         return val * u.Unit('sr-1')
138 
139 
140 class SkyDisk(SkySpatialModel):
141     r"""Constant radial disk model.
142 
143     .. math::
144 
145         \phi(lon, lat) = \frac{1}{2 \pi (1 - \cos{r}) } \cdot
146                 \begin{cases}
147                     1 & \text{for } \theta \leq r_0 \\
148                     0 & \text{for } \theta < r_0
149                 \end{cases}
150 
151     where :math:`\theta` is the sky separation
152 
153     Parameters
154     ----------
155     lon_0 : `~astropy.coordinates.Longitude`
156         :math:`lon_0`
157     lat_0 : `~astropy.coordinates.Latitude`
158         :math:`lat_0`
159     r_0 : `~astropy.coordinates.Angle`
160         :math:`r_0`
161     """
162 
163     def __init__(self, lon_0, lat_0, r_0):
164         self.parameters = ParameterList([
165             Parameter('lon_0', Longitude(lon_0)),
166             Parameter('lat_0', Latitude(lat_0)),
167             Parameter('r_0', Angle(r_0))
168         ])
169 
170     @staticmethod
171     def evaluate(lon, lat, lon_0, lat_0, r_0):
172         """Evaluate the model (static function)."""
173         sep = angular_separation(lon, lat, lon_0, lat_0)
174         sep = sep.to('rad').value
175         r_0 = r_0.to('rad').value
176 
177         norm = 1. / (2 * np.pi * (1 - np.cos(r_0)))
178         val = np.where(sep <= r_0, norm, 0)
179 
180         return val * u.Unit('sr-1')
181 
182 
183 class SkyShell(SkySpatialModel):
184     r"""Shell model
185 
186     .. math::
187 
188         \phi(lon, lat) = \frac{3}{2 \pi (r_{out}^3 - r_{in}^3)} \cdot
189                 \begin{cases}
190                     \sqrt{r_{out}^2 - \theta^2} - \sqrt{r_{in}^2 - \theta^2} &
191                                  \text{for } \theta \lt r_{in} \\
192                     \sqrt{r_{out}^2 - \theta^2} & 
193                                  \text{for } r_{in} \leq \theta \lt r_{out} \\
194                     0 & \text{for } \theta > r_{out}
195                 \end{cases}
196 
197     where :math:`\theta` is the sky separation and :math:`r_out = r_in` + width
198 
199     Note that the normalization is a small angle approximation,
200     although that approximation is still very good even for 10 deg radius shells.
201 
202     Parameters
203     ----------
204     lon_0 : `~astropy.coordinates.Longitude`
205         :math:`lon_0`
206     lat_0 : `~astropy.coordinates.Latitude`
207         :math:`lat_0`
208     radius : `~astropy.coordinates.Angle`
209         Inner radius, :math:`r_{in}`
210     width : `~astropy.coordinates.Angle`
211         Shell width
212     """
213 
214     def __init__(self, lon_0, lat_0, radius, width):
215         self.parameters = ParameterList([
216             Parameter('lon_0', Longitude(lon_0)),
217             Parameter('lat_0', Latitude(lat_0)),
218             Parameter('radius', Angle(radius)),
219             Parameter('width', Angle(width))
220         ])
221 
222     @staticmethod
223     def evaluate(lon, lat, lon_0, lat_0, radius, width):
224         """Evaluate the model (static function)."""
225         sep = angular_separation(lon, lat, lon_0, lat_0)
226         sep = sep.to('rad').value
227         r_i = radius.to('rad').value
228         r_o = (radius + width).to('rad').value
229 
230         norm = 3 / (2 * np.pi * (r_o ** 3 - r_i ** 3))
231 
232         with np.errstate(invalid='ignore'):
233             val_out = np.sqrt(r_o ** 2 - sep ** 2)
234             val_in = val_out - np.sqrt(r_i ** 2 - sep ** 2)
235             val = np.select([sep < r_i, sep < r_o], [val_in, val_out])
236 
237         return norm * val * u.Unit('sr-1')
238 
239 
240 class SkyDiffuseConstant(SkySpatialModel):
241     """Spatially constant (isotropic) spatial model.
242 
243     Parameters
244     ----------
245     value : `~astropy.units.Quantity`
246         Value
247     """
248 
249     def __init__(self, value=1):
250         self.parameters = ParameterList([
251             Parameter('value', value),
252         ])
253 
254     @staticmethod
255     def evaluate(lon, lat, value):
256         # TODO: try fitting this -> probably the interface doesn't work?!
257         return value
258 
259 
260 class SkyDiffuseMap(SkySpatialModel):
261     """Spatial sky map template model.
262 
263     At the moment only support 2D maps.
264     TODO: support maps with an energy axis here or in a separate class?
265     TODO: should we cache some interpolator object for efficiency?
266 
267     Parameters
268     ----------
269     map : `~gammapy.map.Map`
270         Map template
271     norm : `~astropy.units.Quantity`
272         Norm parameter (multiplied with map values)
273     meta : dict, optional
274         Meta information, meta['filename'] will be used for serialization
275     """
276 
277     def __init__(self, map, norm=1, meta=None):
278         self._map = map
279         self.parameters = ParameterList([
280             Parameter('norm', norm),
281         ])
282         self.meta = dict() if meta is None else meta
283 
284     @classmethod
285     def read(cls, filename, **kwargs):
286         """Read spatial template model from FITS image.
287 
288         Parameters
289         ----------
290         filename : str
291             FITS image filename.
292         """
293         template = Map.read(filename, **kwargs)
294         return cls(template)
295 
296     def evaluate(self, lon, lat, norm):
297         coord = dict(
298             lon=lon.to('deg').value,
299             lat=lat.to('deg').value,
300         )
301         val = self._map.interp_by_coord(coord, fill_value=0)
302         # TODO: use map unit? self._map.unit
303         return norm * val * u.Unit('sr-1')
304
```
--- END FILES ---

Please first localize the bug based on the issue statement, and then generate a patch according to the `git diff` format fenced by three backticks.

Here is an example:

```diff
diff --git a/examples/server_async.py b/examples/server_async.py
--- a/examples/server_async.py
+++ b/examples/server_async.py
@@ -313,4 +313,4 @@
 
 
 if __name__ == "__main__":
-    asyncio.run(run_async_server("."), debug=True)
+    asyncio.run(run_async_server(), debug=True)
diff --git a/examples/server_sync.py b/examples/server_sync.py
--- a/examples/server_sync.py
+++ b/examples/server_sync.py
@@ -313,5 +313,5 @@
 
 
 if __name__ == "__main__":
-    server = run_sync_server(".")
+    server = run_sync_server()
     server.shutdown()

```

diff --git a/gammapy/image/models/new.py b/gammapy/image/models/new.py
--- a/gammapy/image/models/new.py
+++ b/gammapy/image/models/new.py
@@ -82,7 +82,7 @@
     def evaluate(lon, lat, lon_0, lat_0):
         """Evaluate the model (static function)."""
 
-        wrapval=lon_0 + 180*u.deg
+        wrapval = lon_0 + 180 * u.deg
         lon = Angle(lon).wrap_at(wrapval)
 
         _, grad_lon = np.gradient(lon)
@@ -90,10 +90,10 @@
         lon_diff = np.abs((lon - lon_0) / grad_lon)
         lat_diff = np.abs((lat - lat_0) / grad_lat)
 
-        lon_val = np.select([lon_diff < 1], [1 - lon_diff], 0)
-        lat_val = np.select([lat_diff < 1], [1 - lat_diff], 0)
+        lon_val = np.select([lon_diff < 1], [1 - lon_diff], 0) / np.abs(grad_lon)
+        lat_val = np.select([lat_diff < 1], [1 - lat_diff], 0) / np.abs(grad_lat)
         val = lon_val * lat_val
-        return val * u.Unit('sr-1')
+        return val.to('sr-1')
 
 
 class SkyGaussian(SkySpatialModel):
@@ -189,7 +189,7 @@
                 \begin{cases}
                     \sqrt{r_{out}^2 - \theta^2} - \sqrt{r_{in}^2 - \theta^2} &
                                  \text{for } \theta \lt r_{in} \\
-                    \sqrt{r_{out}^2 - \theta^2} & 
+                    \sqrt{r_{out}^2 - \theta^2} &
                                  \text{for } r_{in} \leq \theta \lt r_{out} \\
                     0 & \text{for } \theta > r_{out}
                 \end{cases}

{"golden_diff": "diff --git a/gammapy/image/models/new.py b/gammapy/image/models/new.py\n--- a/gammapy/image/models/new.py\n+++ b/gammapy/image/models/new.py\n@@ -82,7 +82,7 @@\n     def evaluate(lon, lat, lon_0, lat_0):\n         \"\"\"Evaluate the model (static function).\"\"\"\n \n-        wrapval=lon_0 + 180*u.deg\n+        wrapval = lon_0 + 180 * u.deg\n         lon = Angle(lon).wrap_at(wrapval)\n \n         _, grad_lon = np.gradient(lon)\n@@ -90,10 +90,10 @@\n         lon_diff = np.abs((lon - lon_0) / grad_lon)\n         lat_diff = np.abs((lat - lat_0) / grad_lat)\n \n-        lon_val = np.select([lon_diff < 1], [1 - lon_diff], 0)\n-        lat_val = np.select([lat_diff < 1], [1 - lat_diff], 0)\n+        lon_val = np.select([lon_diff < 1], [1 - lon_diff], 0) / np.abs(grad_lon)\n+        lat_val = np.select([lat_diff < 1], [1 - lat_diff], 0) / np.abs(grad_lat)\n         val = lon_val * lat_val\n-        return val * u.Unit('sr-1')\n+        return val.to('sr-1')\n \n \n class SkyGaussian(SkySpatialModel):\n@@ -189,7 +189,7 @@\n                 \\begin{cases}\n                     \\sqrt{r_{out}^2 - \\theta^2} - \\sqrt{r_{in}^2 - \\theta^2} &\n                                  \\text{for } \\theta \\lt r_{in} \\\\\n-                    \\sqrt{r_{out}^2 - \\theta^2} & \n+                    \\sqrt{r_{out}^2 - \\theta^2} &\n                                  \\text{for } r_{in} \\leq \\theta \\lt r_{out} \\\\\n                     0 & \\text{for } \\theta > r_{out}\n                 \\end{cases}\n", "issue": "Incorrect results for SkyPointSource\nThe SkyPointSource spatial model gives a wrong npred.\r\n\r\neg:\r\n\r\n`spatial_model = SkyGaussian(\r\n    lon_0='0.2 deg',\r\n    lat_0='0.1 deg',\r\n    sigma='0.01 deg',\r\n)\r\nspectral_model = PowerLaw(\r\n    index=3,\r\n    amplitude='1e-11 cm-2 s-1 TeV-1',\r\n    reference='1 TeV',\r\n)\r\nsky_model = SkyModel(\r\n    spatial_model=spatial_model,\r\n    spectral_model=spectral_model,\r\n)'\r\n'evaluator = MapEvaluator(sky_model, exposure=exposure_map)'\r\n'npred = evaluator.compute_npred()'\r\n'print(npred.sum())`\r\n\r\ngives  `3846` predicted events in total.\r\n\r\n`spatial_model1 = SkyPointSource(\r\n    lon_0='0.2 deg',\r\n    lat_0='0.1 deg',\r\n)\r\nspectral_model = PowerLaw(\r\n    index=3,\r\n    amplitude='1e-11 cm-2 s-1 TeV-1',\r\n    reference='1 TeV',\r\n)\r\nsky_model1 = SkyModel(\r\n    spatial_model=spatial_model1,\r\n    spectral_model=spectral_model,\r\n)'\r\n`evaluator = MapEvaluator(sky_model1, exposure=exposure_map)`\r\n`npred = evaluator.compute_npred()`\r\n`print(npred.sum())`\r\n\r\ngives `0.0004824` predicted events in total\r\n\r\n\r\n\r\n(for the detailed code, see https://gist.github.com/AtreyeeS/d5161ebe30f29c1a19b34757a1916301)\r\n\r\n\r\n\r\n\n", "before_files": [{"content": "# Licensed under a 3-clause BSD style license - see LICENSE.rst\n\"\"\"\nSpatial models for astrophysical gamma-ray sources.\n\"\"\"\nfrom __future__ import absolute_import, division, print_function, unicode_literals\nimport copy\nimport abc\nimport numpy as np\nimport astropy.units as u\nfrom astropy.coordinates.angle_utilities import angular_separation\nfrom astropy.coordinates import Angle, Longitude, Latitude\nfrom ...extern import six\nfrom ...utils.modeling import Parameter, ParameterList\nfrom ...maps import Map\n\n__all__ = [\n    'SkySpatialModel',\n    'SkyPointSource',\n    'SkyGaussian',\n    'SkyDisk',\n    'SkyShell',\n    'SkyDiffuseConstant',\n    'SkyDiffuseMap',\n]\n\n\n@six.add_metaclass(abc.ABCMeta)\nclass SkySpatialModel(object):\n    \"\"\"SkySpatial model base class.\n    \"\"\"\n\n    def __str__(self):\n        ss = self.__class__.__name__\n        ss += '\\n\\nParameters: \\n\\n\\t'\n\n        table = self.parameters.to_table()\n        ss += '\\n\\t'.join(table.pformat())\n\n        if self.parameters.covariance is not None:\n            ss += '\\n\\nCovariance: \\n\\n\\t'\n            covar = self.parameters.covariance_to_table()\n            ss += '\\n\\t'.join(covar.pformat())\n        return ss\n\n    def __call__(self, lon, lat):\n        \"\"\"Call evaluate method\"\"\"\n        kwargs = dict()\n        for par in self.parameters.parameters:\n            kwargs[par.name] = par.quantity\n\n        return self.evaluate(lon, lat, **kwargs)\n\n    def copy(self):\n        \"\"\"A deep copy.\"\"\"\n        return copy.deepcopy(self)\n\n\nclass SkyPointSource(SkySpatialModel):\n    r\"\"\"Point Source.\n\n    .. math::\n\n        \\phi(lon, lat) = \\delta{(lon - lon_0, lat - lat_0)}\n\n    A tolerance of 1 arcsecond is accepted for numerical stability\n\n    Parameters\n    ----------\n    lon_0 : `~astropy.coordinates.Longitude`\n        :math:`lon_0`\n    lat_0 : `~astropy.coordinates.Latitude`\n        :math:`lat_0`\n    \"\"\"\n\n    def __init__(self, lon_0, lat_0):\n        self.parameters = ParameterList([\n            Parameter('lon_0', Longitude(lon_0)),\n            Parameter('lat_0', Latitude(lat_0))\n        ])\n\n    @staticmethod\n    def evaluate(lon, lat, lon_0, lat_0):\n        \"\"\"Evaluate the model (static function).\"\"\"\n\n        wrapval=lon_0 + 180*u.deg\n        lon = Angle(lon).wrap_at(wrapval)\n\n        _, grad_lon = np.gradient(lon)\n        grad_lat, _ = np.gradient(lat)\n        lon_diff = np.abs((lon - lon_0) / grad_lon)\n        lat_diff = np.abs((lat - lat_0) / grad_lat)\n\n        lon_val = np.select([lon_diff < 1], [1 - lon_diff], 0)\n        lat_val = np.select([lat_diff < 1], [1 - lat_diff], 0)\n        val = lon_val * lat_val\n        return val * u.Unit('sr-1')\n\n\nclass SkyGaussian(SkySpatialModel):\n    r\"\"\"Two-dimensional symmetric Gaussian model.\n\n    .. math::\n\n        \\phi(lon, lat) = \\frac{1}{2\\pi\\sigma^2} \\exp{\\left(-\\frac{1}{2}\n            \\frac{\\theta^2}{\\sigma^2}\\right)}\n\n    where :math:`\\theta` is the sky separation\n\n    Parameters\n    ----------\n    lon_0 : `~astropy.coordinates.Longitude`\n        :math:`lon_0`\n    lat_0 : `~astropy.coordinates.Latitude`\n        :math:`lat_0`\n    sigma : `~astropy.coordinates.Angle`\n        :math:`\\sigma`\n    \"\"\"\n\n    def __init__(self, lon_0, lat_0, sigma):\n        self.parameters = ParameterList([\n            Parameter('lon_0', Longitude(lon_0)),\n            Parameter('lat_0', Latitude(lat_0)),\n            Parameter('sigma', Angle(sigma))\n        ])\n\n    @staticmethod\n    def evaluate(lon, lat, lon_0, lat_0, sigma):\n        \"\"\"Evaluate the model (static function).\"\"\"\n        sep = angular_separation(lon, lat, lon_0, lat_0)\n        sep = sep.to('rad').value\n        sigma = sigma.to('rad').value\n\n        norm = 1 / (2 * np.pi * sigma ** 2)\n        exponent = -0.5 * (sep / sigma) ** 2\n        val = norm * np.exp(exponent)\n\n        return val * u.Unit('sr-1')\n\n\nclass SkyDisk(SkySpatialModel):\n    r\"\"\"Constant radial disk model.\n\n    .. math::\n\n        \\phi(lon, lat) = \\frac{1}{2 \\pi (1 - \\cos{r}) } \\cdot\n                \\begin{cases}\n                    1 & \\text{for } \\theta \\leq r_0 \\\\\n                    0 & \\text{for } \\theta < r_0\n                \\end{cases}\n\n    where :math:`\\theta` is the sky separation\n\n    Parameters\n    ----------\n    lon_0 : `~astropy.coordinates.Longitude`\n        :math:`lon_0`\n    lat_0 : `~astropy.coordinates.Latitude`\n        :math:`lat_0`\n    r_0 : `~astropy.coordinates.Angle`\n        :math:`r_0`\n    \"\"\"\n\n    def __init__(self, lon_0, lat_0, r_0):\n        self.parameters = ParameterList([\n            Parameter('lon_0', Longitude(lon_0)),\n            Parameter('lat_0', Latitude(lat_0)),\n            Parameter('r_0', Angle(r_0))\n        ])\n\n    @staticmethod\n    def evaluate(lon, lat, lon_0, lat_0, r_0):\n        \"\"\"Evaluate the model (static function).\"\"\"\n        sep = angular_separation(lon, lat, lon_0, lat_0)\n        sep = sep.to('rad').value\n        r_0 = r_0.to('rad').value\n\n        norm = 1. / (2 * np.pi * (1 - np.cos(r_0)))\n        val = np.where(sep <= r_0, norm, 0)\n\n        return val * u.Unit('sr-1')\n\n\nclass SkyShell(SkySpatialModel):\n    r\"\"\"Shell model\n\n    .. math::\n\n        \\phi(lon, lat) = \\frac{3}{2 \\pi (r_{out}^3 - r_{in}^3)} \\cdot\n                \\begin{cases}\n                    \\sqrt{r_{out}^2 - \\theta^2} - \\sqrt{r_{in}^2 - \\theta^2} &\n                                 \\text{for } \\theta \\lt r_{in} \\\\\n                    \\sqrt{r_{out}^2 - \\theta^2} & \n                                 \\text{for } r_{in} \\leq \\theta \\lt r_{out} \\\\\n                    0 & \\text{for } \\theta > r_{out}\n                \\end{cases}\n\n    where :math:`\\theta` is the sky separation and :math:`r_out = r_in` + width\n\n    Note that the normalization is a small angle approximation,\n    although that approximation is still very good even for 10 deg radius shells.\n\n    Parameters\n    ----------\n    lon_0 : `~astropy.coordinates.Longitude`\n        :math:`lon_0`\n    lat_0 : `~astropy.coordinates.Latitude`\n        :math:`lat_0`\n    radius : `~astropy.coordinates.Angle`\n        Inner radius, :math:`r_{in}`\n    width : `~astropy.coordinates.Angle`\n        Shell width\n    \"\"\"\n\n    def __init__(self, lon_0, lat_0, radius, width):\n        self.parameters = ParameterList([\n            Parameter('lon_0', Longitude(lon_0)),\n            Parameter('lat_0', Latitude(lat_0)),\n            Parameter('radius', Angle(radius)),\n            Parameter('width', Angle(width))\n        ])\n\n    @staticmethod\n    def evaluate(lon, lat, lon_0, lat_0, radius, width):\n        \"\"\"Evaluate the model (static function).\"\"\"\n        sep = angular_separation(lon, lat, lon_0, lat_0)\n        sep = sep.to('rad').value\n        r_i = radius.to('rad').value\n        r_o = (radius + width).to('rad').value\n\n        norm = 3 / (2 * np.pi * (r_o ** 3 - r_i ** 3))\n\n        with np.errstate(invalid='ignore'):\n            val_out = np.sqrt(r_o ** 2 - sep ** 2)\n            val_in = val_out - np.sqrt(r_i ** 2 - sep ** 2)\n            val = np.select([sep < r_i, sep < r_o], [val_in, val_out])\n\n        return norm * val * u.Unit('sr-1')\n\n\nclass SkyDiffuseConstant(SkySpatialModel):\n    \"\"\"Spatially constant (isotropic) spatial model.\n\n    Parameters\n    ----------\n    value : `~astropy.units.Quantity`\n        Value\n    \"\"\"\n\n    def __init__(self, value=1):\n        self.parameters = ParameterList([\n            Parameter('value', value),\n        ])\n\n    @staticmethod\n    def evaluate(lon, lat, value):\n        # TODO: try fitting this -> probably the interface doesn't work?!\n        return value\n\n\nclass SkyDiffuseMap(SkySpatialModel):\n    \"\"\"Spatial sky map template model.\n\n    At the moment only support 2D maps.\n    TODO: support maps with an energy axis here or in a separate class?\n    TODO: should we cache some interpolator object for efficiency?\n\n    Parameters\n    ----------\n    map : `~gammapy.map.Map`\n        Map template\n    norm : `~astropy.units.Quantity`\n        Norm parameter (multiplied with map values)\n    meta : dict, optional\n        Meta information, meta['filename'] will be used for serialization\n    \"\"\"\n\n    def __init__(self, map, norm=1, meta=None):\n        self._map = map\n        self.parameters = ParameterList([\n            Parameter('norm', norm),\n        ])\n        self.meta = dict() if meta is None else meta\n\n    @classmethod\n    def read(cls, filename, **kwargs):\n        \"\"\"Read spatial template model from FITS image.\n\n        Parameters\n        ----------\n        filename : str\n            FITS image filename.\n        \"\"\"\n        template = Map.read(filename, **kwargs)\n        return cls(template)\n\n    def evaluate(self, lon, lat, norm):\n        coord = dict(\n            lon=lon.to('deg').value,\n            lat=lat.to('deg').value,\n        )\n        val = self._map.interp_by_coord(coord, fill_value=0)\n        # TODO: use map unit? self._map.unit\n        return norm * val * u.Unit('sr-1')\n", "path": "gammapy/image/models/new.py"}], "after_files": [{"content": "# Licensed under a 3-clause BSD style license - see LICENSE.rst\n\"\"\"\nSpatial models for astrophysical gamma-ray sources.\n\"\"\"\nfrom __future__ import absolute_import, division, print_function, unicode_literals\nimport copy\nimport abc\nimport numpy as np\nimport astropy.units as u\nfrom astropy.coordinates.angle_utilities import angular_separation\nfrom astropy.coordinates import Angle, Longitude, Latitude\nfrom ...extern import six\nfrom ...utils.modeling import Parameter, ParameterList\nfrom ...maps import Map\n\n__all__ = [\n    'SkySpatialModel',\n    'SkyPointSource',\n    'SkyGaussian',\n    'SkyDisk',\n    'SkyShell',\n    'SkyDiffuseConstant',\n    'SkyDiffuseMap',\n]\n\n\n@six.add_metaclass(abc.ABCMeta)\nclass SkySpatialModel(object):\n    \"\"\"SkySpatial model base class.\n    \"\"\"\n\n    def __str__(self):\n        ss = self.__class__.__name__\n        ss += '\\n\\nParameters: \\n\\n\\t'\n\n        table = self.parameters.to_table()\n        ss += '\\n\\t'.join(table.pformat())\n\n        if self.parameters.covariance is not None:\n            ss += '\\n\\nCovariance: \\n\\n\\t'\n            covar = self.parameters.covariance_to_table()\n            ss += '\\n\\t'.join(covar.pformat())\n        return ss\n\n    def __call__(self, lon, lat):\n        \"\"\"Call evaluate method\"\"\"\n        kwargs = dict()\n        for par in self.parameters.parameters:\n            kwargs[par.name] = par.quantity\n\n        return self.evaluate(lon, lat, **kwargs)\n\n    def copy(self):\n        \"\"\"A deep copy.\"\"\"\n        return copy.deepcopy(self)\n\n\nclass SkyPointSource(SkySpatialModel):\n    r\"\"\"Point Source.\n\n    .. math::\n\n        \\phi(lon, lat) = \\delta{(lon - lon_0, lat - lat_0)}\n\n    A tolerance of 1 arcsecond is accepted for numerical stability\n\n    Parameters\n    ----------\n    lon_0 : `~astropy.coordinates.Longitude`\n        :math:`lon_0`\n    lat_0 : `~astropy.coordinates.Latitude`\n        :math:`lat_0`\n    \"\"\"\n\n    def __init__(self, lon_0, lat_0):\n        self.parameters = ParameterList([\n            Parameter('lon_0', Longitude(lon_0)),\n            Parameter('lat_0', Latitude(lat_0))\n        ])\n\n    @staticmethod\n    def evaluate(lon, lat, lon_0, lat_0):\n        \"\"\"Evaluate the model (static function).\"\"\"\n\n        wrapval = lon_0 + 180 * u.deg\n        lon = Angle(lon).wrap_at(wrapval)\n\n        _, grad_lon = np.gradient(lon)\n        grad_lat, _ = np.gradient(lat)\n        lon_diff = np.abs((lon - lon_0) / grad_lon)\n        lat_diff = np.abs((lat - lat_0) / grad_lat)\n\n        lon_val = np.select([lon_diff < 1], [1 - lon_diff], 0) / np.abs(grad_lon)\n        lat_val = np.select([lat_diff < 1], [1 - lat_diff], 0) / np.abs(grad_lat)\n        val = lon_val * lat_val\n        return val.to('sr-1')\n\n\nclass SkyGaussian(SkySpatialModel):\n    r\"\"\"Two-dimensional symmetric Gaussian model.\n\n    .. math::\n\n        \\phi(lon, lat) = \\frac{1}{2\\pi\\sigma^2} \\exp{\\left(-\\frac{1}{2}\n            \\frac{\\theta^2}{\\sigma^2}\\right)}\n\n    where :math:`\\theta` is the sky separation\n\n    Parameters\n    ----------\n    lon_0 : `~astropy.coordinates.Longitude`\n        :math:`lon_0`\n    lat_0 : `~astropy.coordinates.Latitude`\n        :math:`lat_0`\n    sigma : `~astropy.coordinates.Angle`\n        :math:`\\sigma`\n    \"\"\"\n\n    def __init__(self, lon_0, lat_0, sigma):\n        self.parameters = ParameterList([\n            Parameter('lon_0', Longitude(lon_0)),\n            Parameter('lat_0', Latitude(lat_0)),\n            Parameter('sigma', Angle(sigma))\n        ])\n\n    @staticmethod\n    def evaluate(lon, lat, lon_0, lat_0, sigma):\n        \"\"\"Evaluate the model (static function).\"\"\"\n        sep = angular_separation(lon, lat, lon_0, lat_0)\n        sep = sep.to('rad').value\n        sigma = sigma.to('rad').value\n\n        norm = 1 / (2 * np.pi * sigma ** 2)\n        exponent = -0.5 * (sep / sigma) ** 2\n        val = norm * np.exp(exponent)\n\n        return val * u.Unit('sr-1')\n\n\nclass SkyDisk(SkySpatialModel):\n    r\"\"\"Constant radial disk model.\n\n    .. math::\n\n        \\phi(lon, lat) = \\frac{1}{2 \\pi (1 - \\cos{r}) } \\cdot\n                \\begin{cases}\n                    1 & \\text{for } \\theta \\leq r_0 \\\\\n                    0 & \\text{for } \\theta < r_0\n                \\end{cases}\n\n    where :math:`\\theta` is the sky separation\n\n    Parameters\n    ----------\n    lon_0 : `~astropy.coordinates.Longitude`\n        :math:`lon_0`\n    lat_0 : `~astropy.coordinates.Latitude`\n        :math:`lat_0`\n    r_0 : `~astropy.coordinates.Angle`\n        :math:`r_0`\n    \"\"\"\n\n    def __init__(self, lon_0, lat_0, r_0):\n        self.parameters = ParameterList([\n            Parameter('lon_0', Longitude(lon_0)),\n            Parameter('lat_0', Latitude(lat_0)),\n            Parameter('r_0', Angle(r_0))\n        ])\n\n    @staticmethod\n    def evaluate(lon, lat, lon_0, lat_0, r_0):\n        \"\"\"Evaluate the model (static function).\"\"\"\n        sep = angular_separation(lon, lat, lon_0, lat_0)\n        sep = sep.to('rad').value\n        r_0 = r_0.to('rad').value\n\n        norm = 1. / (2 * np.pi * (1 - np.cos(r_0)))\n        val = np.where(sep <= r_0, norm, 0)\n\n        return val * u.Unit('sr-1')\n\n\nclass SkyShell(SkySpatialModel):\n    r\"\"\"Shell model\n\n    .. math::\n\n        \\phi(lon, lat) = \\frac{3}{2 \\pi (r_{out}^3 - r_{in}^3)} \\cdot\n                \\begin{cases}\n                    \\sqrt{r_{out}^2 - \\theta^2} - \\sqrt{r_{in}^2 - \\theta^2} &\n                                 \\text{for } \\theta \\lt r_{in} \\\\\n                    \\sqrt{r_{out}^2 - \\theta^2} &\n                                 \\text{for } r_{in} \\leq \\theta \\lt r_{out} \\\\\n                    0 & \\text{for } \\theta > r_{out}\n                \\end{cases}\n\n    where :math:`\\theta` is the sky separation and :math:`r_out = r_in` + width\n\n    Note that the normalization is a small angle approximation,\n    although that approximation is still very good even for 10 deg radius shells.\n\n    Parameters\n    ----------\n    lon_0 : `~astropy.coordinates.Longitude`\n        :math:`lon_0`\n    lat_0 : `~astropy.coordinates.Latitude`\n        :math:`lat_0`\n    radius : `~astropy.coordinates.Angle`\n        Inner radius, :math:`r_{in}`\n    width : `~astropy.coordinates.Angle`\n        Shell width\n    \"\"\"\n\n    def __init__(self, lon_0, lat_0, radius, width):\n        self.parameters = ParameterList([\n            Parameter('lon_0', Longitude(lon_0)),\n            Parameter('lat_0', Latitude(lat_0)),\n            Parameter('radius', Angle(radius)),\n            Parameter('width', Angle(width))\n        ])\n\n    @staticmethod\n    def evaluate(lon, lat, lon_0, lat_0, radius, width):\n        \"\"\"Evaluate the model (static function).\"\"\"\n        sep = angular_separation(lon, lat, lon_0, lat_0)\n        sep = sep.to('rad').value\n        r_i = radius.to('rad').value\n        r_o = (radius + width).to('rad').value\n\n        norm = 3 / (2 * np.pi * (r_o ** 3 - r_i ** 3))\n\n        with np.errstate(invalid='ignore'):\n            val_out = np.sqrt(r_o ** 2 - sep ** 2)\n            val_in = val_out - np.sqrt(r_i ** 2 - sep ** 2)\n            val = np.select([sep < r_i, sep < r_o], [val_in, val_out])\n\n        return norm * val * u.Unit('sr-1')\n\n\nclass SkyDiffuseConstant(SkySpatialModel):\n    \"\"\"Spatially constant (isotropic) spatial model.\n\n    Parameters\n    ----------\n    value : `~astropy.units.Quantity`\n        Value\n    \"\"\"\n\n    def __init__(self, value=1):\n        self.parameters = ParameterList([\n            Parameter('value', value),\n        ])\n\n    @staticmethod\n    def evaluate(lon, lat, value):\n        # TODO: try fitting this -> probably the interface doesn't work?!\n        return value\n\n\nclass SkyDiffuseMap(SkySpatialModel):\n    \"\"\"Spatial sky map template model.\n\n    At the moment only support 2D maps.\n    TODO: support maps with an energy axis here or in a separate class?\n    TODO: should we cache some interpolator object for efficiency?\n\n    Parameters\n    ----------\n    map : `~gammapy.map.Map`\n        Map template\n    norm : `~astropy.units.Quantity`\n        Norm parameter (multiplied with map values)\n    meta : dict, optional\n        Meta information, meta['filename'] will be used for serialization\n    \"\"\"\n\n    def __init__(self, map, norm=1, meta=None):\n        self._map = map\n        self.parameters = ParameterList([\n            Parameter('norm', norm),\n        ])\n        self.meta = dict() if meta is None else meta\n\n    @classmethod\n    def read(cls, filename, **kwargs):\n        \"\"\"Read spatial template model from FITS image.\n\n        Parameters\n        ----------\n        filename : str\n            FITS image filename.\n        \"\"\"\n        template = Map.read(filename, **kwargs)\n        return cls(template)\n\n    def evaluate(self, lon, lat, norm):\n        coord = dict(\n            lon=lon.to('deg').value,\n            lat=lat.to('deg').value,\n        )\n        val = self._map.interp_by_coord(coord, fill_value=0)\n        # TODO: use map unit? self._map.unit\n        return norm * val * u.Unit('sr-1')\n", "path": "gammapy/image/models/new.py"}]}