From 741a77da41ddb7943f7aa6cf4af336364605ee2e Mon Sep 17 00:00:00 2001
From: Mikhail Svechnikov <m.svechnikov@fz-juelich.de>
Date: Thu, 29 Aug 2024 16:02:43 +0200
Subject: [PATCH] cache rms

---
 GUI/View/Realspace/RealspaceBuilder.cpp       |   2 +-
 Resample/Processed/ReSample.cpp               |  16 +--
 Resample/Slice/ProfileHelper.cpp              |  19 ++--
 Resample/Slice/Slice.cpp                      |  18 +--
 Resample/Slice/Slice.h                        |  10 +-
 Resample/Slice/SliceStack.cpp                 |  17 ++-
 Resample/Slice/SliceStack.h                   |   5 +-
 Resample/Specular/ComputeFluxMagnetic.cpp     |   8 +-
 Resample/Specular/ComputeFluxScalar.cpp       |  18 +--
 Sample/Interface/AutocorrelationModels.cpp    |   6 -
 Sample/Interface/AutocorrelationModels.h      |   8 +-
 Sample/Interface/LayerRoughness.h             |   2 -
 Sample/Interface/RoughnessMap.cpp             |   4 +-
 .../RoughMultiLayerContribution.cpp           |   7 +-
 Tests/Unit/Sample/MultiLayerTest.cpp          |  20 ++--
 auto/Wrap/libBornAgainSample.py               |  12 --
 auto/Wrap/libBornAgainSample_wrap.cpp         | 105 ------------------
 17 files changed, 72 insertions(+), 205 deletions(-)

diff --git a/GUI/View/Realspace/RealspaceBuilder.cpp b/GUI/View/Realspace/RealspaceBuilder.cpp
index 2a512a8bc30..568b074bece 100644
--- a/GUI/View/Realspace/RealspaceBuilder.cpp
+++ b/GUI/View/Realspace/RealspaceBuilder.cpp
@@ -105,7 +105,7 @@ std::unique_ptr<const double2d_t> layerRoughnessMap(const LayerItem& layerItem,
         roughItem->certainInterlayerModel()->createModel();
 
     auto roughness = LayerRoughness(autocorrelation.get(), interlayer.get());
-    if (roughness.rms() == 0)
+    if (true /*roughness.rms() == 0*/)
         return result;
 
     int n = sceneGeometry.numRoughnessPointsAlongAxis;
diff --git a/Resample/Processed/ReSample.cpp b/Resample/Processed/ReSample.cpp
index 2d03aa05aaf..9a8a71b01b5 100644
--- a/Resample/Processed/ReSample.cpp
+++ b/Resample/Processed/ReSample.cpp
@@ -114,7 +114,6 @@ SliceStack slicify(const MultiLayer& sample, bool useAvgMaterials)
             const Layer* const layer = sample.layer(i);
             double tl = layer->thickness();
             const Material* const material = layer->material();
-            const LayerRoughness* roughness = layerTopRoughness(sample, i);
 
             if (i == 0) {
                 ASSERT(tl == 0);
@@ -122,7 +121,9 @@ SliceStack slicify(const MultiLayer& sample, bool useAvgMaterials)
             } else {
                 if (i == nLayers - 1)
                     ASSERT(tl == 0);
-                result.addSlice(tl, *material, roughness);
+                const LayerRoughness* roughness = layerTopRoughness(sample, i);
+                const double rms = sample.layerInterfaceRMS(i - 1);
+                result.addSlice(tl, *material, roughness, rms);
             }
         }
         return result;
@@ -134,8 +135,9 @@ SliceStack slicify(const MultiLayer& sample, bool useAvgMaterials)
         const Layer* const layer = sample.layer(i);
         double tl = layer->thickness();
         const Material* const material = layer->material();
-        const LayerRoughness* roughness = layerTopRoughness(sample, i);
         const ZLimits& particle_span = particle_spans[i];
+        const LayerRoughness* roughness = layerTopRoughness(sample, i);
+        const double rms = (i >= 1) ? sample.layerInterfaceRMS(i - 1) : 0.0;
 
         // if no slicing is needed, create single slice for the layer
         if (!particle_span.isFinite()) { // also if layer contains no particles
@@ -144,7 +146,7 @@ SliceStack slicify(const MultiLayer& sample, bool useAvgMaterials)
             if (i == 0)
                 result.addTopSlice(tl, *material);
             else
-                result.addSlice(tl, *material, roughness);
+                result.addSlice(tl, *material, roughness, rms);
             continue;
         }
 
@@ -165,10 +167,10 @@ SliceStack slicify(const MultiLayer& sample, bool useAvgMaterials)
         else {
             ASSERT(top <= 0);
             if (top < 0) {
-                result.addSlice(-top, *material, roughness);
+                result.addSlice(-top, *material, roughness, rms);
                 result.addNSlices(nSlices, top - bottom, *material);
             } else { // top == 0
-                result.addNSlices(nSlices, top - bottom, *material, roughness);
+                result.addNSlices(nSlices, top - bottom, *material, roughness, rms);
             }
             // middle layer
             if (i < nLayers - 1 && bottom > -tl)
@@ -402,7 +404,7 @@ bool ReSample::polarizing() const
 bool ReSample::hasRoughness() const
 {
     for (const auto& slice : m_stack)
-        if (slice.topRoughness() && slice.topRoughness()->rms() > 0)
+        if (slice.topRMS() > 0)
             return true;
     return false;
 }
diff --git a/Resample/Slice/ProfileHelper.cpp b/Resample/Slice/ProfileHelper.cpp
index c969445e83c..0d3701625f4 100644
--- a/Resample/Slice/ProfileHelper.cpp
+++ b/Resample/Slice/ProfileHelper.cpp
@@ -22,12 +22,12 @@ using std::numbers::pi;
 
 namespace {
 
-double Transition(double x, const LayerRoughness* roughness)
+double Transition(double x, const LayerRoughness* roughness, double rms)
 {
     if (!roughness)
         return Math::Heaviside(x);
 
-    return roughness->interlayerModel()->transient(x, roughness->rms());
+    return roughness->interlayerModel()->transient(x, rms);
 }
 
 const std::string SLD = "SLD";
@@ -77,7 +77,7 @@ std::vector<complex_t> ProfileHelper::profile(const std::vector<double>& z_value
             quantity(slice.material(), component) - quantity(sliceAbove.material(), component);
         for (size_t j = 0; j < z_values.size(); ++j) {
             const double arg = (slice.hig() - z_values[j]);
-            const double t = Transition(arg, slice.topRoughness());
+            const double t = Transition(arg, slice.topRoughness(), slice.topRMS());
             result[j] += diff * t;
         }
     }
@@ -109,14 +109,11 @@ std::pair<double, double> ProfileHelper::defaultLimits() const
         return {0.0, 0.0};
     double interface_span = m_stack.front().low() - m_stack.back().hig();
     double default_margin = interface_span > 0.0 ? interface_span / 20.0 : 10.0;
-    const LayerRoughness* topRoughness = m_stack.at(1).topRoughness();
-    const LayerRoughness* bottomRoughness = m_stack.back().topRoughness();
-
-    double top_margin =
-        topRoughness && topRoughness->rms() > 0 ? 5.0 * topRoughness->rms() : default_margin;
-    double bottom_margin = bottomRoughness && bottomRoughness->rms() > 0
-                               ? 5.0 * bottomRoughness->rms()
-                               : default_margin;
+    const double top_rms = m_stack.at(1).topRMS();
+    const double bottom_rms = m_stack.back().topRMS();
+
+    double top_margin = top_rms > 0 ? 5.0 * top_rms : default_margin;
+    double bottom_margin = bottom_rms > 0 ? 5.0 * bottom_rms : default_margin;
     double z_min = m_stack.back().hig() - bottom_margin;
     double z_max = m_stack.front().low() + top_margin;
     return {z_min, z_max};
diff --git a/Resample/Slice/Slice.cpp b/Resample/Slice/Slice.cpp
index ef2e4838f0d..6200188998a 100644
--- a/Resample/Slice/Slice.cpp
+++ b/Resample/Slice/Slice.cpp
@@ -23,26 +23,17 @@
 using std::numbers::pi;
 
 Slice::Slice(const ZLimits& zRange, const Material& material, const R3& B_field,
-             const LayerRoughness* const roughness)
+             const LayerRoughness* const roughness, double rms)
     : m_z_range(zRange)
     , m_material(material)
     , m_B_field(B_field)
     , m_top_roughness(roughness)
+    , m_top_rms(rms)
 {
 }
 
 Slice::~Slice() = default;
 
-void Slice::setMaterial(const Material& material)
-{
-    m_material = material;
-}
-
-const Material& Slice::material() const
-{
-    return m_material;
-}
-
 double Slice::low() const
 {
     return m_z_range.low();
@@ -63,11 +54,6 @@ double Slice::thicknessOr0() const
     return m_z_range.thicknessOr0();
 }
 
-const LayerRoughness* Slice::topRoughness() const
-{
-    return m_top_roughness;
-}
-
 complex_t Slice::scalarReducedPotential(R3 k, double n_ref) const
 {
     complex_t n = m_material.refractiveIndex((2 * pi) / k.mag());
diff --git a/Resample/Slice/Slice.h b/Resample/Slice/Slice.h
index c956cfad780..d18b5f0cc50 100644
--- a/Resample/Slice/Slice.h
+++ b/Resample/Slice/Slice.h
@@ -29,18 +29,19 @@ class LayerRoughness;
 class Slice {
 public:
     Slice(const ZLimits& zRange, const Material& material, const R3& B_field,
-          const LayerRoughness* roughness);
+          const LayerRoughness* roughness, double rms);
     ~Slice();
 
-    void setMaterial(const Material& material);
-    const Material& material() const;
+    void setMaterial(const Material& material) { m_material = material; }
+    const Material& material() const { return m_material; }
 
     double low() const;
     double hig() const;
     double higOr0() const;
     const ZLimits& span() const { return m_z_range; }
     double thicknessOr0() const;
-    const LayerRoughness* topRoughness() const;
+    const LayerRoughness* topRoughness() const { return m_top_roughness; }
+    double topRMS() const { return m_top_rms; }
 
     //! Return the potential term that is used in the one-dimensional Fresnel calculations
     complex_t scalarReducedPotential(R3 k, double n_ref) const;
@@ -60,6 +61,7 @@ private:
     Material m_material;
     R3 m_B_field; //!< cached value of magnetic induction
     const LayerRoughness* const m_top_roughness;
+    double m_top_rms;
 };
 
 #endif // BORNAGAIN_RESAMPLE_SLICE_SLICE_H
diff --git a/Resample/Slice/SliceStack.cpp b/Resample/Slice/SliceStack.cpp
index 11ca87bc6db..9f07a0e432b 100644
--- a/Resample/Slice/SliceStack.cpp
+++ b/Resample/Slice/SliceStack.cpp
@@ -18,11 +18,11 @@
 
 void SliceStack::addTopSlice(double zbottom, const Material& material)
 {
-    this->emplace_back(Slice(ZLimits(zbottom, ZLimits::inf), material, {}, nullptr));
+    this->emplace_back(Slice(ZLimits(zbottom, ZLimits::inf), material, {}, nullptr, 0.0));
 }
 
 void SliceStack::addSlice(double thickness, const Material& material,
-                          const LayerRoughness* const roughness)
+                          const LayerRoughness* const roughness, double rms)
 {
     ASSERT(!this->empty());
     double top = this->back().low();
@@ -32,18 +32,18 @@ void SliceStack::addSlice(double thickness, const Material& material,
         zRange = std::make_unique<ZLimits>(-ZLimits::inf, top);
     else
         zRange = std::make_unique<ZLimits>(top - thickness, top);
-    this->emplace_back(Slice(*zRange, material, {}, roughness));
+    this->emplace_back(Slice(*zRange, material, {}, roughness, rms));
 }
 
 //! Adds n times the same slice to the stack.
 
 void SliceStack::addNSlices(size_t n, double thickness, const Material& material,
-                            const LayerRoughness* const roughness)
+                            const LayerRoughness* const roughness, double rms)
 {
     ASSERT(thickness > 0);
     ASSERT(n > 0);
     const double slice_thickness = thickness / n;
-    addSlice(slice_thickness, material, roughness);
+    addSlice(slice_thickness, material, roughness, rms);
     for (size_t i = 1; i < n; ++i)
         addSlice(slice_thickness, material);
 }
@@ -73,3 +73,10 @@ const LayerRoughness* SliceStack::bottomRoughness(size_t i_slice) const
         return (*this)[i_slice + 1].topRoughness();
     return nullptr;
 }
+
+double SliceStack::bottomRMS(size_t i_slice) const
+{
+    if (i_slice + 1 < size())
+        return (*this)[i_slice + 1].topRMS();
+    return 0.0;
+}
diff --git a/Resample/Slice/SliceStack.h b/Resample/Slice/SliceStack.h
index 09c24c934f0..f765a70ed09 100644
--- a/Resample/Slice/SliceStack.h
+++ b/Resample/Slice/SliceStack.h
@@ -34,14 +34,15 @@ public:
 
     void addTopSlice(double zbottom, const Material& material);
     void addSlice(double thickness, const Material& material,
-                  const LayerRoughness* roughness = nullptr);
+                  const LayerRoughness* roughness = nullptr, double rms = 0);
     void addNSlices(size_t n, double thickness, const Material& material,
-                    const LayerRoughness* roughness = nullptr);
+                    const LayerRoughness* roughness = nullptr, double rms = 0);
 
     SliceStack setBField(const R3& externalField);
 
     // Top roughness of the slice below
     const LayerRoughness* bottomRoughness(size_t i_slice) const;
+    double bottomRMS(size_t i_slice) const;
 };
 
 #endif // BORNAGAIN_RESAMPLE_SLICE_SLICESTACK_H
diff --git a/Resample/Specular/ComputeFluxMagnetic.cpp b/Resample/Specular/ComputeFluxMagnetic.cpp
index 63b2380f5c5..1d02e028171 100644
--- a/Resample/Specular/ComputeFluxMagnetic.cpp
+++ b/Resample/Specular/ComputeFluxMagnetic.cpp
@@ -182,11 +182,11 @@ std::vector<MatrixFlux> computeTR(const SliceStack& slices, const std::vector<co
 
     for (size_t i = N - 1; i > 0; --i) {
         const auto* roughness = slices.bottomRoughness(i);
-        const double sigma = roughness ? roughness->rms() : 0.;
+        const double rms = slices.bottomRMS(i);
         const auto* r_model = roughness ? roughness->interlayerModel() : nullptr;
 
         // compute the 2x2 blocks of the transfer matrix
-        const auto [sp, sm] = refractionMatrixBlocks(TR[i - 1], TR[i], sigma, r_model);
+        const auto [sp, sm] = refractionMatrixBlocks(TR[i - 1], TR[i], rms, r_model);
         const SpinMatrix delta = TR[i - 1].computeDeltaMatrix(slices[i - 1].thicknessOr0());
 
         // compute the rotation matrix
@@ -288,11 +288,11 @@ SpinMatrix Compute::polarizedReflectivity(const SliceStack& slices,
     for (size_t i = N - 1; i > 0; --i) {
         MatrixFlux tr_i = createCoeff(i - 1);
         const auto* roughness = slices.bottomRoughness(i - 1);
-        const double sigma = roughness ? roughness->rms() : 0.;
+        const double rms = slices.bottomRMS(i - 1);
         const auto* r_model = roughness ? roughness->interlayerModel() : nullptr;
 
         // compute the 2x2 blocks of the transfer matrix
-        const auto [sp, sm] = ::refractionMatrixBlocks(tr_i, tr_i1, sigma, r_model);
+        const auto [sp, sm] = ::refractionMatrixBlocks(tr_i, tr_i1, rms, r_model);
         const SpinMatrix delta = tr_i.computeDeltaMatrix(slices[i - 1].thicknessOr0());
 
         // compute the rotation matrix
diff --git a/Resample/Specular/ComputeFluxScalar.cpp b/Resample/Specular/ComputeFluxScalar.cpp
index 8d40208f837..07025f78b41 100644
--- a/Resample/Specular/ComputeFluxScalar.cpp
+++ b/Resample/Specular/ComputeFluxScalar.cpp
@@ -27,15 +27,15 @@ using std::numbers::pi;
 namespace {
 
 //! See PhysRef, chapter "Scattering by rough interfaces", section "Interface with tanh profile".
-std::pair<complex_t, complex_t> transition(complex_t kzi, complex_t kzi1, double sigma,
+std::pair<complex_t, complex_t> transition(complex_t kzi, complex_t kzi1, double rms,
                                            const InterlayerModel* r_model)
 {
     const complex_t kz_ratio = kzi1 / kzi;
 
-    if (sigma == 0)
+    if (rms == 0)
         return {1. + kz_ratio, 1. - kz_ratio};
 
-    ASSERT(sigma > 0);
+    ASSERT(rms > 0);
     ASSERT(r_model);
 
     if (dynamic_cast<const ErfInterlayer*>(r_model))
@@ -43,14 +43,14 @@ std::pair<complex_t, complex_t> transition(complex_t kzi, complex_t kzi1, double
         // reflection coefficients.
         // Implementation follows A. Gibaud and G. Vignaud, in X-ray and Neutron Reflectivity,
         // edited by J. Daillant and A. Gibaud, volume 770 of Lecture Notes in Physics (2009)
-        return {(1. + kz_ratio) * std::exp(-pow((kzi1 - kzi) * sigma, 2) / 2.),
-                (1. - kz_ratio) * std::exp(-pow((kzi1 + kzi) * sigma, 2) / 2.)};
+        return {(1. + kz_ratio) * std::exp(-pow((kzi1 - kzi) * rms, 2) / 2.),
+                (1. - kz_ratio) * std::exp(-pow((kzi1 + kzi) * rms, 2) / 2.)};
 
     // TANH model assumed
     // Roughness is modelled by tanh profile
     // [e.g. Bahr, Press, et al, Phys. Rev. B, vol. 47 (8), p. 4385 (1993)].
     // but the scale factor is adjusted to give rms == sigma.
-    const double sigeff = std::sqrt(3.0) * sigma;
+    const double sigeff = std::sqrt(3.0) * rms;
     const complex_t roughness = std::sqrt(Math::tanhc(sigeff * kzi1) / Math::tanhc(sigeff * kzi));
 
     return {1. / roughness + kz_ratio * roughness, 1. / roughness - kz_ratio * roughness};
@@ -84,7 +84,7 @@ std::vector<Spinor> computeTR(const SliceStack& slices, const std::vector<comple
         const size_t jthis = X[i - 1];
         const size_t jlast = X[i];
         const auto* roughness = slices.bottomRoughness(jthis);
-        const double rms = roughness ? roughness->rms() : 0.;
+        const double rms = slices.bottomRMS(jthis);
         const auto* r_model = roughness ? roughness->interlayerModel() : nullptr;
 
         const auto [slp, slm] = transition(kz[jthis], kz[jlast], rms, r_model);
@@ -138,10 +138,10 @@ complex_t Compute::scalarReflectivity(const SliceStack& slices, const std::vecto
 
     for (size_t i = N - 1; i > 0; i--) {
         const auto* roughness = slices.bottomRoughness(i - 1);
-        const double sigma = roughness ? roughness->rms() : 0.;
+        const double rms = slices.bottomRMS(i - 1);
         const auto* r_model = roughness ? roughness->interlayerModel() : nullptr;
 
-        const auto [sp, sm] = ::transition(kz[i - 1], kz[i], sigma, r_model);
+        const auto [sp, sm] = ::transition(kz[i - 1], kz[i], rms, r_model);
 
         const complex_t delta = exp_I(kz[i - 1] * slices[i - 1].thicknessOr0());
 
diff --git a/Sample/Interface/AutocorrelationModels.cpp b/Sample/Interface/AutocorrelationModels.cpp
index 580a70f0a03..00dddd74f62 100644
--- a/Sample/Interface/AutocorrelationModels.cpp
+++ b/Sample/Interface/AutocorrelationModels.cpp
@@ -218,9 +218,3 @@ double LinearGrowthModel::spectralFunction(double) const
     ASSERT(m_validated);
     return 0;
 }
-
-double LinearGrowthModel::rms() const
-{
-    // TODO
-    return 0.001;
-}
diff --git a/Sample/Interface/AutocorrelationModels.h b/Sample/Interface/AutocorrelationModels.h
index e43d5982743..45a0b81e79f 100644
--- a/Sample/Interface/AutocorrelationModels.h
+++ b/Sample/Interface/AutocorrelationModels.h
@@ -24,9 +24,6 @@ class AutocorrelationModel : public ICloneable, public INode {
 public:
     AutocorrelationModel* clone() const override = 0;
 
-    //! Returns rms of roughness
-    virtual double rms() const = 0;
-
     //! Returns power spectral density of the surface roughness
     virtual double spectralFunction(double spatial_f) const = 0;
 
@@ -58,7 +55,9 @@ public:
     std::vector<ParaMeta> parDefs() const final;
     std::string validate() const override;
     double spectralFunction(double spatial_f) const override;
-    double rms() const override;
+
+    //! Returns rms of roughness
+    double rms() const;
 
     void setSigma(double sigma) { m_sigma = sigma; }
     double sigma() const { return m_sigma; }
@@ -95,7 +94,6 @@ public:
     //! Returns power spectral density of the surface roughness
     double spectralFunction(double spectrum_below, double thickness, double spatial_f) const;
     double spectralFunction(double spatial_f) const override;
-    double rms() const override;
 
     double crosscorrSpectrum(double spectrum_below, double thickness, double spatial_f) const;
 
diff --git a/Sample/Interface/LayerRoughness.h b/Sample/Interface/LayerRoughness.h
index 1105566edb2..ae08b1f8b0b 100644
--- a/Sample/Interface/LayerRoughness.h
+++ b/Sample/Interface/LayerRoughness.h
@@ -35,8 +35,6 @@ public:
     std::string className() const final { return "LayerRoughness"; }
     std::vector<const INode*> nodeChildren() const override;
 
-    double rms() const { return m_autocorrelation_model->rms(); }
-
     const AutocorrelationModel* autocorrelationModel() const
     {
         return m_autocorrelation_model.get();
diff --git a/Sample/Interface/RoughnessMap.cpp b/Sample/Interface/RoughnessMap.cpp
index f419de63a6d..030a760e85e 100644
--- a/Sample/Interface/RoughnessMap.cpp
+++ b/Sample/Interface/RoughnessMap.cpp
@@ -70,7 +70,7 @@ double2d_t RoughnessMap::generateMap()
 double2d_t RoughnessMap::mapFromHeights() const
 {
     const size_t z_steps = 3999;
-    const double rms = m_layerRoughness->rms();
+    const double rms = 0; // m_layerRoughness->rms(); TODO
     const double sigma_factor = m_layerRoughness->interlayerModel()->sigmaRange();
     const double z_limit = rms * sigma_factor;
     const double step = 2 * z_limit / (z_steps - 1);
@@ -196,7 +196,7 @@ double2d_t RoughnessMap::applyHeightsToSpectrum(const double2d_t& h_map,
 
 void RoughnessMap::createMap()
 {
-    if (m_layerRoughness->rms() < 1e-10) {
+    if (true /*m_layerRoughness->rms() < 1e-10*/) {
         m_rough_map = double2d_t(m_y_points, std::vector<double>(m_x_points));
         return;
     }
diff --git a/Sim/Computation/RoughMultiLayerContribution.cpp b/Sim/Computation/RoughMultiLayerContribution.cpp
index 8148825e87d..e780d49c09a 100644
--- a/Sim/Computation/RoughMultiLayerContribution.cpp
+++ b/Sim/Computation/RoughMultiLayerContribution.cpp
@@ -86,8 +86,7 @@ complex_t get_sum8terms(const ReSample& re_sample, size_t i_layer, const Diffuse
     const complex_t qz3_B = -qz2_B;
     const complex_t qz4_B = -qz1_B;
 
-    const LayerRoughness* roughness = re_sample.averageSlices().bottomRoughness(i_layer);
-    const double rms = roughness ? roughness->rms() : 0.;
+    const double rms = re_sample.averageSlices().bottomRMS(i_layer);
     const complex_t term1 = T_i_A * T_f_A * ::h_above(qz1_A * rms);
     const complex_t term2 = T_i_A * R_f_A * ::h_above(qz2_A * rms);
     const complex_t term3 = R_i_A * T_f_A * ::h_above(qz3_A * rms);
@@ -153,8 +152,8 @@ double Compute::roughMultiLayerContribution(const ReSample& re_sample, const Dif
         const double thickness = roughStack[j].hig() - roughStack[j + 1].hig();
         if (auto* spat_freq_cc = dynamic_cast<const SpatialFrequencyCrosscorrelation*>(
                 roughStack[j].topRoughness()->crosscorrelationModel())) {
-            const double rms_up = roughStack[j].topRoughness()->rms();
-            const double rms_low = roughStack[j + 1].topRoughness()->rms();
+            const double rms_up = roughStack[j].topRMS();
+            const double rms_low = roughStack[j + 1].topRMS();
             crosscorr_with_interface_below[j] = spat_freq_cc->crosscorrSpectrum(
                 rms_up, rms_low, spectrum[j], spectrum[j + 1], thickness, spatial_f);
         } else if (auto* lin_growth_autocorr = dynamic_cast<const LinearGrowthModel*>(
diff --git a/Tests/Unit/Sample/MultiLayerTest.cpp b/Tests/Unit/Sample/MultiLayerTest.cpp
index 2b6248f0ec7..224772063ce 100644
--- a/Tests/Unit/Sample/MultiLayerTest.cpp
+++ b/Tests/Unit/Sample/MultiLayerTest.cpp
@@ -87,17 +87,17 @@ TEST_F(MultiLayerTest, LayerInterfaces)
     const LayerInterface* interface0 = mLayer.layerInterface(0);
     EXPECT_TRUE(nullptr != interface0);
     EXPECT_TRUE(nullptr != interface0->roughness());
-    EXPECT_EQ(interface0->roughness()->rms(), 0.0);
+    EXPECT_EQ(mLayer.layerInterfaceRMS(0), 0.0);
 
     const LayerInterface* interface1 = mLayer.layerInterface(1);
     EXPECT_TRUE(nullptr != interface1);
     EXPECT_TRUE(nullptr != interface1->roughness());
-    EXPECT_EQ(interface1->roughness()->rms(), 0.0);
+    EXPECT_EQ(mLayer.layerInterfaceRMS(1), 0.0);
 
     const LayerInterface* interface2 = mLayer.layerInterface(2);
     EXPECT_TRUE(nullptr != interface2);
     EXPECT_TRUE(nullptr != interface2->roughness());
-    EXPECT_EQ(interface2->roughness()->rms(), 0.0);
+    EXPECT_EQ(mLayer.layerInterfaceRMS(2), 0.0);
 }
 
 TEST_F(MultiLayerTest, Clone)
@@ -119,19 +119,19 @@ TEST_F(MultiLayerTest, Clone)
     const LayerInterface* interface0 = mLayerClone->layerInterface(0);
     EXPECT_TRUE(nullptr != interface0);
     EXPECT_TRUE(nullptr != interface0->roughness());
-    EXPECT_EQ(interface0->roughness()->rms(), 0.0);
+    EXPECT_EQ(mLayerClone->layerInterfaceRMS(0), 0.0);
     EXPECT_TRUE(nullptr == interface0->roughness()->crosscorrelationModel());
 
     const LayerInterface* interface1 = mLayerClone->layerInterface(1);
     EXPECT_TRUE(nullptr != interface1);
     EXPECT_TRUE(nullptr != interface1->roughness());
-    EXPECT_EQ(interface1->roughness()->rms(), 0.0);
+    EXPECT_EQ(mLayerClone->layerInterfaceRMS(1), 0.0);
     EXPECT_TRUE(nullptr == interface1->roughness()->crosscorrelationModel());
 
     const LayerInterface* interface2 = mLayerClone->layerInterface(2);
     EXPECT_TRUE(nullptr != interface2);
     EXPECT_TRUE(nullptr != interface2->roughness());
-    EXPECT_EQ(interface2->roughness()->rms(), 0.0);
+    EXPECT_EQ(mLayerClone->layerInterfaceRMS(2), 0.0);
     EXPECT_TRUE(nullptr == interface2->roughness()->crosscorrelationModel());
 
     delete mLayerClone;
@@ -157,12 +157,12 @@ TEST_F(MultiLayerTest, WithRoughness)
         dynamic_cast<const K_CorrelationModel*>(roughness0->autocorrelationModel());
     EXPECT_TRUE(roughness0_AC);
 
-    EXPECT_EQ(1.1, roughness0->rms());
+    EXPECT_EQ(1.1, mLayer.layerInterfaceRMS(0));
     EXPECT_EQ(.3, roughness0_AC->hurst());
     EXPECT_EQ(0.1, roughness0_AC->lateralCorrLength());
 
     EXPECT_TRUE(roughness1);
-    EXPECT_EQ(roughness1->rms(), 0.0);
+    EXPECT_EQ(mLayer.layerInterfaceRMS(1), 0.0);
 }
 
 TEST_F(MultiLayerTest, CloneWithRoughness)
@@ -196,7 +196,7 @@ TEST_F(MultiLayerTest, CloneWithRoughness)
         dynamic_cast<const K_CorrelationModel*>(roughness0->autocorrelationModel());
     EXPECT_TRUE(roughness0_AC);
 
-    EXPECT_EQ(2.1, roughness0->rms());
+    EXPECT_EQ(2.1, mLayerClone->layerInterfaceRMS(0));
     EXPECT_EQ(.3, roughness0_AC->hurst());
     EXPECT_EQ(12.1, roughness0_AC->lateralCorrLength());
 
@@ -204,7 +204,7 @@ TEST_F(MultiLayerTest, CloneWithRoughness)
         dynamic_cast<const K_CorrelationModel*>(roughness1->autocorrelationModel());
     EXPECT_TRUE(roughness1_AC);
 
-    EXPECT_EQ(1.1, roughness1->rms());
+    EXPECT_EQ(1.1, mLayerClone->layerInterfaceRMS(1));
     EXPECT_EQ(.3, roughness1_AC->hurst());
     EXPECT_EQ(0.1, roughness1_AC->lateralCorrLength());
 
diff --git a/auto/Wrap/libBornAgainSample.py b/auto/Wrap/libBornAgainSample.py
index ab139647663..b20e0cda6ab 100644
--- a/auto/Wrap/libBornAgainSample.py
+++ b/auto/Wrap/libBornAgainSample.py
@@ -4191,10 +4191,6 @@ class LayerRoughness(ISampleNode):
         r"""nodeChildren(LayerRoughness self) -> swig_dummy_type_const_inode_vector"""
         return _libBornAgainSample.LayerRoughness_nodeChildren(self)
 
-    def rms(self):
-        r"""rms(LayerRoughness self) -> double"""
-        return _libBornAgainSample.LayerRoughness_rms(self)
-
     def autocorrelationModel(self):
         r"""autocorrelationModel(LayerRoughness self) -> AutocorrelationModel"""
         return _libBornAgainSample.LayerRoughness_autocorrelationModel(self)
@@ -4322,10 +4318,6 @@ class AutocorrelationModel(libBornAgainBase.ICloneable, libBornAgainParam.INode)
         r"""clone(AutocorrelationModel self) -> AutocorrelationModel"""
         return _libBornAgainSample.AutocorrelationModel_clone(self)
 
-    def rms(self):
-        r"""rms(AutocorrelationModel self) -> double"""
-        return _libBornAgainSample.AutocorrelationModel_rms(self)
-
     def spectralFunction(self, spatial_f):
         r"""spectralFunction(AutocorrelationModel self, double spatial_f) -> double"""
         return _libBornAgainSample.AutocorrelationModel_spectralFunction(self, spatial_f)
@@ -4439,10 +4431,6 @@ class LinearGrowthModel(AutocorrelationModel):
         """
         return _libBornAgainSample.LinearGrowthModel_spectralFunction(self, *args)
 
-    def rms(self):
-        r"""rms(LinearGrowthModel self) -> double"""
-        return _libBornAgainSample.LinearGrowthModel_rms(self)
-
     def crosscorrSpectrum(self, spectrum_below, thickness, spatial_f):
         r"""crosscorrSpectrum(LinearGrowthModel self, double spectrum_below, double thickness, double spatial_f) -> double"""
         return _libBornAgainSample.LinearGrowthModel_crosscorrSpectrum(self, spectrum_below, thickness, spatial_f)
diff --git a/auto/Wrap/libBornAgainSample_wrap.cpp b/auto/Wrap/libBornAgainSample_wrap.cpp
index d1d5c85c24f..d054eb2fe13 100644
--- a/auto/Wrap/libBornAgainSample_wrap.cpp
+++ b/auto/Wrap/libBornAgainSample_wrap.cpp
@@ -50739,40 +50739,6 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_LayerRoughness_rms(PyObject *self, PyObject *args) {
-  PyObject *resultobj = 0;
-  LayerRoughness *arg1 = (LayerRoughness *) 0 ;
-  void *argp1 = 0 ;
-  int res1 = 0 ;
-  PyObject *swig_obj[1] ;
-  double result;
-  
-  (void)self;
-  if (!args) SWIG_fail;
-  swig_obj[0] = args;
-  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_LayerRoughness, 0 |  0 );
-  if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "LayerRoughness_rms" "', argument " "1"" of type '" "LayerRoughness const *""'"); 
-  }
-  arg1 = reinterpret_cast< LayerRoughness * >(argp1);
-  {
-    try {
-      result = (double)((LayerRoughness const *)arg1)->rms();
-    } catch (const std::exception& ex) {
-      // message shown in the Python interpreter
-      const std::string msg {
-        "BornAgain C++ Exception: " + std::string(ex.what())
-      };
-      SWIG_exception(SWIG_RuntimeError, msg.c_str());
-    }
-  }
-  resultobj = SWIG_From_double(static_cast< double >(result));
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
 SWIGINTERN PyObject *_wrap_LayerRoughness_autocorrelationModel(PyObject *self, PyObject *args) {
   PyObject *resultobj = 0;
   LayerRoughness *arg1 = (LayerRoughness *) 0 ;
@@ -52051,40 +52017,6 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_AutocorrelationModel_rms(PyObject *self, PyObject *args) {
-  PyObject *resultobj = 0;
-  AutocorrelationModel *arg1 = (AutocorrelationModel *) 0 ;
-  void *argp1 = 0 ;
-  int res1 = 0 ;
-  PyObject *swig_obj[1] ;
-  double result;
-  
-  (void)self;
-  if (!args) SWIG_fail;
-  swig_obj[0] = args;
-  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_AutocorrelationModel, 0 |  0 );
-  if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "AutocorrelationModel_rms" "', argument " "1"" of type '" "AutocorrelationModel const *""'"); 
-  }
-  arg1 = reinterpret_cast< AutocorrelationModel * >(argp1);
-  {
-    try {
-      result = (double)((AutocorrelationModel const *)arg1)->rms();
-    } catch (const std::exception& ex) {
-      // message shown in the Python interpreter
-      const std::string msg {
-        "BornAgain C++ Exception: " + std::string(ex.what())
-      };
-      SWIG_exception(SWIG_RuntimeError, msg.c_str());
-    }
-  }
-  resultobj = SWIG_From_double(static_cast< double >(result));
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
 SWIGINTERN PyObject *_wrap_AutocorrelationModel_spectralFunction(PyObject *self, PyObject *args) {
   PyObject *resultobj = 0;
   AutocorrelationModel *arg1 = (AutocorrelationModel *) 0 ;
@@ -53571,40 +53503,6 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_LinearGrowthModel_rms(PyObject *self, PyObject *args) {
-  PyObject *resultobj = 0;
-  LinearGrowthModel *arg1 = (LinearGrowthModel *) 0 ;
-  void *argp1 = 0 ;
-  int res1 = 0 ;
-  PyObject *swig_obj[1] ;
-  double result;
-  
-  (void)self;
-  if (!args) SWIG_fail;
-  swig_obj[0] = args;
-  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_LinearGrowthModel, 0 |  0 );
-  if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "LinearGrowthModel_rms" "', argument " "1"" of type '" "LinearGrowthModel const *""'"); 
-  }
-  arg1 = reinterpret_cast< LinearGrowthModel * >(argp1);
-  {
-    try {
-      result = (double)((LinearGrowthModel const *)arg1)->rms();
-    } catch (const std::exception& ex) {
-      // message shown in the Python interpreter
-      const std::string msg {
-        "BornAgain C++ Exception: " + std::string(ex.what())
-      };
-      SWIG_exception(SWIG_RuntimeError, msg.c_str());
-    }
-  }
-  resultobj = SWIG_From_double(static_cast< double >(result));
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
 SWIGINTERN PyObject *_wrap_LinearGrowthModel_crosscorrSpectrum(PyObject *self, PyObject *args) {
   PyObject *resultobj = 0;
   LinearGrowthModel *arg1 = (LinearGrowthModel *) 0 ;
@@ -75468,7 +75366,6 @@ static PyMethodDef SwigMethods[] = {
 	 { "LayerRoughness_clone", _wrap_LayerRoughness_clone, METH_O, "LayerRoughness_clone(LayerRoughness self) -> LayerRoughness"},
 	 { "LayerRoughness_className", _wrap_LayerRoughness_className, METH_O, "LayerRoughness_className(LayerRoughness self) -> std::string"},
 	 { "LayerRoughness_nodeChildren", _wrap_LayerRoughness_nodeChildren, METH_O, "LayerRoughness_nodeChildren(LayerRoughness self) -> swig_dummy_type_const_inode_vector"},
-	 { "LayerRoughness_rms", _wrap_LayerRoughness_rms, METH_O, "LayerRoughness_rms(LayerRoughness self) -> double"},
 	 { "LayerRoughness_autocorrelationModel", _wrap_LayerRoughness_autocorrelationModel, METH_O, "LayerRoughness_autocorrelationModel(LayerRoughness self) -> AutocorrelationModel"},
 	 { "LayerRoughness_interlayerModel", _wrap_LayerRoughness_interlayerModel, METH_O, "LayerRoughness_interlayerModel(LayerRoughness self) -> InterlayerModel"},
 	 { "LayerRoughness_crosscorrelationModel", _wrap_LayerRoughness_crosscorrelationModel, METH_O, "LayerRoughness_crosscorrelationModel(LayerRoughness self) -> CrosscorrelationModel"},
@@ -75503,7 +75400,6 @@ static PyMethodDef SwigMethods[] = {
 	 { "MultiLayer_swigregister", MultiLayer_swigregister, METH_O, NULL},
 	 { "MultiLayer_swiginit", MultiLayer_swiginit, METH_VARARGS, NULL},
 	 { "AutocorrelationModel_clone", _wrap_AutocorrelationModel_clone, METH_O, "AutocorrelationModel_clone(AutocorrelationModel self) -> AutocorrelationModel"},
-	 { "AutocorrelationModel_rms", _wrap_AutocorrelationModel_rms, METH_O, "AutocorrelationModel_rms(AutocorrelationModel self) -> double"},
 	 { "AutocorrelationModel_spectralFunction", _wrap_AutocorrelationModel_spectralFunction, METH_VARARGS, "AutocorrelationModel_spectralFunction(AutocorrelationModel self, double spatial_f) -> double"},
 	 { "AutocorrelationModel_setMaxSpatialFrequency", _wrap_AutocorrelationModel_setMaxSpatialFrequency, METH_VARARGS, "AutocorrelationModel_setMaxSpatialFrequency(AutocorrelationModel self, double val)"},
 	 { "AutocorrelationModel_maxSpatialFrequency", _wrap_AutocorrelationModel_maxSpatialFrequency, METH_O, "AutocorrelationModel_maxSpatialFrequency(AutocorrelationModel self) -> double"},
@@ -75535,7 +75431,6 @@ static PyMethodDef SwigMethods[] = {
 		"LinearGrowthModel_spectralFunction(LinearGrowthModel self, double spectrum_below, double thickness, double spatial_f) -> double\n"
 		"LinearGrowthModel_spectralFunction(LinearGrowthModel self, double spatial_f) -> double\n"
 		""},
-	 { "LinearGrowthModel_rms", _wrap_LinearGrowthModel_rms, METH_O, "LinearGrowthModel_rms(LinearGrowthModel self) -> double"},
 	 { "LinearGrowthModel_crosscorrSpectrum", _wrap_LinearGrowthModel_crosscorrSpectrum, METH_VARARGS, "LinearGrowthModel_crosscorrSpectrum(LinearGrowthModel self, double spectrum_below, double thickness, double spatial_f) -> double"},
 	 { "LinearGrowthModel_setClusterVolume", _wrap_LinearGrowthModel_setClusterVolume, METH_VARARGS, "LinearGrowthModel_setClusterVolume(LinearGrowthModel self, double particle_volume)"},
 	 { "LinearGrowthModel_clusterVolume", _wrap_LinearGrowthModel_clusterVolume, METH_O, "LinearGrowthModel_clusterVolume(LinearGrowthModel self) -> double"},
-- 
GitLab