diff --git a/Sim/Simulation/DepthprobeSimulation.h b/Sim/Simulation/DepthprobeSimulation.h
index 0790eaedca126c402c84da5079965389296324b0..515721e6997a361afa32609c753504c83109ccf2 100644
--- a/Sim/Simulation/DepthprobeSimulation.h
+++ b/Sim/Simulation/DepthprobeSimulation.h
@@ -19,13 +19,8 @@
 #include "Device/Data/Datafield.h"
 #include "Sim/Simulation/ISimulation.h"
 
-#include <memory>
-#include <vector>
-
 class Beam;
 class IAxis;
-class IComputation;
-class ICoordSystem;
 class IFootprintFactor;
 
 //! Simulation of radiation depth profile.
@@ -86,9 +81,6 @@ private:
     void initElementVector() override;
 
     //! Generate a single threaded computation for a given range of simulation elements
-    //! @param re_sample Preprocessed version of our sample
-    //! @param start Index of the first element to include into computation
-    //! @param n_elements Number of elements to process
     std::unique_ptr<IComputation> createComputation(const ReSample& re_sample, size_t start,
                                                     size_t n_elements) override;
 
@@ -99,8 +91,6 @@ private:
     void validateParametrization(const ParameterDistribution& par_distr) const override;
 
     //! Normalize the detector counts to beam intensity, to solid angle, and to exposure angle.
-    //! @param start Index of the first element to operate on
-    //! @param n_elements Number of elements to process
     void normalize(size_t start, size_t n_elements) override;
 
     void addBackgroundIntensity(size_t start, size_t n_elements) override;
@@ -117,7 +107,7 @@ private:
     std::unique_ptr<IAxis> m_z_axis;
     std::vector<std::valarray<double>> m_cache;
     std::vector<DepthprobeElement> m_depth_eles;
-#endif
+#endif // SWIG
 };
 
 #endif // BORNAGAIN_SIM_SIMULATION_DEPTHPROBESIMULATION_H
diff --git a/Sim/Simulation/ISimulation.h b/Sim/Simulation/ISimulation.h
index 85cfdc97d9e3af9aae01ea176a00a532e82823c3..95d94535d2d7adf623c546585af9de9110452836 100644
--- a/Sim/Simulation/ISimulation.h
+++ b/Sim/Simulation/ISimulation.h
@@ -86,26 +86,19 @@ protected:
     DistributionHandler& distributionHandler();
 
 private:
-    virtual void initDistributionHandler() {}
-
-    //! Force polarized computation even in absence of sample magnetization or external fields
-    virtual bool force_polarized() const = 0;
-
+    //... Executor:
     void runSingleSimulation(const ReSample& re_sample, size_t batch_start, size_t batch_size,
                              double weight = 1.0);
 
+    //... Virtual executors:
+    virtual void initDistributionHandler() {}
+
     //! Put into a clean state for running a simulation.
     virtual void prepareSimulation() = 0;
 
-    //! Returns simulation result, based on intensity held in elements vector.
-    virtual SimulationResult packResult() const = 0;
-
     //! Initializes the vector of ISimulation elements.
     virtual void initElementVector() = 0;
 
-    //! Gets the number of elements this simulation needs to calculate.
-    virtual size_t numberOfElements() const = 0;
-
     //! Generate a single threaded computation for a given range of simulation elements.
     //! @param re_sample Preprocessed version of our sample
     //! @param start Index of the first element to include into computation
@@ -113,20 +106,28 @@ private:
     virtual std::unique_ptr<IComputation> createComputation(const ReSample& re_sample, size_t start,
                                                             size_t n_elements) = 0;
 
-    //! Checks the distribution validity for simulation.
-    virtual void validateParametrization(const ParameterDistribution&) const {}
-
-    virtual void addBackgroundIntensity(size_t start_ind, size_t n_elements) = 0;
-
     //! Normalize the detector counts to beam intensity, to solid angle, and to exposure angle.
-    //! @param start_ind Index of the first element to operate on
-    //! @param n_elements Number of elements to process
-    virtual void normalize(size_t start_ind, size_t n_elements) = 0;
+    virtual void normalize(size_t start, size_t n_elements) = 0;
+
+    virtual void addBackgroundIntensity(size_t start, size_t n_elements) = 0;
 
     virtual void addDataToCache(double weight) = 0;
 
     virtual void moveDataFromCache() = 0;
 
+    //... Virtual getters:
+    //! Force polarized computation even in absence of sample magnetization or external fields
+    virtual bool force_polarized() const = 0;
+
+    //! Gets the number of elements this simulation needs to calculate.
+    virtual size_t numberOfElements() const = 0;
+
+    //! Returns simulation result, based on intensity held in elements vector.
+    virtual SimulationResult packResult() const = 0;
+
+    //! Checks the distribution validity for simulation.
+    virtual void validateParametrization(const ParameterDistribution&) const {}
+
     //... Simulation model:
     std::unique_ptr<const MultiLayer> m_sample;
     std::unique_ptr<const IBackground> m_background;
diff --git a/Sim/Simulation/ScatteringSimulation.cpp b/Sim/Simulation/ScatteringSimulation.cpp
index 9011aecda31a9a51054ee0dd2c18d964ac2afe2a..079c70b982219863a8425ea7ee885133fe6173d3 100644
--- a/Sim/Simulation/ScatteringSimulation.cpp
+++ b/Sim/Simulation/ScatteringSimulation.cpp
@@ -39,33 +39,6 @@ ScatteringSimulation::ScatteringSimulation(const Beam& beam, const MultiLayer& s
 
 ScatteringSimulation::~ScatteringSimulation() = default;
 
-SimulationResult ScatteringSimulation::packResult() const
-{
-    auto* const coordsys = simCoordSystem();
-    std::unique_ptr<const Datafield> data(m_detector->createDetectorIntensity(m_eles));
-    return {*data, *coordsys};
-}
-
-const ICoordSystem* ScatteringSimulation::simCoordSystem() const
-{
-    auto* const det2D = dynamic_cast<const IDetector*>(getDetector());
-    ASSERT(det2D);
-    return det2D->scatteringCoords(beam());
-}
-
-void ScatteringSimulation::initElementVector()
-{
-    m_eles = generateElements(beam());
-
-    if (m_cache.empty())
-        m_cache.resize(m_eles.size(), 0.0);
-}
-
-bool ScatteringSimulation::force_polarized() const
-{
-    return m_detector->analyzer().analyzerDirection() != R3{};
-}
-
 void ScatteringSimulation::addMask(const IShape2D& shape, bool mask_value)
 {
     m_detector->addMask(shape, mask_value);
@@ -81,12 +54,11 @@ void ScatteringSimulation::setRegionOfInterest(double xlow, double ylow, double
     m_detector->setRegionOfInterest(xlow, ylow, xup, yup);
 }
 
-void ScatteringSimulation::prepareSimulation()
+const ICoordSystem* ScatteringSimulation::simCoordSystem() const
 {
-    m_active_indices = m_detector->active_indices();
-    m_pixels.reserve(m_active_indices.size());
-    for (auto detector_index : m_active_indices)
-        m_pixels.emplace_back(m_detector->createPixel(detector_index));
+    auto* const det2D = dynamic_cast<const IDetector*>(getDetector());
+    ASSERT(det2D);
+    return det2D->scatteringCoords(beam());
 }
 
 //! init callbacks for setting the parameter values
@@ -114,6 +86,22 @@ void ScatteringSimulation::initDistributionHandler()
     }
 }
 
+void ScatteringSimulation::prepareSimulation()
+{
+    m_active_indices = m_detector->active_indices();
+    m_pixels.reserve(m_active_indices.size());
+    for (auto detector_index : m_active_indices)
+        m_pixels.emplace_back(m_detector->createPixel(detector_index));
+}
+
+void ScatteringSimulation::initElementVector()
+{
+    m_eles = generateElements(beam());
+
+    if (m_cache.empty())
+        m_cache.resize(m_eles.size(), 0.0);
+}
+
 std::unique_ptr<IComputation>
 ScatteringSimulation::createComputation(const ReSample& re_sample, size_t start, size_t n_elements)
 {
@@ -123,25 +111,6 @@ ScatteringSimulation::createComputation(const ReSample& re_sample, size_t start,
                                              begin + static_cast<long>(n_elements));
 }
 
-std::vector<DiffuseElement> ScatteringSimulation::generateElements(const Beam& beam)
-{
-    const double wavelength = beam.wavelength();
-    const double alpha_i = beam.alpha_i();
-    const double phi_i = beam.phi_i();
-    const SpinMatrix beam_polMatrices = beam.polMatrix();
-    const SpinMatrix analyzer_operator = m_detector->analyzer().matrix();
-    const size_t i_specular = m_detector->indexOfSpecular(beam);
-    const size_t N = m_active_indices.size();
-
-    std::vector<DiffuseElement> result;
-    result.reserve(N);
-    for (size_t i = 0; i < N; ++i)
-        result.emplace_back(DiffuseElement(
-            wavelength, alpha_i, phi_i, std::unique_ptr<IPixel>(m_pixels[i]->clone()),
-            beam_polMatrices, analyzer_operator, m_active_indices[i] == i_specular));
-    return result;
-}
-
 void ScatteringSimulation::normalize(size_t start, size_t n_elements)
 {
     const double beam_intensity = m_beam->intensity();
@@ -182,7 +151,38 @@ void ScatteringSimulation::moveDataFromCache()
     m_cache.clear();
 }
 
+bool ScatteringSimulation::force_polarized() const
+{
+    return m_detector->analyzer().analyzerDirection() != R3{};
+}
+
 size_t ScatteringSimulation::numberOfElements() const
 {
     return m_active_indices.size();
 }
+
+SimulationResult ScatteringSimulation::packResult() const
+{
+    auto* const coordsys = simCoordSystem();
+    std::unique_ptr<const Datafield> data(m_detector->createDetectorIntensity(m_eles));
+    return {*data, *coordsys};
+}
+
+std::vector<DiffuseElement> ScatteringSimulation::generateElements(const Beam& beam)
+{
+    const double wavelength = beam.wavelength();
+    const double alpha_i = beam.alpha_i();
+    const double phi_i = beam.phi_i();
+    const SpinMatrix beam_polMatrices = beam.polMatrix();
+    const SpinMatrix analyzer_operator = m_detector->analyzer().matrix();
+    const size_t i_specular = m_detector->indexOfSpecular(beam);
+    const size_t N = m_active_indices.size();
+
+    std::vector<DiffuseElement> result;
+    result.reserve(N);
+    for (size_t i = 0; i < N; ++i)
+        result.emplace_back(DiffuseElement(
+            wavelength, alpha_i, phi_i, std::unique_ptr<IPixel>(m_pixels[i]->clone()),
+            beam_polMatrices, analyzer_operator, m_active_indices[i] == i_specular));
+    return result;
+}
diff --git a/Sim/Simulation/ScatteringSimulation.h b/Sim/Simulation/ScatteringSimulation.h
index 59d8e79abd67c9d7f3a837a9e88b4092c2c29226..800fda0f33688fe43b4af8031184a6322fb1a29c 100644
--- a/Sim/Simulation/ScatteringSimulation.h
+++ b/Sim/Simulation/ScatteringSimulation.h
@@ -70,26 +70,18 @@ public:
     }
 
 private:
-    bool force_polarized() const override;
+    //... Overridden executors:
+    void initDistributionHandler() override;
 
-    //! Put into a clean state for running a simulation
     void prepareSimulation() override;
 
-    void initDistributionHandler() override;
+    void initElementVector() override;
 
     //! Generate a single threaded computation for a given range of simulation elements
-    //! @param re_sample Preprocessed version of our sample
-    //! @param start Index of the first element to include into computation
-    //! @param n_elements Number of elements to process
     std::unique_ptr<IComputation> createComputation(const ReSample& re_sample, size_t start,
                                                     size_t n_elements) override;
 
-    //! Generate simulation elements for given beam
-    std::vector<DiffuseElement> generateElements(const Beam& beam);
-
     //! Normalize the detector counts to beam intensity, to solid angle, and to exposure angle.
-    //! @param start_ind Index of the first element to operate on
-    //! @param n_elements Number of elements to process
     void normalize(size_t start, size_t n_elements) override;
 
     void addBackgroundIntensity(size_t start, size_t n_elements) override;
@@ -98,14 +90,18 @@ private:
 
     void moveDataFromCache() override;
 
-    SimulationResult packResult() const override;
-
-    //! Initializes the vector of ISimulation elements
-    void initElementVector() override;
+    //... Overridden getters:
+    bool force_polarized() const override;
 
     //! Gets the number of elements this simulation needs to calculate
     size_t numberOfElements() const override;
 
+    SimulationResult packResult() const override;
+
+    //... Local function:
+    //! Generate simulation elements for given beam
+    std::vector<DiffuseElement> generateElements(const Beam& beam);
+
     std::shared_ptr<Beam> m_beam;
     std::unique_ptr<IDetector> m_detector;
     std::vector<double> m_cache;