diff --git a/ChangeLog b/ChangeLog
index 592a65a..2494534 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,8 +1,10 @@
-2018-02-08
+2018-02-13
 
         * Bug fix version 1.5.1 released
         ================================
 
+        Correct sky map smoothing.
+
 
 2018-01-23
 
diff --git a/NEWS b/NEWS
index 2cf4bbe..5f6facb 100644
--- a/NEWS
+++ b/NEWS
@@ -1,6 +1,6 @@
 New Features and Important Changes in GammaLib 1.5.1
 
-8 February 2018
+13 February 2018
 
 
 1. Introduction
@@ -83,7 +83,9 @@ None
 
 12. Sky module
 --------------
-None
+The smoothing kernel covered before only two sectors, now all four sectors
+are covered. The smoothing now also works correctly in case that more than
+a single sky map is comprised in a GSkyMap object.
 
 
 13. Support module
diff --git a/doc/source/admin/release_history/1.5.1.rst b/doc/source/admin/release_history/1.5.1.rst
index 956769b..fb83d44 100644
--- a/doc/source/admin/release_history/1.5.1.rst
+++ b/doc/source/admin/release_history/1.5.1.rst
@@ -10,7 +10,8 @@ GammaLib 1.5.1 is a bug fix release for GammaLib release 1.5.
 
 The following bugs were fixed:
 
-* ...
+* Correct sky map smoothing
+
 
 The following improvements were made:
 
diff --git a/src/sky/GSkyMap.cpp b/src/sky/GSkyMap.cpp
index cc2ee9b..25194ad 100644
--- a/src/sky/GSkyMap.cpp
+++ b/src/sky/GSkyMap.cpp
@@ -1866,11 +1866,11 @@ void GSkyMap::smooth(const std::string& kernel, const double& par)
         GFft fft_kernel(smooth_kernel(kernel, par));
 
         // Loop over all sky maps
-        for (int i = 0; i < m_num_maps; ++i) {
+        for (int k = 0; k < m_num_maps; ++k) {
 
             // Extract sky map
             GNdarray array(m_num_x, m_num_y);
-            const double *src = m_pixels.data() + i*m_num_pixels;
+            const double *src = m_pixels.data() + k*m_num_pixels;
             double       *dst = array.data();
             for (int i = 0; i < m_num_pixels; ++i) {
                 *dst++ = *src++;
@@ -1882,12 +1882,12 @@ void GSkyMap::smooth(const std::string& kernel, const double& par)
             // Multiply FFT of sky map with FFT of kernel
             GFft fft_smooth = fft_array * fft_kernel;
 
-            // Background transform sky map
+            // Backward transform sky map
             GNdarray smooth = fft_smooth.backward();
 
             // Insert sky map
             src = smooth.data();
-            dst = m_pixels.data() + i*m_num_pixels;
+            dst = m_pixels.data() + k*m_num_pixels;
             for (int i = 0; i < m_num_pixels; ++i) {
                 *dst++ = *src++;
             }
@@ -3187,7 +3187,7 @@ GNdarray GSkyMap::smooth_kernel(const std::string& kernel,
     const GWcs* wcs = dynamic_cast<const GWcs*>(m_proj);
     if (wcs == NULL) {
         std::string msg = "Sky map is not a WCS projection. Method is only "
-                          "valid for WCS projects.";
+                          "valid for WCS projections.";
         throw GException::invalid_argument(G_SMOOTH_KERNEL, msg);
     }
 
@@ -3213,6 +3213,14 @@ GNdarray GSkyMap::smooth_kernel(const std::string& kernel,
                 if (std::sqrt(xqs + y*y) <= par) {
                     kern(ix1,iy1) += 1.0;
                     sum           += 1.0;
+                    if (ix2 < m_num_x) {
+                        kern(ix2,iy1) += 1.0;
+                        sum           += 1.0;
+                    }
+                    if (iy2 < m_num_y) {
+                        kern(ix1,iy2) += 1.0;
+                        sum           += 1.0;
+                    }
                     if ((ix2 < m_num_x) && (iy2 < m_num_y)) {
                         kern(ix2,iy2) += 1.0;
                         sum           += 1.0;
@@ -3238,6 +3246,14 @@ GNdarray GSkyMap::smooth_kernel(const std::string& kernel,
                 double value   = std::exp(norm * (xqs + y*y));
                 kern(ix1,iy1) += value;
                 sum           += value;
+                if (ix2 < m_num_x) {
+                    kern(ix2,iy1) += value;
+                    sum           += value;
+                }
+                if (iy2 < m_num_y) {
+                    kern(ix1,iy2) += value;
+                    sum           += value;
+                }
                 if ((ix2 < m_num_x) && (iy2 < m_num_y)) {
                     kern(ix2,iy2) += value;
                     sum           += value;
diff --git a/test/test_GSky.cpp b/test/test_GSky.cpp
index fc20424..683a84e 100644
--- a/test/test_GSky.cpp
+++ b/test/test_GSky.cpp
@@ -1212,6 +1212,7 @@ void TestGSky::test_GSkyMap(void)
     // Smooth map using DISK kernel
     GSkyMap map_smooth1 = map_smooth;
     map_smooth1.smooth("DISK", 0.2);
+    map_smooth1.save("test_map_smooth_disk.fits", true);
 
     // Compute sum of smoothed sky map
     double sum_smooth1 = 0.0;
@@ -1223,6 +1224,7 @@ void TestGSky::test_GSkyMap(void)
     // Smooth map using GAUSSIAN kernel
     GSkyMap map_smooth2 = map_smooth;
     map_smooth2.smooth("GAUSSIAN", 0.2);
+    map_smooth2.save("test_map_smooth_gaussian.fits", true);
 
     // Compute sum of smoothed sky map
     double sum_smooth2 = 0.0;