diff --git a/Core/Export/SampleToPython.cpp b/Core/Export/SampleToPython.cpp
index 607b37c131b483d5579127858b2750e57f6ae0a6..9ecf2397e546d4a0fea6f572125a52d746eee6ec 100644
--- a/Core/Export/SampleToPython.cpp
+++ b/Core/Export/SampleToPython.cpp
@@ -556,17 +556,17 @@ std::string SampleToPython::defineMultiLayers() const
         if (numberOfLayers) {
             result << indent() << key << ".addLayer(" << m_objs->obj2key(s->layer(0)) << ")\n";
 
-            size_t iLayer = 1;
-            while (iLayer != numberOfLayers) {
-                const LayerInterface* layerInterface = s->layerInterface(iLayer - 1);
+            size_t i_layer = 1;
+            while (i_layer != numberOfLayers) {
+                const LayerInterface* layerInterface = s->layerInterface(i_layer - 1);
                 if (const LayerRoughness* rough = layerInterface->getRoughness())
                     result << indent() << key << ".addLayerWithTopRoughness("
-                           << m_objs->obj2key(s->layer(iLayer)) << ", "
-                           << m_objs->obj2key(rough) << ")\n";
+                           << m_objs->obj2key(s->layer(i_layer)) << ", " << m_objs->obj2key(rough)
+                           << ")\n";
                 else
-                    result << indent() << key << ".addLayer("
-                           << m_objs->obj2key(s->layer(iLayer)) << ")\n";
-                iLayer++;
+                    result << indent() << key << ".addLayer(" << m_objs->obj2key(s->layer(i_layer))
+                           << ")\n";
+                i_layer++;
             }
         }
         result << "\n" << indent() << "return " << key << "\n";
diff --git a/Device/Coord/Axes.h b/Device/Coord/Axes.h
index 891dd25f3c694ce900ac3b40290cf2bfab66b120..1491f809541cbcb4c1217dbe32d34b81a8d34624 100644
--- a/Device/Coord/Axes.h
+++ b/Device/Coord/Axes.h
@@ -12,8 +12,8 @@
 //
 //  ************************************************************************************************
 
-#ifndef BORNAGAIN_DEVICE_UNIT_AXES_H
-#define BORNAGAIN_DEVICE_UNIT_AXES_H
+#ifndef BORNAGAIN_DEVICE_COORD_AXES_H
+#define BORNAGAIN_DEVICE_COORD_AXES_H
 
 //! Wrapper around enum, required for clean Swig interface.
 //! Swig will expose Axes.UNDEFINED etc, i.e. "Coords" will be omitted.
@@ -24,4 +24,4 @@ public:
     enum Coords { UNDEFINED, NBINS, RADIANS, DEGREES, MM, QSPACE, QXQY, RQ4 };
 };
 
-#endif // BORNAGAIN_DEVICE_UNIT_AXES_H
+#endif // BORNAGAIN_DEVICE_COORD_AXES_H
diff --git a/Device/Coord/AxisNames.h b/Device/Coord/AxisNames.h
index 960244f8ba40fee0f7957bcb77f92d87f1ea7dbd..94d2f3663e9f44317e12e1be7920ec4ef136d480 100644
--- a/Device/Coord/AxisNames.h
+++ b/Device/Coord/AxisNames.h
@@ -17,8 +17,8 @@
 #endif
 
 #ifndef USER_API
-#ifndef BORNAGAIN_DEVICE_UNIT_AXISNAMES_H
-#define BORNAGAIN_DEVICE_UNIT_AXISNAMES_H
+#ifndef BORNAGAIN_DEVICE_COORD_AXISNAMES_H
+#define BORNAGAIN_DEVICE_COORD_AXISNAMES_H
 
 #include "Device/Coord/Axes.h"
 #include "Wrap/WinDllMacros.h"
@@ -49,5 +49,5 @@ extern BA_DEVICE_API_ const std::map<Axes::Coords, std::string> specAxisQ;
 extern BA_DEVICE_API_ const std::map<Axes::Coords, std::string> sampleDepthAxis;
 } // namespace DataUtils::AxisNames
 
-#endif // BORNAGAIN_DEVICE_UNIT_AXISNAMES_H
+#endif // BORNAGAIN_DEVICE_COORD_AXISNAMES_H
 #endif // USER_API
diff --git a/Device/Coord/CoordSystem1D.h b/Device/Coord/CoordSystem1D.h
index 4b29bb571e0e6f9f1d6e79d2220221ef427898f8..2a5820b789f7bebfae349ebc2bb79b1d30ecea11 100644
--- a/Device/Coord/CoordSystem1D.h
+++ b/Device/Coord/CoordSystem1D.h
@@ -17,8 +17,8 @@
 #endif
 
 #ifndef USER_API
-#ifndef BORNAGAIN_CORE_SCAN_COORDSYSTEM1D_H
-#define BORNAGAIN_CORE_SCAN_COORDSYSTEM1D_H
+#ifndef BORNAGAIN_DEVICE_COORD_COORDSYSTEM1D_H
+#define BORNAGAIN_DEVICE_COORD_COORDSYSTEM1D_H
 
 #include "Device/Coord/ICoordSystem.h"
 
@@ -121,5 +121,5 @@ private:
     std::function<double(double)> getTraslatorTo(Axes::Coords units) const override;
 };
 
-#endif // BORNAGAIN_CORE_SCAN_COORDSYSTEM1D_H
+#endif // BORNAGAIN_DEVICE_COORD_COORDSYSTEM1D_H
 #endif // USER_API
diff --git a/Device/Coord/ICoordSystem.h b/Device/Coord/ICoordSystem.h
index 0b553502d29e7a20a279ea89b30460f7d65407f6..0e7551a0f864c3c91d055968ec9a8b1d18ca7229 100644
--- a/Device/Coord/ICoordSystem.h
+++ b/Device/Coord/ICoordSystem.h
@@ -17,8 +17,8 @@
 #endif
 
 #ifndef USER_API
-#ifndef BORNAGAIN_DEVICE_UNIT_ICOORDSYSTEM_H
-#define BORNAGAIN_DEVICE_UNIT_ICOORDSYSTEM_H
+#ifndef BORNAGAIN_DEVICE_COORD_ICOORDSYSTEM_H
+#define BORNAGAIN_DEVICE_COORD_ICOORDSYSTEM_H
 
 #include "Base/Types/ICloneable.h"
 #include "Device/Coord/Axes.h"
@@ -79,5 +79,5 @@ private:
     virtual std::vector<std::map<Axes::Coords, std::string>> createNameMaps() const = 0;
 };
 
-#endif // BORNAGAIN_DEVICE_UNIT_ICOORDSYSTEM_H
+#endif // BORNAGAIN_DEVICE_COORD_ICOORDSYSTEM_H
 #endif // USER_API
diff --git a/GUI/Application/ApplicationSettings.cpp b/GUI/Application/ApplicationSettings.cpp
index 0742f96864a8fb41a36ca8a3a72977709e8447de..ee23dc861d9bc997f61f45c09c1e1eba5456dfad 100644
--- a/GUI/Application/ApplicationSettings.cpp
+++ b/GUI/Application/ApplicationSettings.cpp
@@ -2,7 +2,7 @@
 //
 //  BornAgain: simulate and fit reflection and scattering
 //
-//! @file      GUI/mainwindow/projectmanager.cpp
+//! @file      GUI/Application/ApplicationSettings.cpp
 //! @brief     Implements class ProjectManager
 //!
 //! @homepage  http://www.bornagainproject.org
diff --git a/GUI/Models/GUIDomainSampleVisitor.cpp b/GUI/Models/GUIDomainSampleVisitor.cpp
index 9914eeb5dbddb2228afffe6aa87ae77ccae32df9..4343c3c0c0135a42c6829a3a5c691daa46c4721b 100644
--- a/GUI/Models/GUIDomainSampleVisitor.cpp
+++ b/GUI/Models/GUIDomainSampleVisitor.cpp
@@ -107,9 +107,9 @@ void GUIDomainSampleVisitor::visit(const Layer* sample)
 
     const auto* multilayer = dynamic_cast<const MultiLayer*>(m_itemToSample[parent]);
     ASSERT(multilayer);
-    size_t iLayer = SampleUtils::Multilayer::IndexOfLayer(*multilayer, sample);
+    size_t i_layer = SampleUtils::Multilayer::IndexOfLayer(*multilayer, sample);
     const LayerInterface* top_interface =
-        iLayer == 0 ? nullptr : multilayer->layerInterface(iLayer - 1);
+        i_layer == 0 ? nullptr : multilayer->layerInterface(i_layer - 1);
 
     LayerItem* layer_item = m_sampleModel->insertItem<LayerItem>(parent);
     layer_item->setMaterial(createMaterialFromDomain(sample->material()));
diff --git a/GUI/Models/GUIObjectBuilder.cpp b/GUI/Models/GUIObjectBuilder.cpp
index 4aeee32117343a9b3a4c99c905539cadd010ed5a..3044380c246977102fde0a633b5eca485934ecd8 100644
--- a/GUI/Models/GUIObjectBuilder.cpp
+++ b/GUI/Models/GUIObjectBuilder.cpp
@@ -2,7 +2,7 @@
 //
 //  BornAgain: simulate and fit reflection and scattering
 //
-//! @file      GUI/Models/GUI::Model::ObjectBuilder.cpp
+//! @file      GUI/Models/GUIObjectBuilder.cpp
 //! @brief     Implements GUI::Model::ObjectBuilder namespace
 //!
 //! @homepage  http://www.bornagainproject.org
diff --git a/GUI/Views/CommonWidgets/StyleUtils.cpp b/GUI/Views/CommonWidgets/StyleUtils.cpp
index 552c12d23144f62a8a82527c822def163407a5ea..1f8060527d6b1850672e7a564ea66155ce0f2628 100644
--- a/GUI/Views/CommonWidgets/StyleUtils.cpp
+++ b/GUI/Views/CommonWidgets/StyleUtils.cpp
@@ -2,7 +2,7 @@
 //
 //  BornAgain: simulate and fit reflection and scattering
 //
-//! @file      GUI/Views/CommonWidgets/GUI::StyleUtils.cpp
+//! @file      GUI/Views/CommonWidgets/StyleUtils.cpp
 //! @brief     Defines GUI::StyleUtils namespace
 //!
 //! @homepage  http://www.bornagainproject.org
diff --git a/GUI/Views/MaterialEditor/MaterialEditorUtils.cpp b/GUI/Views/MaterialEditor/MaterialEditorUtils.cpp
index 91a8e5053ad19b946e2c69a4d8165e764d206f69..1992b14fbee928bae8d4887a9cfc1f76edd4ae8b 100644
--- a/GUI/Views/MaterialEditor/MaterialEditorUtils.cpp
+++ b/GUI/Views/MaterialEditor/MaterialEditorUtils.cpp
@@ -2,7 +2,7 @@
 //
 //  BornAgain: simulate and fit reflection and scattering
 //
-//! @file      GUI/Views/MaterialEditor/MaterialItemUtils.cpp
+//! @file      GUI/Views/MaterialEditor/MaterialEditorUtils.cpp
 //! @brief     Implements class MaterialItemUtils
 //!
 //! @homepage  http://www.bornagainproject.org
diff --git a/GUI/Views/MaterialEditor/MaterialEditorUtils.h b/GUI/Views/MaterialEditor/MaterialEditorUtils.h
index 4de9eee0da42a07ec4813dfa5d8af4d690850c5f..d8aba65bd1d6aee6d4e1ce60b819b305a690379e 100644
--- a/GUI/Views/MaterialEditor/MaterialEditorUtils.h
+++ b/GUI/Views/MaterialEditor/MaterialEditorUtils.h
@@ -2,7 +2,7 @@
 //
 //  BornAgain: simulate and fit reflection and scattering
 //
-//! @file      GUI/Views/MaterialEditor/MaterialItemUtils.h
+//! @file      GUI/Views/MaterialEditor/MaterialEditorUtils.h
 //! @brief     Defines namespace GUI::View::MaterialItemUtils
 //!
 //! @homepage  http://www.bornagainproject.org
diff --git a/Resample/Coherence/FFTerm.cpp b/Resample/Coherence/FFTerm.cpp
index 540490ba58aa3c4f2222f92d01a9c6f1c9102852..8ad08d91dbccd8c8f97bd6a8a4f9e3c8c12f8b1e 100644
--- a/Resample/Coherence/FFTerm.cpp
+++ b/Resample/Coherence/FFTerm.cpp
@@ -13,12 +13,11 @@
 //  ************************************************************************************************
 
 #include "Resample/Coherence/FFTerm.h"
-#include "Resample/FFCompute/IComputeFF.h"
+#include "Base/Utils/Assert.h"
+#include "Resample/FFCompute/IComputeScalar.h"
+#include "Resample/FFCompute/IComputePol.h"
 
-CoherentFFTerm::CoherentFFTerm(IComputeFF* computer)
-    : m_computer(computer)
-{
-}
+CoherentFFTerm::CoherentFFTerm(IComputeFF* computer) : m_computer(computer) {}
 
 CoherentFFTerm::CoherentFFTerm(const CoherentFFTerm& other)
     : m_computer(std::unique_ptr<IComputeFF>(other.m_computer->clone()))
@@ -29,12 +28,16 @@ CoherentFFTerm::~CoherentFFTerm() = default;
 
 complex_t CoherentFFTerm::coherentFF(const SimulationElement& ele) const
 {
-    return m_computer->coherentFF(ele);
+    const auto* computer = dynamic_cast<const IComputeScalar*>(m_computer.get());
+    ASSERT(computer);
+    return computer->coherentFF(ele);
 }
 
 Eigen::Matrix2cd CoherentFFTerm::coherentPolFF(const SimulationElement& ele) const
 {
-    return m_computer->coherentPolFF(ele);
+    const auto* computer = dynamic_cast<const IComputePol*>(m_computer.get());
+    ASSERT(computer);
+    return computer->coherentPolFF(ele);
 }
 
 double CoherentFFTerm::radialExtension() const
diff --git a/Resample/Coherence/FFTerm.h b/Resample/Coherence/FFTerm.h
index 7d6e3770eea4691ac82d58fff991589810f9f418..cd90ee8837a9d321a50364f2a71c6a0f49262fd5 100644
--- a/Resample/Coherence/FFTerm.h
+++ b/Resample/Coherence/FFTerm.h
@@ -43,7 +43,7 @@ public:
     double radialExtension() const;
 
 private:
-    const std::unique_ptr<IComputeFF> m_computer;
+    const std::unique_ptr<const IComputeFF> m_computer;
 };
 
 #endif // BORNAGAIN_RESAMPLE_COHERENCE_FFTERM_H
diff --git a/Resample/FFCompute/ComputeBA.cpp b/Resample/FFCompute/ComputeBA.cpp
index 028b556ebddbd4707e1f54eac922eec8f7c94e2d..86b19402b8272ae3b3e83cb4e6a54733a9e1991d 100644
--- a/Resample/FFCompute/ComputeBA.cpp
+++ b/Resample/FFCompute/ComputeBA.cpp
@@ -16,16 +16,18 @@
 #include "Sample/Material/WavevectorInfo.h"
 #include "Sample/Scattering/IFormFactor.h"
 
-ComputeBA::ComputeBA(const IFormFactor& ff, size_t iLayer) : IComputeFF(ff, iLayer) {}
+ComputeBA::ComputeBA(const IFormFactor& ff, size_t i_layer) : IComputeScalar(ff, i_layer) {}
 
 ComputeBA::~ComputeBA() = default;
 
 ComputeBA* ComputeBA::clone() const
 {
-    return new ComputeBA(*m_ff, m_iLayer);
+    return new ComputeBA(*m_ff, m_i_layer);
 }
 
-complex_t ComputeBA::computeFF(const WavevectorInfo& wavevectors) const
+complex_t ComputeBA::computeFF(const WavevectorInfo& wavevectors,
+                        const std::unique_ptr<const IFlux>&,
+                        const std::unique_ptr<const IFlux>&) const
 {
     return m_ff->theFF(wavevectors);
 }
diff --git a/Resample/FFCompute/ComputeBA.h b/Resample/FFCompute/ComputeBA.h
index 7f8ccc572f3e907ad6a2a31c48e290f0da66851c..ed5e2d4da938671c43ecdf76a9aeee569a8efd53 100644
--- a/Resample/FFCompute/ComputeBA.h
+++ b/Resample/FFCompute/ComputeBA.h
@@ -20,22 +20,24 @@
 #ifndef BORNAGAIN_RESAMPLE_FFCOMPUTE_COMPUTEBA_H
 #define BORNAGAIN_RESAMPLE_FFCOMPUTE_COMPUTEBA_H
 
-#include "Resample/FFCompute/IComputeFF.h"
+#include "Resample/FFCompute/IComputeScalar.h"
 #include <memory>
 
 //! Provides scalar form factor evaluation in Born Approximation for given IFormFactor.
 
 //! @ingroup formfactors_internal
 
-class ComputeBA : public IComputeFF {
+class ComputeBA : public IComputeScalar {
 public:
-    ComputeBA(const IFormFactor& ff, size_t iLayer);
+    ComputeBA(const IFormFactor& ff, size_t i_layer);
     ~ComputeBA() override;
 
     ComputeBA* clone() const override;
 
     //! Calculates and returns a form factor calculation in BA
-    complex_t computeFF(const WavevectorInfo& wavevectors) const override;
+    complex_t computeFF(const WavevectorInfo& wavevectors,
+                        const std::unique_ptr<const IFlux>& inFlux,
+                        const std::unique_ptr<const IFlux>& outFlux) const override;
 };
 
 #endif // BORNAGAIN_RESAMPLE_FFCOMPUTE_COMPUTEBA_H
diff --git a/Resample/FFCompute/ComputeBAPol.cpp b/Resample/FFCompute/ComputeBAPol.cpp
index 13ccfa43e250bf672d1058a8022bbb26581327c2..eeb87dfb666f897eb85541ef64cf19474fe8ed9d 100644
--- a/Resample/FFCompute/ComputeBAPol.cpp
+++ b/Resample/FFCompute/ComputeBAPol.cpp
@@ -14,24 +14,20 @@
 
 #include "Resample/FFCompute/ComputeBAPol.h"
 #include "Sample/Material/WavevectorInfo.h"
-#include <stdexcept>
+#include "Sample/Scattering/IFormFactor.h"
 
-ComputeBAPol::ComputeBAPol(const IFormFactor& ff, size_t iLayer) : IComputeFF(ff, iLayer) {}
+ComputeBAPol::ComputeBAPol(const IFormFactor& ff, size_t i_layer) : IComputePol(ff, i_layer) {}
 
 ComputeBAPol::~ComputeBAPol() = default;
 
 ComputeBAPol* ComputeBAPol::clone() const
 {
-    return new ComputeBAPol(*m_ff, m_iLayer);
+    return new ComputeBAPol(*m_ff, m_i_layer);
 }
 
-complex_t ComputeBAPol::computeFF(const WavevectorInfo&) const
-{
-    throw std::runtime_error("ComputeBAPol::evaluate: "
-                             "should never be called for matrix interactions");
-}
-
-Eigen::Matrix2cd ComputeBAPol::computePolFF(const WavevectorInfo& wavevectors) const
+Eigen::Matrix2cd ComputeBAPol::computePolFF(const WavevectorInfo& wavevectors,
+                        const std::unique_ptr<const IFlux>&,
+                        const std::unique_ptr<const IFlux>&) const
 {
     Eigen::Matrix2cd ff_BA = m_ff->thePolFF(wavevectors);
     Eigen::Matrix2cd result;
diff --git a/Resample/FFCompute/ComputeBAPol.h b/Resample/FFCompute/ComputeBAPol.h
index 13b27e643060b953c837cd935dca2a1fee099dc8..9411f825dc54bf1cca3ed7e6c9d1a30877279f01 100644
--- a/Resample/FFCompute/ComputeBAPol.h
+++ b/Resample/FFCompute/ComputeBAPol.h
@@ -20,26 +20,27 @@
 #ifndef BORNAGAIN_RESAMPLE_FFCOMPUTE_COMPUTEBAPOL_H
 #define BORNAGAIN_RESAMPLE_FFCOMPUTE_COMPUTEBAPOL_H
 
-#include "Resample/FFCompute/IComputeFF.h"
-#include "Sample/Scattering/IFormFactor.h"
+#include "Resample/FFCompute/IComputePol.h"
 #include <memory>
 
+class SimulationElement;
+class WavevectorInfo;
+
 //! Provides polarized form factor evaluation in Born Approximation for given IFormFactor.
 
 //! @ingroup formfactors_internal
 
-class ComputeBAPol : public IComputeFF {
+class ComputeBAPol : public IComputePol {
 public:
-    ComputeBAPol(const IFormFactor& ff, size_t iLayer);
+    ComputeBAPol(const IFormFactor& ff, size_t i_layer);
     ~ComputeBAPol() override;
 
     ComputeBAPol* clone() const override;
 
-    //! Throws not-implemented exception
-    complex_t computeFF(const WavevectorInfo& wavevectors) const override;
-
     //! Calculates and returns a polarized form factor calculation in BA
-    Eigen::Matrix2cd computePolFF(const WavevectorInfo& wavevectors) const override;
+    Eigen::Matrix2cd computePolFF(const WavevectorInfo& wavevectors,
+                        const std::unique_ptr<const IFlux>&,
+                        const std::unique_ptr<const IFlux>&) const override;
 };
 
 #endif // BORNAGAIN_RESAMPLE_FFCOMPUTE_COMPUTEBAPOL_H
diff --git a/Resample/FFCompute/ComputeDWBA.cpp b/Resample/FFCompute/ComputeDWBA.cpp
index 83381a5a9801bad2dedeea6e0a4db551d4fb83ff..fb22dbc47550c85d61f96a0b8ed276be73d4852d 100644
--- a/Resample/FFCompute/ComputeDWBA.cpp
+++ b/Resample/FFCompute/ComputeDWBA.cpp
@@ -17,30 +17,28 @@
 #include "Sample/Material/WavevectorInfo.h"
 #include "Sample/Scattering/IFormFactor.h"
 
-ComputeDWBA::ComputeDWBA(const IFormFactor& ff, size_t iLayer) : IComputeFF(ff, iLayer) {}
+ComputeDWBA::ComputeDWBA(const IFormFactor& ff, size_t i_layer) : IComputeScalar(ff, i_layer) {}
 
 ComputeDWBA::~ComputeDWBA() = default;
 
 ComputeDWBA* ComputeDWBA::clone() const
 {
-    ComputeDWBA* result = new ComputeDWBA(*m_ff, m_iLayer);
-    auto in_coefs = std::unique_ptr<const IFlux>(m_inFlux ? m_inFlux->clone() : nullptr);
-    auto out_coefs = std::unique_ptr<const IFlux>(m_outFlux ? m_outFlux->clone() : nullptr);
-    result->setFlux(std::move(in_coefs), std::move(out_coefs));
-    return result;
+    return new ComputeDWBA(*m_ff, m_i_layer);
 }
 
-complex_t ComputeDWBA::computeFF(const WavevectorInfo& wavevectors) const
+complex_t ComputeDWBA::computeFF(const WavevectorInfo& wavevectors,
+                        const std::unique_ptr<const IFlux>& inFlux,
+                        const std::unique_ptr<const IFlux>& outFlux) const
 {
     // Retrieve the two different incoming wavevectors in the layer
     cvector_t k_i_T = wavevectors.getKi();
-    k_i_T.setZ(-m_inFlux->getScalarKz());
+    k_i_T.setZ(-inFlux->getScalarKz());
     cvector_t k_i_R = k_i_T;
     k_i_R.setZ(-k_i_T.z());
 
     // Retrieve the two different outgoing wavevector bins in the layer
     cvector_t k_f_T = wavevectors.getKf();
-    k_f_T.setZ(m_outFlux->getScalarKz());
+    k_f_T.setZ(outFlux->getScalarKz());
     cvector_t k_f_R = k_f_T;
     k_f_R.setZ(-k_f_T.z());
 
@@ -52,10 +50,10 @@ complex_t ComputeDWBA::computeFF(const WavevectorInfo& wavevectors) const
     WavevectorInfo k_RR(k_i_R, k_f_R, wavelength);
 
     // Get the four R,T coefficients
-    complex_t T_in = m_inFlux->getScalarT();
-    complex_t R_in = m_inFlux->getScalarR();
-    complex_t T_out = m_outFlux->getScalarT();
-    complex_t R_out = m_outFlux->getScalarR();
+    complex_t T_in = inFlux->getScalarT();
+    complex_t R_in = inFlux->getScalarR();
+    complex_t T_out = outFlux->getScalarT();
+    complex_t R_out = outFlux->getScalarR();
 
     // The four different scattering contributions; S stands for scattering
     // off the particle, R for reflection off the layer interface
@@ -66,9 +64,3 @@ complex_t ComputeDWBA::computeFF(const WavevectorInfo& wavevectors) const
 
     return term_S + term_RS + term_SR + term_RSR;
 }
-
-void ComputeDWBA::setFlux(std::unique_ptr<const IFlux> inFlux, std::unique_ptr<const IFlux> outFlux)
-{
-    m_inFlux = std::move(inFlux);
-    m_outFlux = std::move(outFlux);
-}
diff --git a/Resample/FFCompute/ComputeDWBA.h b/Resample/FFCompute/ComputeDWBA.h
index bf6910ffb2e5daeff17c6f299402a9890cea549b..2cf39177ee11a6ee291b782bd679a2a40b362b84 100644
--- a/Resample/FFCompute/ComputeDWBA.h
+++ b/Resample/FFCompute/ComputeDWBA.h
@@ -20,7 +20,7 @@
 #ifndef BORNAGAIN_RESAMPLE_FFCOMPUTE_COMPUTEDWBA_H
 #define BORNAGAIN_RESAMPLE_FFCOMPUTE_COMPUTEDWBA_H
 
-#include "Resample/FFCompute/IComputeFF.h"
+#include "Resample/FFCompute/IComputeScalar.h"
 #include <memory>
 
 class IFlux;
@@ -29,24 +29,19 @@ class IFlux;
 
 //! @ingroup formfactors_internal
 
-class ComputeDWBA : public IComputeFF {
+class ComputeDWBA : public IComputeScalar {
 public:
-    ComputeDWBA(const IFormFactor& ff, size_t iLayer);
+    ComputeDWBA(const IFormFactor& ff, size_t i_layer);
     ~ComputeDWBA() override;
 
     ComputeDWBA* clone() const override;
 
     //! Returns the coherent sum of the four DWBA terms for scalar scattering.
-    complex_t computeFF(const WavevectorInfo& wavevectors) const override;
-
-    void setFlux(std::unique_ptr<const IFlux> inFlux,
-                 std::unique_ptr<const IFlux> outFlux) override;
+    complex_t computeFF(const WavevectorInfo& wavevectors,
+                        const std::unique_ptr<const IFlux>& inFlux,
+                        const std::unique_ptr<const IFlux>& outFlux) const override;
 
     friend class TestPolarizedDWBATerms;
-
-private:
-    std::unique_ptr<const IFlux> m_inFlux;
-    std::unique_ptr<const IFlux> m_outFlux;
 };
 
 #endif // BORNAGAIN_RESAMPLE_FFCOMPUTE_COMPUTEDWBA_H
diff --git a/Resample/FFCompute/ComputeDWBAPol.cpp b/Resample/FFCompute/ComputeDWBAPol.cpp
index 9de9bc273d97f4d0ae318ab0a43fb4edac2258f4..1cd7a983cfd9a642cab17a993a5c700375857c46 100644
--- a/Resample/FFCompute/ComputeDWBAPol.cpp
+++ b/Resample/FFCompute/ComputeDWBAPol.cpp
@@ -28,45 +28,36 @@ complex_t VecMatVecProduct(const Eigen::Vector2cd& vec1, const Eigen::Matrix2cd&
 } // namespace
 
 
-ComputeDWBAPol::ComputeDWBAPol(const IFormFactor& ff, size_t iLayer) : IComputeFF(ff, iLayer) {}
+ComputeDWBAPol::ComputeDWBAPol(const IFormFactor& ff, size_t i_layer) : IComputePol(ff, i_layer) {}
 
 ComputeDWBAPol::~ComputeDWBAPol() = default;
 
 ComputeDWBAPol* ComputeDWBAPol::clone() const
 {
-    ComputeDWBAPol* result = new ComputeDWBAPol(*m_ff, m_iLayer);
-    std::unique_ptr<const IFlux> in_coefs =
-        m_inFlux ? std::unique_ptr<const IFlux>(m_inFlux->clone()) : nullptr;
-    std::unique_ptr<const IFlux> out_coefs =
-        m_outFlux ? std::unique_ptr<const IFlux>(m_outFlux->clone()) : nullptr;
-    result->setFlux(std::move(in_coefs), std::move(out_coefs));
-    return result;
+    return new ComputeDWBAPol(*m_ff, m_i_layer);
 }
 
-complex_t ComputeDWBAPol::computeFF(const WavevectorInfo&) const
-{
-    throw std::runtime_error("Bug: forbidden call of ComputeDWBAPol::evaluate");
-}
-
-Eigen::Matrix2cd ComputeDWBAPol::computePolFF(const WavevectorInfo& wavevectors) const
+Eigen::Matrix2cd ComputeDWBAPol::computePolFF(const WavevectorInfo& wavevectors,
+                                const std::unique_ptr<const IFlux>& inFlux,
+                                const std::unique_ptr<const IFlux>& outFlux) const
 {
     // the required wavevectors inside the layer for
     // different eigenmodes and in- and outgoing wavevector;
     complex_t kix = wavevectors.getKi().x();
     complex_t kiy = wavevectors.getKi().y();
-    cvector_t ki_1R(kix, kiy, m_inFlux->getKz()(0));
-    cvector_t ki_1T(kix, kiy, -m_inFlux->getKz()(0));
-    cvector_t ki_2R(kix, kiy, m_inFlux->getKz()(1));
-    cvector_t ki_2T(kix, kiy, -m_inFlux->getKz()(1));
+    cvector_t ki_1R(kix, kiy, inFlux->getKz()(0));
+    cvector_t ki_1T(kix, kiy, -inFlux->getKz()(0));
+    cvector_t ki_2R(kix, kiy, inFlux->getKz()(1));
+    cvector_t ki_2T(kix, kiy, -inFlux->getKz()(1));
 
     cvector_t kf_1R = wavevectors.getKf();
-    kf_1R.setZ(-(complex_t)m_outFlux->getKz()(0));
+    kf_1R.setZ(-(complex_t)outFlux->getKz()(0));
     cvector_t kf_1T = kf_1R;
-    kf_1T.setZ((complex_t)m_outFlux->getKz()(0));
+    kf_1T.setZ((complex_t)outFlux->getKz()(0));
     cvector_t kf_2R = kf_1R;
-    kf_2R.setZ(-(complex_t)m_outFlux->getKz()(1));
+    kf_2R.setZ(-(complex_t)outFlux->getKz()(1));
     cvector_t kf_2T = kf_1R;
-    kf_2T.setZ((complex_t)m_outFlux->getKz()(1));
+    kf_2T.setZ((complex_t)outFlux->getKz()(1));
     // now each of the 16 matrix terms of the polarized DWBA is calculated:
     // NOTE: when the underlying reflection/transmission coefficients are
     // scalar, the eigenmodes have identical eigenvalues and spin polarization
@@ -91,111 +82,104 @@ Eigen::Matrix2cd ComputeDWBAPol::computePolFF(const WavevectorInfo& wavevectors)
 
     // eigenmode 1 -> eigenmode 1: direct scattering
     ff_BA = m_ff->thePolFF({ki_1T, kf_1T, wavelength});
-    M11_S(0, 0) = -VecMatVecProduct(m_outFlux->T1min(), ff_BA, m_inFlux->T1plus());
-    M11_S(0, 1) = VecMatVecProduct(m_outFlux->T1plus(), ff_BA, m_inFlux->T1plus());
-    M11_S(1, 0) = -VecMatVecProduct(m_outFlux->T1min(), ff_BA, m_inFlux->T1min());
-    M11_S(1, 1) = VecMatVecProduct(m_outFlux->T1plus(), ff_BA, m_inFlux->T1min());
+    M11_S(0, 0) = -VecMatVecProduct(outFlux->T1min(), ff_BA, inFlux->T1plus());
+    M11_S(0, 1) = VecMatVecProduct(outFlux->T1plus(), ff_BA, inFlux->T1plus());
+    M11_S(1, 0) = -VecMatVecProduct(outFlux->T1min(), ff_BA, inFlux->T1min());
+    M11_S(1, 1) = VecMatVecProduct(outFlux->T1plus(), ff_BA, inFlux->T1min());
     // eigenmode 1 -> eigenmode 1: reflection and then scattering
     ff_BA = m_ff->thePolFF({ki_1R, kf_1T, wavelength});
-    M11_RS(0, 0) = -VecMatVecProduct(m_outFlux->T1min(), ff_BA, m_inFlux->R1plus());
-    M11_RS(0, 1) = VecMatVecProduct(m_outFlux->T1plus(), ff_BA, m_inFlux->R1plus());
-    M11_RS(1, 0) = -VecMatVecProduct(m_outFlux->T1min(), ff_BA, m_inFlux->R1min());
-    M11_RS(1, 1) = VecMatVecProduct(m_outFlux->T1plus(), ff_BA, m_inFlux->R1min());
+    M11_RS(0, 0) = -VecMatVecProduct(outFlux->T1min(), ff_BA, inFlux->R1plus());
+    M11_RS(0, 1) = VecMatVecProduct(outFlux->T1plus(), ff_BA, inFlux->R1plus());
+    M11_RS(1, 0) = -VecMatVecProduct(outFlux->T1min(), ff_BA, inFlux->R1min());
+    M11_RS(1, 1) = VecMatVecProduct(outFlux->T1plus(), ff_BA, inFlux->R1min());
     // eigenmode 1 -> eigenmode 1: scattering and then reflection
     ff_BA = m_ff->thePolFF({ki_1T, kf_1R, wavelength});
-    M11_SR(0, 0) = -VecMatVecProduct(m_outFlux->R1min(), ff_BA, m_inFlux->T1plus());
-    M11_SR(0, 1) = VecMatVecProduct(m_outFlux->R1plus(), ff_BA, m_inFlux->T1plus());
-    M11_SR(1, 0) = -VecMatVecProduct(m_outFlux->R1min(), ff_BA, m_inFlux->T1min());
-    M11_SR(1, 1) = VecMatVecProduct(m_outFlux->R1plus(), ff_BA, m_inFlux->T1min());
+    M11_SR(0, 0) = -VecMatVecProduct(outFlux->R1min(), ff_BA, inFlux->T1plus());
+    M11_SR(0, 1) = VecMatVecProduct(outFlux->R1plus(), ff_BA, inFlux->T1plus());
+    M11_SR(1, 0) = -VecMatVecProduct(outFlux->R1min(), ff_BA, inFlux->T1min());
+    M11_SR(1, 1) = VecMatVecProduct(outFlux->R1plus(), ff_BA, inFlux->T1min());
     // eigenmode 1 -> eigenmode 1: reflection, scattering and again reflection
     ff_BA = m_ff->thePolFF({ki_1R, kf_1R, wavelength});
-    M11_RSR(0, 0) = -VecMatVecProduct(m_outFlux->R1min(), ff_BA, m_inFlux->R1plus());
-    M11_RSR(0, 1) = VecMatVecProduct(m_outFlux->R1plus(), ff_BA, m_inFlux->R1plus());
-    M11_RSR(1, 0) = -VecMatVecProduct(m_outFlux->R1min(), ff_BA, m_inFlux->R1min());
-    M11_RSR(1, 1) = VecMatVecProduct(m_outFlux->R1plus(), ff_BA, m_inFlux->R1min());
+    M11_RSR(0, 0) = -VecMatVecProduct(outFlux->R1min(), ff_BA, inFlux->R1plus());
+    M11_RSR(0, 1) = VecMatVecProduct(outFlux->R1plus(), ff_BA, inFlux->R1plus());
+    M11_RSR(1, 0) = -VecMatVecProduct(outFlux->R1min(), ff_BA, inFlux->R1min());
+    M11_RSR(1, 1) = VecMatVecProduct(outFlux->R1plus(), ff_BA, inFlux->R1min());
 
     // eigenmode 1 -> eigenmode 2: direct scattering
     ff_BA = m_ff->thePolFF({ki_1T, kf_2T, wavelength});
-    M12_S(0, 0) = -VecMatVecProduct(m_outFlux->T2min(), ff_BA, m_inFlux->T1plus());
-    M12_S(0, 1) = VecMatVecProduct(m_outFlux->T2plus(), ff_BA, m_inFlux->T1plus());
-    M12_S(1, 0) = -VecMatVecProduct(m_outFlux->T2min(), ff_BA, m_inFlux->T1min());
-    M12_S(1, 1) = VecMatVecProduct(m_outFlux->T2plus(), ff_BA, m_inFlux->T1min());
+    M12_S(0, 0) = -VecMatVecProduct(outFlux->T2min(), ff_BA, inFlux->T1plus());
+    M12_S(0, 1) = VecMatVecProduct(outFlux->T2plus(), ff_BA, inFlux->T1plus());
+    M12_S(1, 0) = -VecMatVecProduct(outFlux->T2min(), ff_BA, inFlux->T1min());
+    M12_S(1, 1) = VecMatVecProduct(outFlux->T2plus(), ff_BA, inFlux->T1min());
     // eigenmode 1 -> eigenmode 2: reflection and then scattering
     ff_BA = m_ff->thePolFF({ki_1R, kf_2T, wavelength});
-    M12_RS(0, 0) = -VecMatVecProduct(m_outFlux->T2min(), ff_BA, m_inFlux->R1plus());
-    M12_RS(0, 1) = VecMatVecProduct(m_outFlux->T2plus(), ff_BA, m_inFlux->R1plus());
-    M12_RS(1, 0) = -VecMatVecProduct(m_outFlux->T2min(), ff_BA, m_inFlux->R1min());
-    M12_RS(1, 1) = VecMatVecProduct(m_outFlux->T2plus(), ff_BA, m_inFlux->R1min());
+    M12_RS(0, 0) = -VecMatVecProduct(outFlux->T2min(), ff_BA, inFlux->R1plus());
+    M12_RS(0, 1) = VecMatVecProduct(outFlux->T2plus(), ff_BA, inFlux->R1plus());
+    M12_RS(1, 0) = -VecMatVecProduct(outFlux->T2min(), ff_BA, inFlux->R1min());
+    M12_RS(1, 1) = VecMatVecProduct(outFlux->T2plus(), ff_BA, inFlux->R1min());
     // eigenmode 1 -> eigenmode 2: scattering and then reflection
     ff_BA = m_ff->thePolFF({ki_1T, kf_2R, wavelength});
-    M12_SR(0, 0) = -VecMatVecProduct(m_outFlux->R2min(), ff_BA, m_inFlux->T1plus());
-    M12_SR(0, 1) = VecMatVecProduct(m_outFlux->R2plus(), ff_BA, m_inFlux->T1plus());
-    M12_SR(1, 0) = -VecMatVecProduct(m_outFlux->R2min(), ff_BA, m_inFlux->T1min());
-    M12_SR(1, 1) = VecMatVecProduct(m_outFlux->R2plus(), ff_BA, m_inFlux->T1min());
+    M12_SR(0, 0) = -VecMatVecProduct(outFlux->R2min(), ff_BA, inFlux->T1plus());
+    M12_SR(0, 1) = VecMatVecProduct(outFlux->R2plus(), ff_BA, inFlux->T1plus());
+    M12_SR(1, 0) = -VecMatVecProduct(outFlux->R2min(), ff_BA, inFlux->T1min());
+    M12_SR(1, 1) = VecMatVecProduct(outFlux->R2plus(), ff_BA, inFlux->T1min());
     // eigenmode 1 -> eigenmode 2: reflection, scattering and again reflection
     ff_BA = m_ff->thePolFF({ki_1R, kf_2R, wavelength});
-    M12_RSR(0, 0) = -VecMatVecProduct(m_outFlux->R2min(), ff_BA, m_inFlux->R1plus());
-    M12_RSR(0, 1) = VecMatVecProduct(m_outFlux->R2plus(), ff_BA, m_inFlux->R1plus());
-    M12_RSR(1, 0) = -VecMatVecProduct(m_outFlux->R2min(), ff_BA, m_inFlux->R1min());
-    M12_RSR(1, 1) = VecMatVecProduct(m_outFlux->R2plus(), ff_BA, m_inFlux->R1min());
+    M12_RSR(0, 0) = -VecMatVecProduct(outFlux->R2min(), ff_BA, inFlux->R1plus());
+    M12_RSR(0, 1) = VecMatVecProduct(outFlux->R2plus(), ff_BA, inFlux->R1plus());
+    M12_RSR(1, 0) = -VecMatVecProduct(outFlux->R2min(), ff_BA, inFlux->R1min());
+    M12_RSR(1, 1) = VecMatVecProduct(outFlux->R2plus(), ff_BA, inFlux->R1min());
 
     // eigenmode 2 -> eigenmode 1: direct scattering
     ff_BA = m_ff->thePolFF({ki_2T, kf_1T, wavelength});
-    M21_S(0, 0) = -VecMatVecProduct(m_outFlux->T1min(), ff_BA, m_inFlux->T2plus());
-    M21_S(0, 1) = VecMatVecProduct(m_outFlux->T1plus(), ff_BA, m_inFlux->T2plus());
-    M21_S(1, 0) = -VecMatVecProduct(m_outFlux->T1min(), ff_BA, m_inFlux->T2min());
-    M21_S(1, 1) = VecMatVecProduct(m_outFlux->T1plus(), ff_BA, m_inFlux->T2min());
+    M21_S(0, 0) = -VecMatVecProduct(outFlux->T1min(), ff_BA, inFlux->T2plus());
+    M21_S(0, 1) = VecMatVecProduct(outFlux->T1plus(), ff_BA, inFlux->T2plus());
+    M21_S(1, 0) = -VecMatVecProduct(outFlux->T1min(), ff_BA, inFlux->T2min());
+    M21_S(1, 1) = VecMatVecProduct(outFlux->T1plus(), ff_BA, inFlux->T2min());
     // eigenmode 2 -> eigenmode 1: reflection and then scattering
     ff_BA = m_ff->thePolFF({ki_2R, kf_1T, wavelength});
-    M21_RS(0, 0) = -VecMatVecProduct(m_outFlux->T1min(), ff_BA, m_inFlux->R2plus());
-    M21_RS(0, 1) = VecMatVecProduct(m_outFlux->T1plus(), ff_BA, m_inFlux->R2plus());
-    M21_RS(1, 0) = -VecMatVecProduct(m_outFlux->T1min(), ff_BA, m_inFlux->R2min());
-    M21_RS(1, 1) = VecMatVecProduct(m_outFlux->T1plus(), ff_BA, m_inFlux->R2min());
+    M21_RS(0, 0) = -VecMatVecProduct(outFlux->T1min(), ff_BA, inFlux->R2plus());
+    M21_RS(0, 1) = VecMatVecProduct(outFlux->T1plus(), ff_BA, inFlux->R2plus());
+    M21_RS(1, 0) = -VecMatVecProduct(outFlux->T1min(), ff_BA, inFlux->R2min());
+    M21_RS(1, 1) = VecMatVecProduct(outFlux->T1plus(), ff_BA, inFlux->R2min());
     // eigenmode 2 -> eigenmode 1: scattering and then reflection
     ff_BA = m_ff->thePolFF({ki_2T, kf_1R, wavelength});
-    M21_SR(0, 0) = -VecMatVecProduct(m_outFlux->R1min(), ff_BA, m_inFlux->T2plus());
-    M21_SR(0, 1) = VecMatVecProduct(m_outFlux->R1plus(), ff_BA, m_inFlux->T2plus());
-    M21_SR(1, 0) = -VecMatVecProduct(m_outFlux->R1min(), ff_BA, m_inFlux->T2min());
-    M21_SR(1, 1) = VecMatVecProduct(m_outFlux->R1plus(), ff_BA, m_inFlux->T2min());
+    M21_SR(0, 0) = -VecMatVecProduct(outFlux->R1min(), ff_BA, inFlux->T2plus());
+    M21_SR(0, 1) = VecMatVecProduct(outFlux->R1plus(), ff_BA, inFlux->T2plus());
+    M21_SR(1, 0) = -VecMatVecProduct(outFlux->R1min(), ff_BA, inFlux->T2min());
+    M21_SR(1, 1) = VecMatVecProduct(outFlux->R1plus(), ff_BA, inFlux->T2min());
     // eigenmode 2 -> eigenmode 1: reflection, scattering and again reflection
     ff_BA = m_ff->thePolFF({ki_2R, kf_1R, wavelength});
-    M21_RSR(0, 0) = -VecMatVecProduct(m_outFlux->R1min(), ff_BA, m_inFlux->R2plus());
-    M21_RSR(0, 1) = VecMatVecProduct(m_outFlux->R1plus(), ff_BA, m_inFlux->R2plus());
-    M21_RSR(1, 0) = -VecMatVecProduct(m_outFlux->R1min(), ff_BA, m_inFlux->R2min());
-    M21_RSR(1, 1) = VecMatVecProduct(m_outFlux->R1plus(), ff_BA, m_inFlux->R2min());
+    M21_RSR(0, 0) = -VecMatVecProduct(outFlux->R1min(), ff_BA, inFlux->R2plus());
+    M21_RSR(0, 1) = VecMatVecProduct(outFlux->R1plus(), ff_BA, inFlux->R2plus());
+    M21_RSR(1, 0) = -VecMatVecProduct(outFlux->R1min(), ff_BA, inFlux->R2min());
+    M21_RSR(1, 1) = VecMatVecProduct(outFlux->R1plus(), ff_BA, inFlux->R2min());
 
     // eigenmode 2 -> eigenmode 2: direct scattering
     ff_BA = m_ff->thePolFF({ki_2T, kf_2T, wavelength});
-    M22_S(0, 0) = -VecMatVecProduct(m_outFlux->T2min(), ff_BA, m_inFlux->T2plus());
-    M22_S(0, 1) = VecMatVecProduct(m_outFlux->T2plus(), ff_BA, m_inFlux->T2plus());
-    M22_S(1, 0) = -VecMatVecProduct(m_outFlux->T2min(), ff_BA, m_inFlux->T2min());
-    M22_S(1, 1) = VecMatVecProduct(m_outFlux->T2plus(), ff_BA, m_inFlux->T2min());
+    M22_S(0, 0) = -VecMatVecProduct(outFlux->T2min(), ff_BA, inFlux->T2plus());
+    M22_S(0, 1) = VecMatVecProduct(outFlux->T2plus(), ff_BA, inFlux->T2plus());
+    M22_S(1, 0) = -VecMatVecProduct(outFlux->T2min(), ff_BA, inFlux->T2min());
+    M22_S(1, 1) = VecMatVecProduct(outFlux->T2plus(), ff_BA, inFlux->T2min());
     // eigenmode 2 -> eigenmode 2: reflection and then scattering
     ff_BA = m_ff->thePolFF({ki_2R, kf_2T, wavelength});
-    M22_RS(0, 0) = -VecMatVecProduct(m_outFlux->T2min(), ff_BA, m_inFlux->R2plus());
-    M22_RS(0, 1) = VecMatVecProduct(m_outFlux->T2plus(), ff_BA, m_inFlux->R2plus());
-    M22_RS(1, 0) = -VecMatVecProduct(m_outFlux->T2min(), ff_BA, m_inFlux->R2min());
-    M22_RS(1, 1) = VecMatVecProduct(m_outFlux->T2plus(), ff_BA, m_inFlux->R2min());
+    M22_RS(0, 0) = -VecMatVecProduct(outFlux->T2min(), ff_BA, inFlux->R2plus());
+    M22_RS(0, 1) = VecMatVecProduct(outFlux->T2plus(), ff_BA, inFlux->R2plus());
+    M22_RS(1, 0) = -VecMatVecProduct(outFlux->T2min(), ff_BA, inFlux->R2min());
+    M22_RS(1, 1) = VecMatVecProduct(outFlux->T2plus(), ff_BA, inFlux->R2min());
     // eigenmode 2 -> eigenmode 2: scattering and then reflection
     ff_BA = m_ff->thePolFF({ki_2T, kf_2R, wavelength});
-    M22_SR(0, 0) = -VecMatVecProduct(m_outFlux->R2min(), ff_BA, m_inFlux->T2plus());
-    M22_SR(0, 1) = VecMatVecProduct(m_outFlux->R2plus(), ff_BA, m_inFlux->T2plus());
-    M22_SR(1, 0) = -VecMatVecProduct(m_outFlux->R2min(), ff_BA, m_inFlux->T2min());
-    M22_SR(1, 1) = VecMatVecProduct(m_outFlux->R2plus(), ff_BA, m_inFlux->T2min());
+    M22_SR(0, 0) = -VecMatVecProduct(outFlux->R2min(), ff_BA, inFlux->T2plus());
+    M22_SR(0, 1) = VecMatVecProduct(outFlux->R2plus(), ff_BA, inFlux->T2plus());
+    M22_SR(1, 0) = -VecMatVecProduct(outFlux->R2min(), ff_BA, inFlux->T2min());
+    M22_SR(1, 1) = VecMatVecProduct(outFlux->R2plus(), ff_BA, inFlux->T2min());
     // eigenmode 2 -> eigenmode 2: reflection, scattering and again reflection
     ff_BA = m_ff->thePolFF({ki_2R, kf_2R, wavelength});
-    M22_RSR(0, 0) = -VecMatVecProduct(m_outFlux->R2min(), ff_BA, m_inFlux->R2plus());
-    M22_RSR(0, 1) = VecMatVecProduct(m_outFlux->R2plus(), ff_BA, m_inFlux->R2plus());
-    M22_RSR(1, 0) = -VecMatVecProduct(m_outFlux->R2min(), ff_BA, m_inFlux->R2min());
-    M22_RSR(1, 1) = VecMatVecProduct(m_outFlux->R2plus(), ff_BA, m_inFlux->R2min());
+    M22_RSR(0, 0) = -VecMatVecProduct(outFlux->R2min(), ff_BA, inFlux->R2plus());
+    M22_RSR(0, 1) = VecMatVecProduct(outFlux->R2plus(), ff_BA, inFlux->R2plus());
+    M22_RSR(1, 0) = -VecMatVecProduct(outFlux->R2min(), ff_BA, inFlux->R2min());
+    M22_RSR(1, 1) = VecMatVecProduct(outFlux->R2plus(), ff_BA, inFlux->R2min());
 
     return M11_S + M11_RS + M11_SR + M11_RSR + M12_S + M12_RS + M12_SR + M12_RSR + M21_S + M21_RS
            + M21_SR + M21_RSR + M22_S + M22_RS + M22_SR + M22_RSR;
 }
-
-void ComputeDWBAPol::setFlux(std::unique_ptr<const IFlux> inFlux,
-                             std::unique_ptr<const IFlux> outFlux)
-{
-    m_inFlux = std::move(inFlux);
-    m_outFlux = std::move(outFlux);
-}
diff --git a/Resample/FFCompute/ComputeDWBAPol.h b/Resample/FFCompute/ComputeDWBAPol.h
index d8524d080d66fa9d6261adcbb788ba70d7fbd448..8bcd73ddc9228ee4fac7326933bd843293aa79bf 100644
--- a/Resample/FFCompute/ComputeDWBAPol.h
+++ b/Resample/FFCompute/ComputeDWBAPol.h
@@ -20,7 +20,7 @@
 #ifndef BORNAGAIN_RESAMPLE_FFCOMPUTE_COMPUTEDWBAPOL_H
 #define BORNAGAIN_RESAMPLE_FFCOMPUTE_COMPUTEDWBAPOL_H
 
-#include "Resample/FFCompute/IComputeFF.h"
+#include "Resample/FFCompute/IComputePol.h"
 #include <memory>
 
 class IFlux;
@@ -29,27 +29,19 @@ class IFlux;
 
 //! @ingroup formfactors_internal
 
-class ComputeDWBAPol : public IComputeFF {
+class ComputeDWBAPol : public IComputePol {
 public:
-    ComputeDWBAPol(const IFormFactor& ff, size_t iLayer);
+    ComputeDWBAPol(const IFormFactor& ff, size_t i_layer);
     ~ComputeDWBAPol() override;
 
     ComputeDWBAPol* clone() const override;
 
-    //! Throws not-implemented exception.
-    complex_t computeFF(const WavevectorInfo& wavevectors) const override;
-
     //! Returns the coherent sum of the four DWBA terms for polarized scattering.
-    Eigen::Matrix2cd computePolFF(const WavevectorInfo& wavevectors) const override;
-
-    void setFlux(std::unique_ptr<const IFlux> inFlux,
-                 std::unique_ptr<const IFlux> outFlux) override;
+    Eigen::Matrix2cd computePolFF(const WavevectorInfo& wavevectors,
+                                const std::unique_ptr<const IFlux>& inFlux,
+                                const std::unique_ptr<const IFlux>& outFlux) const override;
 
     friend class TestPolarizedDWBATerms;
-
-private:
-    std::unique_ptr<const IFlux> m_inFlux;
-    std::unique_ptr<const IFlux> m_outFlux;
 };
 
 #endif // BORNAGAIN_RESAMPLE_FFCOMPUTE_COMPUTEDWBAPOL_H
diff --git a/Resample/FFCompute/IComputeFF.cpp b/Resample/FFCompute/IComputeFF.cpp
index e13c607f6dbda5d5e46b799c7c0bcbf715532425..fd1c63a36b8c8a350a75d841ded89322c70339f8 100644
--- a/Resample/FFCompute/IComputeFF.cpp
+++ b/Resample/FFCompute/IComputeFF.cpp
@@ -13,24 +13,14 @@
 //  ************************************************************************************************
 
 #include "Resample/FFCompute/IComputeFF.h"
-#include "Base/Pixel/SimulationElement.h"
-#include "Base/Utils/Assert.h"
-#include "Resample/Flux/IFlux.h" // required by VS19 compiler
-#include "Resample/Fresnel/IFresnelMap.h"
-#include "Sample/Material/WavevectorInfo.h"
 #include "Sample/Scattering/IFormFactor.h"
-#include <stdexcept>
 
-IComputeFF::IComputeFF(const IFormFactor& ff, size_t iLayer)
-    : m_ff(ff.clone()), m_iLayer(iLayer) {}
-
-IComputeFF::~IComputeFF() = default;
-
-void IComputeFF::setAmbientMaterial(const Material& material)
+IComputeFF::IComputeFF(const IFormFactor& ff, size_t i_layer) : m_ff(ff.clone()), m_i_layer(i_layer)
 {
-    m_ff->setAmbientMaterial(material);
 }
 
+IComputeFF::~IComputeFF() = default;
+
 double IComputeFF::volume() const
 {
     return m_ff->volume();
@@ -50,30 +40,3 @@ double IComputeFF::topZ(const IRotation& rotation) const
 {
     return m_ff->topZ(rotation);
 }
-
-Eigen::Matrix2cd IComputeFF::computePolFF(const WavevectorInfo&) const
-{
-    throw std::runtime_error("Bug: impossible call to FFCompute::evaluatePol");
-}
-
-void IComputeFF::setFlux(std::unique_ptr<const IFlux>, std::unique_ptr<const IFlux>) {}
-
-complex_t IComputeFF::coherentFF(const SimulationElement& ele)
-{
-    const IFresnelMap* const fresnel_map = ele.fresnelMap();
-    ASSERT(fresnel_map);
-    auto inFlux = fresnel_map->getInFlux(ele.getKi(), m_iLayer);
-    auto outFlux = fresnel_map->getOutFlux(ele.getMeanKf(), m_iLayer);
-    setFlux(std::move(inFlux), std::move(outFlux));
-    return computeFF({ele});
-}
-
-Eigen::Matrix2cd IComputeFF::coherentPolFF(const SimulationElement& ele)
-{
-    const IFresnelMap* const fresnel_map = ele.fresnelMap();
-    ASSERT(fresnel_map);
-    auto inFlux = fresnel_map->getInFlux(ele.getKi(), m_iLayer);
-    auto outFlux = fresnel_map->getOutFlux(ele.getMeanKf(), m_iLayer);
-    setFlux(std::move(inFlux), std::move(outFlux));
-    return computePolFF({ele});
-}
diff --git a/Resample/FFCompute/IComputeFF.h b/Resample/FFCompute/IComputeFF.h
index 23347ae46514b7b20da1cfe3113e3bb63858cc8b..29a3b63e652c283f5d77c518063e900c041265a2 100644
--- a/Resample/FFCompute/IComputeFF.h
+++ b/Resample/FFCompute/IComputeFF.h
@@ -24,12 +24,10 @@
 #include <Eigen/Core>
 #include <memory>
 
-class IFormFactor;
 class IFlux;
+class IFormFactor;
 class IRotation;
 class Material;
-class SimulationElement;
-class WavevectorInfo;
 
 //! Abstract base class for form factor evaluations.
 //!
@@ -42,28 +40,18 @@ public:
     virtual ~IComputeFF();
     virtual IComputeFF* clone() const = 0;
 
-    virtual void setAmbientMaterial(const Material& material);
-
-    size_t iLayer() const { return m_iLayer; }
+    size_t iLayer() const { return m_i_layer; }
 
     virtual double volume() const;
     virtual double radialExtension() const;
     virtual double bottomZ(const IRotation& rotation) const;
     virtual double topZ(const IRotation& rotation) const;
-    virtual complex_t computeFF(const WavevectorInfo& wavevectors) const = 0;
-    //! Returns scattering amplitude for matrix interactions
-    virtual Eigen::Matrix2cd computePolFF(const WavevectorInfo& wavevectors) const;
-    //! Sets reflection/transmission info
-    virtual void setFlux(std::unique_ptr<const IFlux>, std::unique_ptr<const IFlux>);
-
-    complex_t coherentFF(const SimulationElement& sim_element); // TODO make const
-    Eigen::Matrix2cd coherentPolFF(const SimulationElement& sim_element);
 
 protected:
-    IComputeFF(const IFormFactor& ff, size_t iLayer);
+    IComputeFF(const IFormFactor& ff, size_t i_layer);
 
-    std::unique_ptr<IFormFactor> m_ff;
-    const size_t m_iLayer;
+    std::unique_ptr<const IFormFactor> m_ff;
+    const size_t m_i_layer;
 };
 
 #endif // BORNAGAIN_RESAMPLE_FFCOMPUTE_ICOMPUTEFF_H
diff --git a/Resample/FFCompute/IComputePol.cpp b/Resample/FFCompute/IComputePol.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f9e3729b6e359b20c88208dbd2e3437107e39a47
--- /dev/null
+++ b/Resample/FFCompute/IComputePol.cpp
@@ -0,0 +1,30 @@
+//  ************************************************************************************************
+//
+//  BornAgain: simulate and fit reflection and scattering
+//
+//! @file      Resample/FFCompute/IComputePol.cpp
+//! @brief     Implements class ComputeBA.
+//!
+//! @homepage  http://www.bornagainproject.org
+//! @license   GNU General Public License v3 or higher (see COPYING)
+//! @copyright Forschungszentrum Jülich GmbH 2018
+//! @authors   Scientific Computing Group at MLZ (see CITATION, AUTHORS)
+//
+//  ************************************************************************************************
+
+#include "Resample/FFCompute/IComputePol.h"
+#include "Base/Utils/Assert.h"
+#include "Base/Pixel/SimulationElement.h"
+#include "Resample/Flux/IFlux.h" // required by VS19 compiler
+#include "Resample/Fresnel/IFresnelMap.h"
+#include "Sample/Material/WavevectorInfo.h"
+
+
+Eigen::Matrix2cd IComputePol::coherentPolFF(const SimulationElement& ele) const
+{
+    const IFresnelMap* const fresnel_map = ele.fresnelMap();
+    ASSERT(fresnel_map);
+    auto inFlux = fresnel_map->getInFlux(ele.getKi(), m_i_layer);
+    auto outFlux = fresnel_map->getOutFlux(ele.getMeanKf(), m_i_layer);
+    return computePolFF({ele}, inFlux, outFlux);
+}
diff --git a/Resample/FFCompute/IComputePol.h b/Resample/FFCompute/IComputePol.h
new file mode 100644
index 0000000000000000000000000000000000000000..15f9ed5a3d5d0e703a7ac0a46c00e83caa4e2a7d
--- /dev/null
+++ b/Resample/FFCompute/IComputePol.h
@@ -0,0 +1,46 @@
+//  ************************************************************************************************
+//
+//  BornAgain: simulate and fit reflection and scattering
+//
+//! @file      Resample/FFCompute/IComputePol.h
+//! @brief     Defines class ComputeBA.
+//!
+//! @homepage  http://www.bornagainproject.org
+//! @license   GNU General Public License v3 or higher (see COPYING)
+//! @copyright Forschungszentrum Jülich GmbH 2018
+//! @authors   Scientific Computing Group at MLZ (see CITATION, AUTHORS)
+//
+//  ************************************************************************************************
+
+#ifdef SWIG
+#error no need to expose this header to Swig
+#endif
+
+#ifndef USER_API
+#ifndef BORNAGAIN_RESAMPLE_FFCOMPUTE_ICOMPUTEPOL_H
+#define BORNAGAIN_RESAMPLE_FFCOMPUTE_ICOMPUTEPOL_H
+
+#include "Resample/FFCompute/IComputeFF.h"
+
+class SimulationElement;
+class WavevectorInfo;
+
+//! Provides polarized form factor evaluation for given IFormFactor.
+//! Inherited by ComputeBAPol and ComputeDWBAPol.
+
+//! @ingroup formfactors_internal
+
+class IComputePol : public IComputeFF {
+public:
+    Eigen::Matrix2cd coherentPolFF(const SimulationElement& sim_element) const;
+    //! Returns scattering amplitude for matrix interactions
+    virtual Eigen::Matrix2cd computePolFF(const WavevectorInfo& wavevectors,
+                                const std::unique_ptr<const IFlux>& inFlux,
+                                const std::unique_ptr<const IFlux>& outFlux) const = 0;
+
+protected:
+    IComputePol(const IFormFactor& ff, size_t i_layer) : IComputeFF(ff, i_layer) {}
+};
+
+#endif // BORNAGAIN_RESAMPLE_FFCOMPUTE_ICOMPUTEPOL_H
+#endif // USER_API
diff --git a/Resample/FFCompute/IComputeScalar.cpp b/Resample/FFCompute/IComputeScalar.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..78faf173cc993b4c661f64f9232494b8037e821c
--- /dev/null
+++ b/Resample/FFCompute/IComputeScalar.cpp
@@ -0,0 +1,30 @@
+//  ************************************************************************************************
+//
+//  BornAgain: simulate and fit reflection and scattering
+//
+//! @file      Resample/FFCompute/IComputeScalar.cpp
+//! @brief     Implements class ComputeBA.
+//!
+//! @homepage  http://www.bornagainproject.org
+//! @license   GNU General Public License v3 or higher (see COPYING)
+//! @copyright Forschungszentrum Jülich GmbH 2018
+//! @authors   Scientific Computing Group at MLZ (see CITATION, AUTHORS)
+//
+//  ************************************************************************************************
+
+#include "Resample/FFCompute/IComputeScalar.h"
+#include "Base/Utils/Assert.h"
+#include "Base/Pixel/SimulationElement.h"
+#include "Resample/Flux/IFlux.h" // required by VS19 compiler
+#include "Resample/Fresnel/IFresnelMap.h"
+#include "Sample/Material/WavevectorInfo.h"
+
+
+complex_t IComputeScalar::coherentFF(const SimulationElement& ele) const
+{
+    const IFresnelMap* const fresnel_map = ele.fresnelMap();
+    ASSERT(fresnel_map);
+    auto inFlux = fresnel_map->getInFlux(ele.getKi(), m_i_layer);
+    auto outFlux = fresnel_map->getOutFlux(ele.getMeanKf(), m_i_layer);
+    return computeFF({ele}, inFlux, outFlux);
+}
diff --git a/Resample/FFCompute/IComputeScalar.h b/Resample/FFCompute/IComputeScalar.h
new file mode 100644
index 0000000000000000000000000000000000000000..14fd60c75b6be83e2840964f5b5c2e6354dcd3ee
--- /dev/null
+++ b/Resample/FFCompute/IComputeScalar.h
@@ -0,0 +1,44 @@
+//  ************************************************************************************************
+//
+//  BornAgain: simulate and fit reflection and scattering
+//
+//! @file      Resample/FFCompute/IComputeScalar.h
+//! @brief     Defines class ComputeBA.
+//!
+//! @homepage  http://www.bornagainproject.org
+//! @license   GNU General Public License v3 or higher (see COPYING)
+//! @copyright Forschungszentrum Jülich GmbH 2018
+//! @authors   Scientific Computing Group at MLZ (see CITATION, AUTHORS)
+//
+//  ************************************************************************************************
+
+#ifdef SWIG
+#error no need to expose this header to Swig
+#endif
+
+#ifndef USER_API
+#ifndef BORNAGAIN_RESAMPLE_FFCOMPUTE_ICOMPUTESCALAR_H
+#define BORNAGAIN_RESAMPLE_FFCOMPUTE_ICOMPUTESCALAR_H
+
+#include "Resample/FFCompute/IComputeFF.h"
+
+class SimulationElement;
+class WavevectorInfo;
+
+//! Provides scalar form factor evaluation for given IFormFactor.
+//! Inherited by ComputeBA and ComputeDWBA.
+
+//! @ingroup formfactors_internal
+
+class IComputeScalar : public IComputeFF {
+public:
+    virtual complex_t computeFF(const WavevectorInfo& wavevectors,
+                                const std::unique_ptr<const IFlux>& inFlux,
+                                const std::unique_ptr<const IFlux>& outFlux) const = 0;
+    complex_t coherentFF(const SimulationElement& sim_element) const;
+protected:
+    IComputeScalar(const IFormFactor& ff, size_t i_layer) : IComputeFF(ff, i_layer) {}
+};
+
+#endif // BORNAGAIN_RESAMPLE_FFCOMPUTE_ICOMPUTESCALAR_H
+#endif // USER_API
diff --git a/Resample/Fresnel/IFresnelMap.h b/Resample/Fresnel/IFresnelMap.h
index 4e6e26355f7e91b27b7a2a6fc9d458caab7c2ef3..bb6cffe09df5d471f432023107c56a7ad1e583f7 100644
--- a/Resample/Fresnel/IFresnelMap.h
+++ b/Resample/Fresnel/IFresnelMap.h
@@ -40,12 +40,11 @@ public:
     virtual ~IFresnelMap();
 
     //! Retrieves the amplitude coefficients for given wavevector and layer.
-    virtual std::unique_ptr<const IFlux> getInFlux(const kvector_t& kvec,
-                                                   size_t iLayer) const = 0;
+    virtual std::unique_ptr<const IFlux> getInFlux(const kvector_t& kvec, size_t i_layer) const = 0;
 
     //! Retrieves the amplitude coefficients for a (time-reversed) outgoing wavevector.
     virtual std::unique_ptr<const IFlux> getOutFlux(const kvector_t& kvec,
-                                                    size_t iLayer) const = 0;
+                                                    size_t i_layer) const = 0;
 
     const SliceStack& slices() const;
 
diff --git a/Resample/Fresnel/MatrixFresnelMap.cpp b/Resample/Fresnel/MatrixFresnelMap.cpp
index e6c8a5524a8ff334195a99304d223f265abd0e6a..603fc55e33995a7f521303958e537fe380a9378b 100644
--- a/Resample/Fresnel/MatrixFresnelMap.cpp
+++ b/Resample/Fresnel/MatrixFresnelMap.cpp
@@ -39,30 +39,29 @@ size_t MatrixFresnelMap::HashKVector::operator()(const kvector_t& kvec) const no
 }
 
 std::unique_ptr<const IFlux> MatrixFresnelMap::getOutFlux(const kvector_t& kvec,
-                                                          size_t iLayer) const
+                                                          size_t i_layer) const
 {
-    return computeFlux(-kvec, iLayer, m_inverted_slices, m_hash_table_out);
+    return computeFlux(-kvec, i_layer, m_inverted_slices, m_hash_table_out);
 }
 
 std::unique_ptr<const IFlux> MatrixFresnelMap::getInFlux(const kvector_t& kvec,
-                                                         size_t iLayer) const
+                                                         size_t i_layer) const
 {
-    return computeFlux(kvec, iLayer, m_slices, m_hash_table_in);
+    return computeFlux(kvec, i_layer, m_slices, m_hash_table_in);
 }
 
-std::unique_ptr<const IFlux> MatrixFresnelMap::computeFlux(const kvector_t& kvec,
-                                                           size_t iLayer,
+std::unique_ptr<const IFlux> MatrixFresnelMap::computeFlux(const kvector_t& kvec, size_t i_layer,
                                                            const SliceStack& slices,
                                                            CoefficientHash& hash_table) const
 {
     if (!m_use_cache) {
         const auto coeffs = m_strategy->Execute(slices, kvec);
-        return std::unique_ptr<const IFlux>(coeffs[iLayer]->clone());
+        return std::unique_ptr<const IFlux>(coeffs[i_layer]->clone());
     }
     // from cache
     auto it = hash_table.find(kvec);
     if (it == hash_table.end())
         it = hash_table.emplace(kvec, m_strategy->Execute(slices, kvec)).first;
     const auto& coef_vector = it->second;
-    return std::unique_ptr<const IFlux>(coef_vector[iLayer]->clone());
+    return std::unique_ptr<const IFlux>(coef_vector[i_layer]->clone());
 }
diff --git a/Resample/Fresnel/MatrixFresnelMap.h b/Resample/Fresnel/MatrixFresnelMap.h
index 5985232c847b4431b65ba661762518934f36ca77..958a4849b21d431f2678dffa7e645d5805e43b2a 100644
--- a/Resample/Fresnel/MatrixFresnelMap.h
+++ b/Resample/Fresnel/MatrixFresnelMap.h
@@ -48,14 +48,12 @@ private:
         size_t operator()(const kvector_t& kvec) const noexcept;
     };
 
-    std::unique_ptr<const IFlux> getInFlux(const kvector_t& kvec,
-                                           size_t iLayer) const override;
-    std::unique_ptr<const IFlux> getOutFlux(const kvector_t& kvec,
-                                            size_t iLayer) const override;
+    std::unique_ptr<const IFlux> getInFlux(const kvector_t& kvec, size_t i_layer) const override;
+    std::unique_ptr<const IFlux> getOutFlux(const kvector_t& kvec, size_t i_layer) const override;
 
     using CoefficientHash = std::unordered_map<kvector_t, ISpecularStrategy::coeffs_t, HashKVector>;
 
-    std::unique_ptr<const IFlux> computeFlux(const kvector_t& kvec, size_t iLayer,
+    std::unique_ptr<const IFlux> computeFlux(const kvector_t& kvec, size_t i_layer,
                                              const SliceStack& slices,
                                              CoefficientHash& hash_table) const;
     SliceStack m_inverted_slices;
diff --git a/Resample/Fresnel/ScalarFresnelMap.cpp b/Resample/Fresnel/ScalarFresnelMap.cpp
index c3cc474e02cdeff0e6fb10e7e4720af564022c57..948a1e9a9672db01e050eb2335393266a3e9d1e0 100644
--- a/Resample/Fresnel/ScalarFresnelMap.cpp
+++ b/Resample/Fresnel/ScalarFresnelMap.cpp
@@ -32,17 +32,17 @@ ScalarFresnelMap::Hash2Doubles::operator()(const std::pair<double, double>& doub
 }
 
 std::unique_ptr<const IFlux> ScalarFresnelMap::getOutFlux(const kvector_t& kvec,
-                                                          size_t iLayer) const
+                                                          size_t i_layer) const
 {
-    return getInFlux(-kvec, iLayer);
+    return getInFlux(-kvec, i_layer);
 }
 
 std::unique_ptr<const IFlux> ScalarFresnelMap::getInFlux(const kvector_t& kvec,
-                                                         size_t iLayer) const
+                                                         size_t i_layer) const
 {
     if (!m_use_cache) {
         auto coeffs = m_strategy->Execute(m_slices, kvec);
-        return std::unique_ptr<const IFlux>(coeffs[iLayer]->clone());
+        return std::unique_ptr<const IFlux>(coeffs[i_layer]->clone());
     }
     // from cache
     std::pair<double, double> k2_theta(kvec.mag2(), kvec.theta());
@@ -50,5 +50,5 @@ std::unique_ptr<const IFlux> ScalarFresnelMap::getInFlux(const kvector_t& kvec,
     if (it == m_cache.end())
         it = m_cache.emplace(k2_theta, m_strategy->Execute(m_slices, kvec)).first;
     const auto& coef_vector = it->second;
-    return std::unique_ptr<const IFlux>(coef_vector[iLayer]->clone());
+    return std::unique_ptr<const IFlux>(coef_vector[i_layer]->clone());
 }
diff --git a/Resample/Fresnel/ScalarFresnelMap.h b/Resample/Fresnel/ScalarFresnelMap.h
index fb42c647f2b4f945cf8ab6e7d854c055e32efd60..e1a3f8febf5fb4071c65b032291f0c998933c382 100644
--- a/Resample/Fresnel/ScalarFresnelMap.h
+++ b/Resample/Fresnel/ScalarFresnelMap.h
@@ -48,10 +48,8 @@ private:
         size_t operator()(const std::pair<double, double>& doubles) const noexcept;
     };
 
-    std::unique_ptr<const IFlux> getInFlux(const kvector_t& kvec,
-                                           size_t iLayer) const override;
-    std::unique_ptr<const IFlux> getOutFlux(const kvector_t& kvec,
-                                            size_t iLayer) const override;
+    std::unique_ptr<const IFlux> getInFlux(const kvector_t& kvec, size_t i_layer) const override;
+    std::unique_ptr<const IFlux> getOutFlux(const kvector_t& kvec, size_t i_layer) const override;
 
     mutable std::unordered_map<std::pair<double, double>, ISpecularStrategy::coeffs_t, Hash2Doubles>
         m_cache;
diff --git a/Resample/Processed/ParticleRegions.cpp b/Resample/Processed/ParticleRegions.cpp
index 1c91f80a2331b3d18267f0b3bf45601d629148dc..68fd14c37c572507ab898e079f773c2f90667856 100644
--- a/Resample/Processed/ParticleRegions.cpp
+++ b/Resample/Processed/ParticleRegions.cpp
@@ -45,8 +45,8 @@ public:
     std::vector<ZLimits> layerZLimits() const;
 
 private:
-    size_t iLayerTop(double top_z) const;
-    size_t iLayerBottom(double bottom_z) const;
+    size_t i_layerTop(double top_z) const;
+    size_t i_layerBottom(double bottom_z) const;
     void updateLayerLimits(size_t i_layer, ZLimits limits);
     std::vector<double> m_layers_bottomz;
     std::vector<ZLimits> m_layer_fill_limits;
@@ -99,8 +99,8 @@ void LayerFillLimits::update(ZLimits particle_limits, double offset)
         throw std::runtime_error("LayerFillLimits::update: lower_limit > upper_limit.");
     if (bottom == top) // zero-size particle
         return;
-    size_t top_index = iLayerTop(top);
-    size_t bottom_index = iLayerBottom(bottom);
+    size_t top_index = i_layerTop(top);
+    size_t bottom_index = i_layerBottom(bottom);
     for (size_t i_layer = top_index; i_layer < bottom_index + 1; ++i_layer) {
         ZLimits limits(bottom, top);
         updateLayerLimits(i_layer, limits);
@@ -112,7 +112,7 @@ std::vector<ZLimits> LayerFillLimits::layerZLimits() const
     return m_layer_fill_limits;
 }
 
-size_t LayerFillLimits::iLayerTop(double top_z) const
+size_t LayerFillLimits::i_layerTop(double top_z) const
 {
     if (m_layers_bottomz.empty())
         return 0;
@@ -122,7 +122,7 @@ size_t LayerFillLimits::iLayerTop(double top_z) const
     return static_cast<size_t>(m_layers_bottomz.rend() - index_above);
 }
 
-size_t LayerFillLimits::iLayerBottom(double bottom_z) const
+size_t LayerFillLimits::i_layerBottom(double bottom_z) const
 {
     if (m_layers_bottomz.empty())
         return 0;
diff --git a/Resample/Processed/ProcessedLayout.cpp b/Resample/Processed/ProcessedLayout.cpp
index 90b2346d3aa95a49c985a56e41f93ee633d01b06..d648d83ac955366f54e304440f994fd3c40cc2ac 100644
--- a/Resample/Processed/ProcessedLayout.cpp
+++ b/Resample/Processed/ProcessedLayout.cpp
@@ -26,6 +26,7 @@
 #include "Sample/Aggregate/ParticleLayout.h"
 #include "Sample/Material/HomogeneousRegion.h"
 #include "Sample/Particle/IParticle.h"
+#include "Sample/Scattering/IFormFactor.h"
 
 
 ProcessedLayout::ProcessedLayout(const ParticleLayout& layout, const SliceStack& slices,
@@ -37,9 +38,8 @@ ProcessedLayout::ProcessedLayout(const ParticleLayout& layout, const SliceStack&
     const double layout_abundance = layout.getTotalAbundance();
     ASSERT(layout_abundance > 0);
     for (const auto* particle : layout.particles())
-        m_formfactors.emplace_back(
-            processParticle(*particle, slices, z_ref, m_surface_density / layout_abundance,
-                            layout_abundance));
+        m_formfactors.emplace_back(processParticle(
+            *particle, slices, z_ref, m_surface_density / layout_abundance, layout_abundance));
 
     if (const auto* iff = layout.interferenceFunction())
         m_iff.reset(iff->clone());
@@ -97,34 +97,32 @@ CoherentFFSum ProcessedLayout::processParticle(const IParticle& particle, const
         for (auto& region : entry.second)
             region.m_volume *= abundance * intensity_factor;
     for (const auto& entry : region_map) {
-        const size_t iLayer = entry.first;
+        const size_t i_layer = entry.first;
         const auto& regions = entry.second;
-        m_region_map[iLayer].insert(m_region_map[iLayer].begin(), regions.begin(),
-                                         regions.end());
+        m_region_map[i_layer].insert(m_region_map[i_layer].begin(), regions.begin(), regions.end());
     }
 
     std::vector<CoherentFFTerm> terms;
     for (size_t i = 0; i < sliced_ffs.size(); ++i) { // TODO provide slices_ffs.cbegin() etc
         const auto pair = sliced_ffs.getPair(i);
-        const IFormFactor& ff = *pair.first;
-        const size_t iLayer = pair.second;
+        const size_t i_layer = pair.second;
+        const Material& material = slices[i_layer].material();
+        const std::unique_ptr<IFormFactor> ff(pair.first->clone());
+        ff->setAmbientMaterial(material);
+
         std::unique_ptr<IComputeFF> computer;
         if (slices.size() > 1) {
             if (m_polarized)
-                computer = std::make_unique<ComputeDWBAPol>(ff, iLayer);
+                computer = std::make_unique<ComputeDWBAPol>(*ff, i_layer);
             else
-                computer = std::make_unique<ComputeDWBA>(ff, iLayer);
+                computer = std::make_unique<ComputeDWBA>(*ff, i_layer);
         } else {
             // no need for DWBA, use BA
             if (m_polarized)
-                computer = std::make_unique<ComputeBAPol>(ff, iLayer);
+                computer = std::make_unique<ComputeBAPol>(*ff, i_layer);
             else
-                computer = std::make_unique<ComputeBA>(ff, iLayer);
+                computer = std::make_unique<ComputeBA>(*ff, i_layer);
         }
-
-        const Material& slice_material = slices[iLayer].material();
-        computer->setAmbientMaterial(slice_material);
-
         terms.emplace_back(CoherentFFTerm(computer.release()));
     }
     return CoherentFFSum(abundance / total_abundance, terms);
diff --git a/mvvm/tests/testviewmodel/TestToyLayerItem.cpp b/mvvm/tests/testviewmodel/TestToyLayerItem.cpp
index 59a8a758231a063c6265ee4a0941a5cd789c2d93..b10c233b2d5f42b763ef95f0f81321b036942d35 100644
--- a/mvvm/tests/testviewmodel/TestToyLayerItem.cpp
+++ b/mvvm/tests/testviewmodel/TestToyLayerItem.cpp
@@ -67,22 +67,22 @@ TEST_F(ToyLayerItemTest, inViewModel)
     EXPECT_EQ(viewModel.columnCount(), 2);
 
     // accessing to viewItem representing layerItem
-    QModelIndex iLayer = viewModel.index(0, 0);
-    auto viewItem = dynamic_cast<ViewLabelItem*>(viewModel.itemFromIndex(iLayer));
+    QModelIndex i_layer = viewModel.index(0, 0);
+    auto viewItem = dynamic_cast<ViewLabelItem*>(viewModel.itemFromIndex(i_layer));
     EXPECT_TRUE(viewItem != nullptr);
     EXPECT_EQ(viewItem->item(), layerItem);
 
     // it has two rows and two columns, corresponding to our "thickness" and "color" properties
-    EXPECT_EQ(viewModel.rowCount(iLayer), 2);
-    EXPECT_EQ(viewModel.columnCount(iLayer), 2);
+    EXPECT_EQ(viewModel.rowCount(i_layer), 2);
+    EXPECT_EQ(viewModel.columnCount(i_layer), 2);
 
     // accessing to views representing label and value of thickness property
-    QModelIndex thicknessLabelIndex = viewModel.index(0, 0, iLayer);
+    QModelIndex thicknessLabelIndex = viewModel.index(0, 0, i_layer);
     auto thicknessLabelView =
         dynamic_cast<ViewLabelItem*>(viewModel.itemFromIndex(thicknessLabelIndex));
     EXPECT_TRUE(thicknessLabelView != nullptr);
 
-    QModelIndex thicknessValueIndex = viewModel.index(0, 1, iLayer);
+    QModelIndex thicknessValueIndex = viewModel.index(0, 1, i_layer);
     auto thicknessValueView =
         dynamic_cast<ViewDataItem*>(viewModel.itemFromIndex(thicknessValueIndex));
     EXPECT_TRUE(thicknessValueView != nullptr);
@@ -104,8 +104,8 @@ TEST_F(ToyLayerItemTest, layerItemDataChanged)
     // constructing viewModel from sample model
     DefaultViewModel viewModel(&model);
 
-    QModelIndex iLayer = viewModel.index(0, 0);
-    QModelIndex thicknessIndex = viewModel.index(0, 1, iLayer);
+    QModelIndex i_layer = viewModel.index(0, 0);
+    QModelIndex thicknessIndex = viewModel.index(0, 1, i_layer);
 
     QSignalSpy spyDataChanged(&viewModel, &DefaultViewModel::dataChanged);
 
diff --git a/mvvm/tests/testviewmodel/topitemsviewmodel.test.cpp b/mvvm/tests/testviewmodel/topitemsviewmodel.test.cpp
index 03f327569f09208b6e31c2bc1e91b2ac4962b761..fe93c1323404acca0179c8b16b1f09ff972989d5 100644
--- a/mvvm/tests/testviewmodel/topitemsviewmodel.test.cpp
+++ b/mvvm/tests/testviewmodel/topitemsviewmodel.test.cpp
@@ -96,15 +96,15 @@ TEST_F(TopItemsViewModelTest, insertLayerInMultiLayerThenRemove)
 
     // checking their indices
     auto multiiLayer = viewmodel.index(0, 0, QModelIndex());
-    auto iLayer = viewmodel.index(0, 0, multiiLayer);
+    auto i_layer = viewmodel.index(0, 0, multiiLayer);
     EXPECT_EQ(viewmodel.sessionItemFromIndex(multiiLayer), multilayer);
-    EXPECT_EQ(viewmodel.sessionItemFromIndex(iLayer), layer);
+    EXPECT_EQ(viewmodel.sessionItemFromIndex(i_layer), layer);
 
     // checking row and columns
     EXPECT_EQ(viewmodel.rowCount(multiiLayer), 1);
     EXPECT_EQ(viewmodel.columnCount(multiiLayer), 2);
-    EXPECT_EQ(viewmodel.rowCount(iLayer), 0);
-    EXPECT_EQ(viewmodel.columnCount(iLayer), 0);
+    EXPECT_EQ(viewmodel.rowCount(i_layer), 0);
+    EXPECT_EQ(viewmodel.columnCount(i_layer), 0);
 
     // removing layer
     model.removeItem(multilayer, {"", 0});