diff --git a/gammapy/image/models/core.py b/gammapy/image/models/core.py
index c1364f80e3..ee406336e0 100644
--- a/gammapy/image/models/core.py
+++ b/gammapy/image/models/core.py
@@ -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
@@ -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:
@@ -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
     ----------
@@ -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
@@ -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`).
@@ -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):
diff --git a/gammapy/image/models/tests/test_core.py b/gammapy/image/models/tests/test_core.py
index 266d5b1341..4a031e82e3 100644
--- a/gammapy/image/models/tests/test_core.py
+++ b/gammapy/image/models/tests/test_core.py
@@ -8,6 +8,7 @@
 from ..core import (
     SkyPointSource,
     SkyGaussian,
+    SkyGaussianElongated,
     SkyDisk,
     SkyEllipse,
     SkyShell,
@@ -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)