diff --git a/CHANGELOG b/CHANGELOG
index 377e99ce1bcb92ab0733a64b1c844177516041b5..3a62f998bde9dae25d200c8cb73a835018c947b1 100644
--- a/CHANGELOG
+++ b/CHANGELOG
@@ -4,6 +4,7 @@ BornAgain-20.1, released 2023.04.28
     * Correct prefactor for scattering from rough interfaces (k^4, not k^2, #553)
     * Throw error if rough interfaces are combined with polarized simulation (#564)
     * Repair GUI support for angular grazing scans (#574)
+    * Consistently mention K-correlation model which is really used for roughness (#581)
   > Dependencies
     * Adapt to libheinz-1.1.0
 
diff --git a/Sample/Interface/LayerRoughness.cpp b/Sample/Interface/LayerRoughness.cpp
index eb66aec285abb4d7e91ed29df48728f9eb8a42ac..42e026d21d76906c16d49b771611f58806d4437c 100644
--- a/Sample/Interface/LayerRoughness.cpp
+++ b/Sample/Interface/LayerRoughness.cpp
@@ -16,6 +16,7 @@
 #include "Base/Math/Constants.h"
 #include "Base/Py/PyFmt.h"
 #include "Base/Util/Assert.h"
+#include <boost/math/special_functions/bessel.hpp>
 
 //! Constructor of layer roughness.
 //! @param sigma: rms of the roughness in nanometers
@@ -40,8 +41,7 @@ LayerRoughness::LayerRoughness()
 //! Fourier transform of the correlation function of the roughness profile.
 //!
 //! Based on the article
-//! D.K.G. de Boer, Physical review B, Volume 51, Number 8, 15 February 1995
-//! "X-ray reflection and transmission by rough surfaces"
+//! Palasantzas, Phys Rev B, 48, 14472 (1993)
 /* ************************************************************************* */
 double LayerRoughness::spectralFunction(const R3 kvec) const
 {
@@ -54,14 +54,14 @@ double LayerRoughness::spectralFunction(const R3 kvec) const
 }
 
 //! Correlation function of the roughness profile
-
 double LayerRoughness::corrFunction(const R3 k) const
 {
     ASSERT(m_validated);
     double H = m_hurstParameter;
     double clength = m_lateralCorrLength;
     double R = sqrt(k.x() * k.x() + k.y() * k.y());
-    return m_sigma * m_sigma * std::exp(-1.0 * std::pow(R / clength, 2. * H));
+    return m_sigma * m_sigma * std::pow(2., 1. - H) / tgamma(H) * std::pow(R / clength, H)
+           * boost::math::cyl_bessel_k(H, R / clength);
 }
 
 std::string LayerRoughness::pythonConstructor() const
diff --git a/hugo/content/ref/sample/roughness/scattering/index.md b/hugo/content/ref/sample/roughness/scattering/index.md
index 38512c46288ef079bf962c2d8b540c3eb850358d..bde0d2c0d1b5bfbc8ad24b0f407772c759fd4599 100644
--- a/hugo/content/ref/sample/roughness/scattering/index.md
+++ b/hugo/content/ref/sample/roughness/scattering/index.md
@@ -21,7 +21,7 @@ Scattering from a multilayered sample with correlated roughness.
 
 **Note:**
 
-The roughness profile is described by a normally-distributed random function. The roughness correlation function at the jth interface is expressed as: $$ < U\_j (x, y) U\_j (x', y')> = \sigma^2 e^{-\frac{\tau}{ξ}2H}, \tau=[(x-x')^2+(y-y')^2]^{\frac{1}{2}}$$
+The roughness profile is described by a normally-distributed random function. The roughness correlation function at the jth interface is expressed as: $$ < U\_j (x, y) U\_j (x', y')> = \sigma^2 \frac{2^{1-H}}{\Gamma(H)} \left( \frac{\tau}{ξ} \right)^H K_H \left( \frac{\tau}{ξ} \right), \tau=[(x-x')^2+(y-y')^2]^{\frac{1}{2}}$$
 
 * $U\_j(x, y)$ is the height deviation of the jth interface at position $(x, y)$.
 * $\sigma$ gives the rms roughness of the interface. The Hurst parameter $H$, comprised between $0$ and $1$ is connected to the fractal dimension $D=3-H$ of the interface. The smaller $H$ is, the more serrate the surface profile looks. If $H = 1$, the interface has a non fractal nature.