diff --git a/GUI/Model/Sim/InstrumentItems.cpp b/GUI/Model/Sim/InstrumentItems.cpp
index 5f11e824ab637d179202ae6b79f437a5c2b8fff3..086e1871006360c587aef2c80523994d7f4db862 100644
--- a/GUI/Model/Sim/InstrumentItems.cpp
+++ b/GUI/Model/Sim/InstrumentItems.cpp
@@ -72,7 +72,7 @@ const QString ZAxis("ZAxis");
 } // namespace Tag
 
 void setBeamDistribution(ParameterDistribution::WhichParameter which,
-                         const BeamDistributionItem* item, ISimulation* simulation)
+                         const BeamDistributionItem* item, ScatteringSimulation* simulation)
 {
     DistributionItem* di = item->distributionItem();
     ASSERT(di);
diff --git a/GUI/View/Instrument/DepthprobeInstrumentEditor.cpp b/GUI/View/Instrument/DepthprobeInstrumentEditor.cpp
index da38e85fc2f58280e71b96740ce3da35900e38f5..d3c78b3658a0b1845588ad5eedc18120c173928f 100644
--- a/GUI/View/Instrument/DepthprobeInstrumentEditor.cpp
+++ b/GUI/View/Instrument/DepthprobeInstrumentEditor.cpp
@@ -25,7 +25,7 @@ DepthprobeInstrumentEditor::DepthprobeInstrumentEditor(DepthprobeInstrumentItem*
     auto* layout = new QVBoxLayout(this);
     layout->setContentsMargins(0, 0, 0, 0);
 
-    auto* scanEditor = new ScanEditor(this, instrument, instrument->scanItem(), false, false);
+    auto* scanEditor = new ScanEditor(this, instrument, instrument->scanItem(), false, true);
     layout->addWidget(scanEditor);
 
     auto* depthAxisEditor = new AxisForm(this, "Depth axis", &instrument->zAxis(),
diff --git a/GUI/View/Instrument/OffspecInstrumentEditor.cpp b/GUI/View/Instrument/OffspecInstrumentEditor.cpp
index 16908943f0e3bdca0a6cad68e976fae00de201ec..b620ad512240b0151c9ca32fad71e8c70c64e22e 100644
--- a/GUI/View/Instrument/OffspecInstrumentEditor.cpp
+++ b/GUI/View/Instrument/OffspecInstrumentEditor.cpp
@@ -26,7 +26,7 @@ OffspecInstrumentEditor::OffspecInstrumentEditor(OffspecInstrumentItem* instrume
     auto* layout = new QVBoxLayout(this);
     layout->setContentsMargins(0, 0, 0, 0);
 
-    auto* scanEditor = new ScanEditor(this, instrument, instrument->scanItem(), true, false);
+    auto* scanEditor = new ScanEditor(this, instrument, instrument->scanItem(), true, true);
     layout->addWidget(scanEditor);
 
     auto* detectorEditor = new OffspecDetectorEditor(this, instrument);
diff --git a/Resample/Element/SpecularElement.cpp b/Resample/Element/ScanElement.cpp
similarity index 62%
rename from Resample/Element/SpecularElement.cpp
rename to Resample/Element/ScanElement.cpp
index 6fd5f9b1a605c9ab0a19c3c8e745c3ee1bffc516..208ab012efa6402df249ac42e4cae50d5e0f70e3 100644
--- a/Resample/Element/SpecularElement.cpp
+++ b/Resample/Element/ScanElement.cpp
@@ -2,8 +2,8 @@
 //
 //  BornAgain: simulate and fit reflection and scattering
 //
-//! @file      Resample/Element/SpecularElement.cpp
-//! @brief     Implements class the SpecularElement.
+//! @file      Resample/Element/ScanElement.cpp
+//! @brief     Implements class the ScanElement.
 //!
 //! @homepage  http://www.bornagainproject.org
 //! @license   GNU General Public License v3 or higher (see COPYING)
@@ -12,17 +12,19 @@
 //
 //  ************************************************************************************************
 
-#include "Resample/Element/SpecularElement.h"
+#include "Resample/Element/ScanElement.h"
 
-SpecularElement::SpecularElement(size_t i_out, bool computable, double weight, double intensity,
-                                 double footprint, const SpinMatrix& polarizer,
-                                 const SpinMatrix& analyzer, R3 k)
+ScanElement::ScanElement(size_t i_out, bool computable, double weight, double intensity,
+                         double footprint, const SpinMatrix& polarizer, const SpinMatrix& analyzer,
+                         double alpha, double lambda, R3 k)
     : IElement(polarizer, analyzer)
     , m_i_out(i_out)
     , m_computable(computable)
     , m_weight(weight)
     , m_beam_intensity(intensity)
     , m_footprint(footprint)
+    , m_alpha(alpha)
+    , m_lambda(lambda)
     , m_k(k)
 {
 }
diff --git a/Resample/Element/SpecularElement.h b/Resample/Element/ScanElement.h
similarity index 62%
rename from Resample/Element/SpecularElement.h
rename to Resample/Element/ScanElement.h
index 76ca06897f93c4bf2cd469321cff8d51832bde42..39f6e28aacad555387d195c708f6cda4ffb330e4 100644
--- a/Resample/Element/SpecularElement.h
+++ b/Resample/Element/ScanElement.h
@@ -2,8 +2,8 @@
 //
 //  BornAgain: simulate and fit reflection and scattering
 //
-//! @file      Resample/Element/SpecularElement.h
-//! @brief     Declares class the SpecularElement.
+//! @file      Resample/Element/ScanElement.h
+//! @brief     Declares class the ScanElement.
 //!
 //! @homepage  http://www.bornagainproject.org
 //! @license   GNU General Public License v3 or higher (see COPYING)
@@ -15,8 +15,8 @@
 #ifdef SWIG
 #error no need to expose this header to Swig
 #endif // SWIG
-#ifndef BORNAGAIN_RESAMPLE_ELEMENT_SPECULARELEMENT_H
-#define BORNAGAIN_RESAMPLE_ELEMENT_SPECULARELEMENT_H
+#ifndef BORNAGAIN_RESAMPLE_ELEMENT_SCANELEMENT_H
+#define BORNAGAIN_RESAMPLE_ELEMENT_SCANELEMENT_H
 
 #include "Resample/Element/IElement.h"
 #include <functional>
@@ -27,21 +27,23 @@
 
 class SliceStack;
 
-//! Data stucture containing both input and output of a single image pixel for specular simulation.
+//! Data stucture containing input of a single computation point for scan-based simulations.
 
-class SpecularElement : public IElement {
+class ScanElement : public IElement {
 public:
-    SpecularElement(size_t i_out, bool computable, double weight, double intensity,
-                    double footprint, const SpinMatrix& polarizer, const SpinMatrix& analyzer,
-                    R3 k);
+    ScanElement(size_t i_out, bool computable, double weight, double intensity, double footprint,
+                const SpinMatrix& polarizer, const SpinMatrix& analyzer, double alpha,
+                double lambda, R3 k);
 
-    SpecularElement(const SpecularElement& other) = delete;
-    SpecularElement(SpecularElement&& other) noexcept = default;
+    ScanElement(const ScanElement& other) = delete;
+    ScanElement(ScanElement&& other) noexcept = default;
 
     size_t i_out() const { return m_i_out; }
     double weight() const { return m_weight; }
     double beamIntensity() const { return m_beam_intensity; }
     double footprint() const { return m_footprint; }
+    double alpha() const { return m_alpha; }
+    double lambda() const { return m_lambda; }
 
     //! Returns calculation flag (if it's false, zero intensity is assigned to the element)
     bool isCalculated() const { return m_computable; }
@@ -54,7 +56,9 @@ private:
     const double m_weight;
     const double m_beam_intensity;
     const double m_footprint;
+    const double m_alpha;
+    const double m_lambda;
     const R3 m_k;
 };
 
-#endif // BORNAGAIN_RESAMPLE_ELEMENT_SPECULARELEMENT_H
+#endif // BORNAGAIN_RESAMPLE_ELEMENT_SCANELEMENT_H
diff --git a/Sim/Scan/AlphaScan.cpp b/Sim/Scan/AlphaScan.cpp
index 85982d242a24701e6e2396c242353fa9de4a0c2e..cbe73c9f9e50793e7bfbcb15624eac1f4f420e3a 100644
--- a/Sim/Scan/AlphaScan.cpp
+++ b/Sim/Scan/AlphaScan.cpp
@@ -21,26 +21,12 @@
 #include "Device/Beam/Beam.h"
 #include "Device/Beam/IFootprint.h"
 #include "Device/Pol/PolFilter.h"
-#include "Param/Distrib/Distributions.h"
 #include "Param/Distrib/ParameterSample.h"
-#include "Resample/Element/SpecularElement.h"
+#include "Resample/Element/ScanElement.h"
 #include <algorithm> // is_sorted
 
 using PhysConsts::pi;
 
-namespace {
-
-std::vector<ParameterSample> drawDistribution(IDistribution1D* distrib, double point)
-{
-    if (!distrib)
-        return {{point, 1}};
-    distrib->setMean(point);
-    return distrib->distributionSamples();
-}
-
-} // namespace
-
-
 AlphaScan::AlphaScan(const Scale& alpha_axis)
     : PhysicalScan(alpha_axis.clone())
 {
@@ -80,26 +66,26 @@ AlphaScan* AlphaScan::clone() const
     return result;
 }
 
-std::vector<SpecularElement> AlphaScan::generateElements() const
+std::vector<ScanElement> AlphaScan::generateElements() const
 {
-    std::vector<SpecularElement> result;
-    result.reserve(nScan() * nDistributionSamples());
+    std::vector<ScanElement> result;
+    result.reserve(nElements());
 
     for (size_t i = 0; i < m_axis->size(); ++i) {
         const std::vector<ParameterSample> alphaDistrib =
             drawDistribution(m_alpha_distrib.get(), inclinationAt(i) + m_alpha_offset);
         const std::vector<ParameterSample> lambdaDistrib =
             drawDistribution(m_lambda_distrib.get(), wavelengthAt(i));
-        for (auto j : alphaDistrib) {
-            const double alpha = j.value;
-            for (auto k : lambdaDistrib) {
-                const double lambda = k.value;
+        for (auto ad : alphaDistrib) {
+            const double alpha = ad.value;
+            for (auto ld : lambdaDistrib) {
+                const double lambda = ld.value;
                 const bool computable = lambda >= 0 && alpha >= 0 && alpha <= (pi / 2);
-                const double weight = j.weight * k.weight;
+                const double weight = ad.weight * ld.weight;
                 const double footprint = footprintAt(i) ? footprintAt(i)->calculate(alpha) : 1;
                 const R3 veck = vecOfLambdaAlphaPhi(lambda, -alpha);
                 result.emplace_back(i, computable, weight, intensityAt(i), footprint,
-                                    polarizerMatrixAt(i), analyzerMatrix(), veck);
+                                    polarizerMatrixAt(i), analyzerMatrix(), alpha, lambda, veck);
             }
         }
     }
diff --git a/Sim/Scan/AlphaScan.h b/Sim/Scan/AlphaScan.h
index 0404ac39f1a095657462bd46984d2cebe077be65..2b03f00c2440aa069eb09b1adfd67ab3f04356b6 100644
--- a/Sim/Scan/AlphaScan.h
+++ b/Sim/Scan/AlphaScan.h
@@ -17,8 +17,6 @@
 
 #include "Sim/Scan/PhysicalScan.h"
 
-class IDistribution1D;
-
 //! Scan type with inclination angles as coordinate values and a unique wavelength.
 //! Features footprint correction.
 class AlphaScan : public PhysicalScan {
@@ -36,7 +34,7 @@ public:
     AlphaScan* clone() const override;
 
     //! Generates simulation elements for specular simulations
-    std::vector<SpecularElement> generateElements() const override;
+    std::vector<ScanElement> generateElements() const override;
 
     double alphaOffset() const { return m_alpha_offset; }
 
diff --git a/Sim/Scan/BeamScan.cpp b/Sim/Scan/BeamScan.cpp
index 3ac091462495e1b45b388995447482edf0c15321..545de788b27e3d993441c8148f7a7c508e74f40b 100644
--- a/Sim/Scan/BeamScan.cpp
+++ b/Sim/Scan/BeamScan.cpp
@@ -21,6 +21,8 @@
 #include "Device/Beam/Beam.h"
 #include "Device/Beam/IFootprint.h"
 #include "Device/Pol/PolFilter.h"
+#include "Param/Distrib/Distributions.h"
+#include "Param/Distrib/ParameterSample.h"
 #include "Resample/Slice/KzComputation.h"
 
 BeamScan::BeamScan(Scale* axis)
@@ -171,3 +173,12 @@ bool BeamScan::isCommonFootprint() const
     }
     return true;
 }
+
+std::vector<ParameterSample> BeamScan::drawDistribution(IDistribution1D* distrib,
+                                                        double point) const
+{
+    if (!distrib)
+        return {{point, 1}};
+    distrib->setMean(point);
+    return distrib->distributionSamples();
+}
diff --git a/Sim/Scan/BeamScan.h b/Sim/Scan/BeamScan.h
index 5bcd83dc4d36ba10c41a670a2f00f71e1ea59845..b97c46625ae92722e1abe3c2589f6f7d23c22ee2 100644
--- a/Sim/Scan/BeamScan.h
+++ b/Sim/Scan/BeamScan.h
@@ -27,12 +27,14 @@
 
 class Beam;
 class Frame;
+class IDistribution1D;
 class IFootprint;
 class PolFilter;
 class Scale;
+class ScanElement;
 class SliceStack;
-class SpecularElement;
 class SpinMatrix;
+struct ParameterSample;
 
 //! Abstract base class for all types of specular scans.
 
@@ -71,7 +73,7 @@ public:
     SpinMatrix analyzerMatrix() const;
 
     //! Generates simulation elements for specular simulations
-    virtual std::vector<SpecularElement> generateElements() const = 0;
+    virtual std::vector<ScanElement> generateElements() const = 0;
 
     //! Returns coordinate axis assigned to the data holder
     const Scale* coordinateAxis() const { return m_axis.get(); }
@@ -82,6 +84,9 @@ public:
     //! Returns the number of distribution samples
     virtual size_t nDistributionSamples() const = 0;
 
+    //! Returns the number of elements this simulation needs to calculate
+    size_t nElements() const { return nScan() * nDistributionSamples(); }
+
     Frame scanCoordSystem() const;
 
     //! Returns kz values for Abeles computation of reflection/transition coefficients
@@ -93,6 +98,7 @@ public:
 
 protected:
     void copyBeamScan(BeamScan* dest) const; //!< Used by subclass::clone
+    std::vector<ParameterSample> drawDistribution(IDistribution1D* distrib, double point) const;
 
     const std::unique_ptr<Scale> m_axis;
     std::unique_ptr<PolFilter> m_pol_analyzer;
diff --git a/Sim/Scan/LambdaScan.cpp b/Sim/Scan/LambdaScan.cpp
index b14d0f136c447a8c7ffe60eab32a1521961e8594..38e3652a1363ee1c345ed22e36a7de8640155eae 100644
--- a/Sim/Scan/LambdaScan.cpp
+++ b/Sim/Scan/LambdaScan.cpp
@@ -19,7 +19,8 @@
 #include "Base/Vector/GisasDirection.h"
 #include "Device/Beam/Beam.h"
 #include "Device/Beam/IFootprint.h"
-#include "Resample/Element/SpecularElement.h"
+#include "Param/Distrib/ParameterSample.h"
+#include "Resample/Element/ScanElement.h"
 #include <algorithm> // is_sorted
 
 using PhysConsts::pi;
@@ -64,17 +65,29 @@ LambdaScan* LambdaScan::clone() const
     return result;
 }
 
-//! Generates simulation elements for specular simulations
-std::vector<SpecularElement> LambdaScan::generateElements() const
+std::vector<ScanElement> LambdaScan::generateElements() const
 {
-    std::vector<SpecularElement> result;
-    result.reserve(nScan());
+    std::vector<ScanElement> result;
+    result.reserve(nElements());
     for (size_t i = 0; i < m_axis->size(); ++i) {
-        const double kmag = 2 * pi / wavelengthAt(i);
-        const R3 kvec = vecOfKAlphaPhi(kmag, inclinationAt(i), 0);
-        const double footprint = footprintAt(i) ? footprintAt(i)->calculate(inclinationAt(i)) : 1;
-        result.emplace_back(i, kvec.z() <= 0, 1., intensityAt(i), footprint, polarizerMatrixAt(i),
-                            analyzerMatrix(), kvec);
+        const std::vector<ParameterSample> alphaDistrib =
+            drawDistribution(m_alpha_distrib.get(), inclinationAt(i));
+        const std::vector<ParameterSample> lambdaDistrib =
+            drawDistribution(m_lambda_distrib.get(), wavelengthAt(i));
+        for (auto ad : alphaDistrib) {
+            const double alpha = ad.value;
+            for (auto ld : lambdaDistrib) {
+                const double lambda = ld.value;
+                const double kmag = 2 * pi / lambda;
+                const R3 kvec = vecOfKAlphaPhi(kmag, alpha, 0);
+                const bool computable =
+                    lambda >= 0 && alpha >= 0 && alpha <= (pi / 2) && kvec.z() <= 0;
+                const double weight = ad.weight * ld.weight;
+                const double footprint = footprintAt(i) ? footprintAt(i)->calculate(alpha) : 1;
+                result.emplace_back(i, computable, weight, intensityAt(i), footprint,
+                                    polarizerMatrixAt(i), analyzerMatrix(), alpha, lambda, kvec);
+            }
+        }
     }
     return result;
 }
diff --git a/Sim/Scan/LambdaScan.h b/Sim/Scan/LambdaScan.h
index 23ac89d18728c9248c5b50dbbdd99fb0d2dae9c2..f17f25e636297df47ef2f3841141d8e0429568a6 100644
--- a/Sim/Scan/LambdaScan.h
+++ b/Sim/Scan/LambdaScan.h
@@ -36,7 +36,7 @@ public:
     LambdaScan* clone() const override;
 
     //! Generates simulation elements for specular simulations
-    std::vector<SpecularElement> generateElements() const override;
+    std::vector<ScanElement> generateElements() const override;
 
 private:
     LambdaScan(Scale* lambdaScale);
diff --git a/Sim/Scan/QzScan.cpp b/Sim/Scan/QzScan.cpp
index 8c0a17ae74ebd0007674541d065748fff1b60ad4..f4ce4cc16889abaa7c5f6cd105007fbbed24b604 100644
--- a/Sim/Scan/QzScan.cpp
+++ b/Sim/Scan/QzScan.cpp
@@ -21,7 +21,7 @@
 #include "Device/Pol/PolFilter.h"
 #include "Param/Distrib/Distributions.h"
 #include "Param/Distrib/ParameterSample.h"
-#include "Resample/Element/SpecularElement.h"
+#include "Resample/Element/ScanElement.h"
 #include <algorithm> // is_sorted
 
 using PhysConsts::pi;
@@ -88,9 +88,9 @@ std::vector<const INode*> QzScan::nodeChildren() const
 }
 
 //! Generates simulation elements for specular simulations
-std::vector<SpecularElement> QzScan::generateElements() const
+std::vector<ScanElement> QzScan::generateElements() const
 {
-    std::vector<SpecularElement> result;
+    std::vector<ScanElement> result;
     result.reserve(nDistributionSamples());
     for (size_t i = 0; i < m_axis->size(); ++i) {
         const double q0 = m_axis->binCenter(i);
@@ -106,11 +106,11 @@ std::vector<SpecularElement> QzScan::generateElements() const
                 else
                     qz += m_resol_width[0] * j.value;
                 result.emplace_back(i, qz >= 0, j.weight, intensityAt(i), 1., polarizerMatrixAt(i),
-                                    analyzerMatrix(), R3(0, 0, -(qz + m_offset) / 2));
+                                    analyzerMatrix(), 0, 0, R3(0, 0, -(qz + m_offset) / 2));
             }
         } else {
             result.emplace_back(i, q0 >= 0, 1., intensityAt(i), 1., polarizerMatrixAt(i),
-                                analyzerMatrix(), R3(0, 0, -(q0 + m_offset) / 2));
+                                analyzerMatrix(), 0, 0, R3(0, 0, -(q0 + m_offset) / 2));
         }
     }
     return result;
diff --git a/Sim/Scan/QzScan.h b/Sim/Scan/QzScan.h
index fb3684a5051ae4cd1284873f6bfe4b2db5a648b2..4f4bc34edcc47e6be7bbc3c738d91b516eba1f38 100644
--- a/Sim/Scan/QzScan.h
+++ b/Sim/Scan/QzScan.h
@@ -54,7 +54,7 @@ public:
     const IDistribution1D* qzDistribution() const { return m_qz_distrib.get(); }
 
     //! Generates simulation elements for specular simulations
-    std::vector<SpecularElement> generateElements() const override;
+    std::vector<ScanElement> generateElements() const override;
 
     size_t nDistributionSamples() const override;
 
diff --git a/Sim/Simulation/DepthprobeSimulation.cpp b/Sim/Simulation/DepthprobeSimulation.cpp
index bc71d46ce6717c30bc19809bffcebf94ff9c262f..31694f51f251723cf3ef7002b467712789933c79 100644
--- a/Sim/Simulation/DepthprobeSimulation.cpp
+++ b/Sim/Simulation/DepthprobeSimulation.cpp
@@ -24,6 +24,7 @@
 #include "Param/Distrib/DistributionHandler.h"
 #include "Param/Distrib/Distributions.h"
 #include "Resample/Element/IElement.h"
+#include "Resample/Element/ScanElement.h"
 #include "Resample/Flux/ScalarFlux.h"
 #include "Resample/Processed/ReSample.h"
 #include "Sim/Scan/AlphaScan.h"
@@ -60,43 +61,28 @@ std::vector<const INode*> DepthprobeSimulation::nodeChildren() const
 
 //... Overridden executors:
 
-//! init callbacks for setting the parameter values
-void DepthprobeSimulation::initDistributionHandler()
+void DepthprobeSimulation::initScanElementVector()
 {
-    for (const auto& distribution : distributionHandler().paramDistributions()) {
-
-        switch (distribution.whichParameter()) {
-        case ParameterDistribution::BeamInclinationAngle: {
-            distributionHandler().defineCallbackForDistribution(
-                &distribution, [&](double d) { m_scan->setAlphaOffset(d); });
-            break;
-        }
-        case ParameterDistribution::BeamWavelength:
-            distributionHandler().defineCallbackForDistribution(
-                &distribution, [&](double d) { m_scan->setWavelength(d); });
-            break;
-        default:
-            ASSERT_NEVER;
-        }
-    }
+    m_eles = m_scan->generateElements();
 }
 
 void DepthprobeSimulation::runComputation(const ReSample& re_sample, size_t i, double weight)
 {
-    if (m_scan->wavelengthDistribution() || m_scan->alphaDistribution())
-        throw std::runtime_error(
-            "Depthprobe simulation supports neither alpha nor lambda distributions.");
+    const size_t Nz = m_z_axis->size();
+    ASSERT(m_cache.size() / Nz == m_eles.size());
 
-    const size_t n_z = m_z_axis->size();
-    std::valarray<double> intensities; //!< simulated intensity for given z positions
-    intensities.resize(n_z, 0.0);
+    std::valarray<double> z_intensities; //!< simulated intensity for given z positions
+    z_intensities.resize(Nz, 0.0);
 
-    const double result_angle = m_scan->coordinateAxis()->binCenter(i) + m_scan->alphaOffset();
-    if (0 < result_angle && result_angle < (pi / 2)) {
+    ScanElement& scan_ele = *(m_eles.begin() + static_cast<long>(i));
+    const double alpha = scan_ele.alpha();
+    const double lambda = scan_ele.lambda();
+
+    if (scan_ele.isCalculated()) {
         const size_t n_layers = re_sample.numberOfSlices();
-        size_t start_z_ind = n_z;
+        size_t start_z_ind = Nz;
 
-        const R3 ki = vecOfLambdaAlphaPhi(m_scan->wavelengthAt(i), -result_angle);
+        const R3 ki = vecOfLambdaAlphaPhi(lambda, -alpha);
         const Fluxes fluxes = re_sample.fluxesIn(ki);
 
         double z_layer_bottom(0.0);
@@ -130,11 +116,11 @@ void DepthprobeSimulation::runComputation(const ReSample& re_sample, size_t i, d
                 else
                     throw std::runtime_error("Invalid combination of ZDirection flags");
                 if ((m_flags & 12) == WaveProperty_Intensity)
-                    intensities[i_z] = std::norm(psi);
+                    z_intensities[i_z] = std::norm(psi);
                 else if (m_flags & WaveProperty_Modulus)
-                    intensities[i_z] = std::abs(psi);
+                    z_intensities[i_z] = std::abs(psi);
                 else if (m_flags & WaveProperty_Phase)
-                    intensities[i_z] = std::arg(psi);
+                    z_intensities[i_z] = std::arg(psi);
                 else
                     throw std::runtime_error("Invalid combination of WaveProperty flags");
             }
@@ -142,13 +128,11 @@ void DepthprobeSimulation::runComputation(const ReSample& re_sample, size_t i, d
         }
     }
 
-    double intensity_factor = m_scan->intensityAt(i);
-    for (double& v : intensities)
-        v *= intensity_factor;
+    for (double& v : z_intensities)
+        v *= scan_ele.beamIntensity();
 
-    const size_t N1 = m_z_axis->size();
-    for (size_t j = 0; j < N1; ++j)
-        m_cache[j * m_scan->nScan() + i] += intensities[j] * weight;
+    for (size_t j = 0; j < Nz; ++j)
+        m_cache[j * m_scan->nElements() + i] += z_intensities[j] * weight * scan_ele.weight();
 
     progress().incrementDone(1);
 }
@@ -157,7 +141,7 @@ void DepthprobeSimulation::runComputation(const ReSample& re_sample, size_t i, d
 
 size_t DepthprobeSimulation::nElements() const
 {
-    return m_scan->coordinateAxis()->size();
+    return m_scan->nElements();
 }
 
 size_t DepthprobeSimulation::nOutChannels() const
@@ -165,13 +149,23 @@ size_t DepthprobeSimulation::nOutChannels() const
     return nElements() * m_z_axis->size();
 }
 
-Datafield DepthprobeSimulation::packResult()
+Datafield DepthprobeSimulation::packResult() const
 {
-    std::vector<const Scale*> axes{m_scan->coordinateAxis()->clone(), m_z_axis->clone()};
-    auto data = std::make_unique<Datafield>(axes, m_cache);
+    const size_t Nz = m_z_axis->size();
+    const size_t Ns = m_scan->nScan();
+    const size_t Nslong = m_scan->nElements();
+    std::vector<double> out(Ns * Nz, 0.0);
+
+    for (size_t j = 0; j < Nz; ++j)
+        for (size_t i = 0; i < Nslong; ++i) {
+            const ScanElement& ele = m_eles.at(i);
+            const size_t i_real = ele.i_out();
+            out[j * Ns + i_real] += m_cache[j * Nslong + i];
+        }
 
     if (background())
         throw std::runtime_error("nonzero background is not supported by DepthprobeSimulation");
 
-    return {*data};
+    std::vector<const Scale*> axes{m_scan->coordinateAxis()->clone(), m_z_axis->clone()};
+    return {axes, out};
 }
diff --git a/Sim/Simulation/DepthprobeSimulation.h b/Sim/Simulation/DepthprobeSimulation.h
index 26107bd2df425ece66c7e2b0b4c3d37772538ea6..632e1de43ebf0cc446ef54df166f7dd6a0301e61 100644
--- a/Sim/Simulation/DepthprobeSimulation.h
+++ b/Sim/Simulation/DepthprobeSimulation.h
@@ -21,6 +21,7 @@ class AlphaScan;
 class BeamScan;
 class IFootprint;
 class Scale;
+class ScanElement;
 
 extern const int ZDirection_None;
 extern const int ZDirection_Reflected;
@@ -54,9 +55,7 @@ public:
 
 private:
     //... Overridden executors:
-    void initDistributionHandler() override;
-
-    void prepareSimulation() override {}
+    void initScanElementVector() override;
 
     void runComputation(const ReSample& re_sample, size_t iElement, double weight) override;
 
@@ -68,12 +67,15 @@ private:
 
     size_t nOutChannels() const override;
 
-    Datafield packResult() override;
+    Datafield packResult() const override;
 
     //... Model components:
     std::unique_ptr<AlphaScan> m_scan;
     std::unique_ptr<Scale> m_z_axis;
     const int m_flags;
+
+    //... Caches:
+    std::vector<ScanElement> m_eles;
 #endif // SWIG
 };
 
diff --git a/Sim/Simulation/ISimulation.cpp b/Sim/Simulation/ISimulation.cpp
index 189ea708b959e01aeef26aeaf0b49684a5949cf6..31ebc2c5224d2c59ee087f1694d93b285b67cf95 100644
--- a/Sim/Simulation/ISimulation.cpp
+++ b/Sim/Simulation/ISimulation.cpp
@@ -81,13 +81,6 @@ void ISimulation::setBackground(const IBackground& bg)
     m_background.reset(bg.clone());
 }
 
-void ISimulation::addParameterDistribution(ParameterDistribution::WhichParameter whichParameter,
-                                           const IDistribution1D& distribution)
-{
-    ParameterDistribution par_distr(whichParameter, distribution);
-    distributionHandler().addDistribution(par_distr);
-}
-
 void ISimulation::subscribe(const std::function<bool(size_t)>& inform)
 {
     ASSERT(m_progress);
@@ -148,6 +141,7 @@ Datafield ISimulation::simulate()
     if (n_combinations == 1) {
         runSingleSimulation(re_sample, batch_start, batch_size, 1.);
     } else {
+        // only GISAS
         initDistributionHandler();
         for (size_t index = 0; index < n_combinations; ++index) {
             double weight = distributionHandler().setParameterValues(index);
@@ -210,7 +204,7 @@ DistributionHandler& ISimulation::distributionHandler()
 void ISimulation::runSingleSimulation(const ReSample& re_sample, size_t batch_start,
                                       size_t batch_size, double weight)
 {
-    initElementVector();
+    initScanElementVector();
 
     const size_t n_threads = m_options->getNumberOfThreads();
     ASSERT(n_threads > 0);
diff --git a/Sim/Simulation/ISimulation.h b/Sim/Simulation/ISimulation.h
index ff88a4fda61e6b189d7549b2d96961e053e2e7be..71a346a5b467da7346068158f083f1058ecdb323 100644
--- a/Sim/Simulation/ISimulation.h
+++ b/Sim/Simulation/ISimulation.h
@@ -50,9 +50,6 @@ public:
     //... Setters:
     void setBackground(const IBackground& bg);
 
-    void addParameterDistribution(ParameterDistribution::WhichParameter whichParameter,
-                                  const IDistribution1D& distribution);
-
     void setTerminalProgressMonitor();
 
     SimulationOptions& options();
@@ -92,7 +89,7 @@ private:
     virtual void prepareSimulation() {}
 
     //! Initializes the vector of ISimulation elements.
-    virtual void initElementVector() {}
+    virtual void initScanElementVector() {}
 
     //! Runs a single threaded computation for a given range of simulation elements.
     virtual void runComputation(const ReSample& re_sample, size_t iElement, double weight) = 0;
@@ -108,7 +105,7 @@ private:
     virtual size_t nOutChannels() const { return nElements(); }
 
     //! Returns simulation result, based on intensity held in elements vector.
-    virtual Datafield packResult() = 0;
+    virtual Datafield packResult() const = 0;
 
     //... Simulation model:
     std::unique_ptr<const Sample> m_sample;
diff --git a/Sim/Simulation/OffspecSimulation.cpp b/Sim/Simulation/OffspecSimulation.cpp
index db502e2b0a51596e1b5d6cc36c13d1251d25cf0c..a77c320b773bd22cf8353e6b64b0113aac02b32d 100644
--- a/Sim/Simulation/OffspecSimulation.cpp
+++ b/Sim/Simulation/OffspecSimulation.cpp
@@ -24,6 +24,7 @@
 #include "Param/Distrib/DistributionHandler.h"
 #include "Param/Distrib/Distributions.h"
 #include "Resample/Element/DiffuseElement.h"
+#include "Resample/Element/ScanElement.h"
 #include "Sim/Background/IBackground.h"
 #include "Sim/Computation/DWBAComputation.h"
 #include "Sim/Scan/PhysicalScan.h"
@@ -56,59 +57,40 @@ void OffspecSimulation::prepareSimulation()
 
 //... Overridden executors:
 
-//! init callbacks for setting the parameter values
-void OffspecSimulation::initDistributionHandler()
+void OffspecSimulation::initScanElementVector()
 {
-    for (const auto& distribution : distributionHandler().paramDistributions()) {
-
-        switch (distribution.whichParameter()) {
-            /*
-        case ParameterDistribution::BeamWavelength:
-            distributionHandler().defineCallbackForDistribution(
-                &distribution, [&](double d) { m_scan->setWavelength(d); });
-            break;
-            */
-        default:
-            ASSERT_NEVER;
-        }
-    }
+    m_eles = m_scan->generateElements();
 }
 
 void OffspecSimulation::runComputation(const ReSample& re_sample, size_t i, double weight)
 {
-    if (auto* phys_scan = dynamic_cast<PhysicalScan*>(m_scan.get()))
-        if (phys_scan->wavelengthDistribution() || phys_scan->alphaDistribution())
-            throw std::runtime_error(
-                "Offspecular simulation supports neither alpha nor lambda distributions.");
-
-    if (m_cache.empty())
-        m_cache.resize(nElements(), 0.0);
-
     const size_t Na = m_detector->totalSize();
-    size_t j = i / Na; // index in scan
+    size_t j = i / Na; // index in unwrapped scan (including distribution sampling)
     size_t k = i % Na; // index in detector
 
-    const double alpha_i = m_scan->inclinationAt(j);
+    ASSERT(m_cache.size() / Na == m_eles.size());
+
+    ScanElement& scan_ele = *(m_eles.begin() + static_cast<long>(j));
+    const double alpha_i = scan_ele.alpha();
+    const double lambda_i = scan_ele.lambda();
     const double phi_i = 0;
     const bool isSpecular = k == m_detector->indexOfSpecular(alpha_i, phi_i);
 
-    DiffuseElement ele(m_scan->wavelengthAt(j), alpha_i, phi_i, m_pixels[k],
-                       m_scan->polarizerMatrixAt(j), m_detector->analyzer().matrix(), isSpecular);
-
-    double intensity = Compute::scattered_and_reflected(re_sample, options(), ele);
+    DiffuseElement diff_ele(lambda_i, alpha_i, phi_i, m_pixels[k],
+                            m_scan->polarizerMatrixAt(scan_ele.i_out()),
+                            m_detector->analyzer().matrix(), isSpecular);
 
-    if (const auto* footprint = m_scan->footprintAt(j))
-        intensity *= footprint->calculate(alpha_i);
+    double intensity = Compute::scattered_and_reflected(re_sample, options(), diff_ele);
 
     double sin_alpha_i = std::abs(std::sin(alpha_i));
     if (sin_alpha_i == 0.0) {
         intensity = 0;
     } else {
-        const double solid_angle = ele.solidAngle();
-        intensity *= m_scan->intensityAt(j) * solid_angle / sin_alpha_i;
+        const double solid_angle = diff_ele.solidAngle();
+        intensity *= scan_ele.beamIntensity() * scan_ele.footprint() * solid_angle / sin_alpha_i;
     }
 
-    m_cache[i] += intensity * weight;
+    m_cache[i] += intensity * weight * scan_ele.weight();
 
     progress().incrementDone(1);
 }
@@ -122,33 +104,40 @@ bool OffspecSimulation::force_polarized() const
 
 size_t OffspecSimulation::nElements() const
 {
-    return m_detector->totalSize() * m_scan->nScan();
+    return m_detector->totalSize() * m_scan->nElements();
 }
 
-Datafield OffspecSimulation::packResult()
+Datafield OffspecSimulation::packResult() const
 {
     const size_t ns = m_scan->nScan();
+    const size_t nslong = m_scan->nElements();
+
     const size_t nphi = m_detector->axis(0).size();
     const size_t nalp = m_detector->axis(1).size();
     const size_t Ndet = nalp * nphi;
     std::vector<double> out(ns * nalp, 0.);
 
-    // Apply detector resolution and transfer detector image
-    for (size_t j = 0; j < ns; ++j) {
-
-        // TODO restore resolution m_detector->applyDetectorResolution(&detector_image);
-
+    for (size_t j = 0; j < nslong; ++j) {
         for (size_t ia = 0; ia < nalp; ++ia) {
             double val = 0;
+
+            // sum over phi to the distribution point
             for (size_t ip = 0; ip < nphi; ++ip)
                 val += m_cache[j * Ndet + ia * nphi + ip];
-            if (background())
-                val = background()->addBackground(val);
-            out[ia * ns + j] = val;
+
+            // accumulate distribution points to the scan point
+            const ScanElement& ele = m_eles.at(j);
+            const size_t j_real = ele.i_out();
+            out[ia * ns + j_real] += val;
         }
     }
 
-    return {
-        std::vector<const Scale*>{m_scan->coordinateAxis()->clone(), m_detector->axis(1).clone()},
-        out};
+    // TODO restore resolution m_detector->applyDetectorResolution(&detector_image);
+
+    if (background())
+        for (double& o : out)
+            o = background()->addBackground(o);
+
+    std::vector<const Scale*> axes{m_scan->coordinateAxis()->clone(), m_detector->axis(1).clone()};
+    return {axes, out};
 }
diff --git a/Sim/Simulation/OffspecSimulation.h b/Sim/Simulation/OffspecSimulation.h
index 30a8f13d203742adea4fa4dc59b02a4a4709a4bb..d60e770a7641b51c797fd05bd613abf23e5db952 100644
--- a/Sim/Simulation/OffspecSimulation.h
+++ b/Sim/Simulation/OffspecSimulation.h
@@ -22,6 +22,7 @@ class Datafield;
 class OffspecDetector;
 class PhysicalScan;
 class Pixel;
+class ScanElement;
 
 //! Off-specular scattering simulation.
 //!
@@ -44,10 +45,10 @@ public:
 
 private:
     //... Overridden executors:
-    void initDistributionHandler() override;
-
     void prepareSimulation() override;
 
+    void initScanElementVector() override;
+
     void runComputation(const ReSample& re_sample, size_t iElement, double weight) override;
 
     //... Overridden getters:
@@ -55,7 +56,7 @@ private:
 
     size_t nElements() const override;
 
-    Datafield packResult() override;
+    Datafield packResult() const override;
 
     //... Model components:
     std::unique_ptr<PhysicalScan> m_scan;
@@ -63,7 +64,8 @@ private:
 
     //... Caches:
     OwningVector<const Pixel> m_pixels; //!< All unmasked pixels inside ROI.
-#endif                                  // SWIG
+    std::vector<ScanElement> m_eles;
+#endif // SWIG
 };
 
 #endif // BORNAGAIN_SIM_SIMULATION_OFFSPECSIMULATION_H
diff --git a/Sim/Simulation/ScatteringSimulation.cpp b/Sim/Simulation/ScatteringSimulation.cpp
index b4a65572c2213716b270f1298122279cdfe91679..b1aa12b5ce27e0dff529dd7dc7c90c8a7b7b445e 100644
--- a/Sim/Simulation/ScatteringSimulation.cpp
+++ b/Sim/Simulation/ScatteringSimulation.cpp
@@ -38,6 +38,13 @@ ScatteringSimulation::ScatteringSimulation(const Beam& beam, const Sample& sampl
 
 ScatteringSimulation::~ScatteringSimulation() = default;
 
+void ScatteringSimulation::addParameterDistribution(
+    ParameterDistribution::WhichParameter whichParameter, const IDistribution1D& distribution)
+{
+    ParameterDistribution par_distr(whichParameter, distribution);
+    distributionHandler().addDistribution(par_distr);
+}
+
 //... Overridden executors:
 
 //! init callbacks for setting the parameter values
@@ -113,7 +120,7 @@ size_t ScatteringSimulation::nElements() const
     return m_active_indices.size();
 }
 
-Datafield ScatteringSimulation::packResult()
+Datafield ScatteringSimulation::packResult() const
 {
     Datafield result(m_detector->createDetectorMap());
     for (size_t i = 0; i < m_active_indices.size(); ++i)
@@ -121,12 +128,9 @@ Datafield ScatteringSimulation::packResult()
 
     m_detector->applyDetectorResolution(&result);
 
-    if (background()) {
-        for (size_t i : m_active_indices) {
-            double& val = result[i];
-            val = background()->addBackground(val);
-        }
-    }
+    if (background())
+        for (size_t i : m_active_indices)
+            result[i] = background()->addBackground(result[i]);
 
     return result;
 }
diff --git a/Sim/Simulation/ScatteringSimulation.h b/Sim/Simulation/ScatteringSimulation.h
index 2b1be044878b166ad8ce08fc511ad531a1023f4d..bd5b55fcf420621c154467e601ffa69c05259252 100644
--- a/Sim/Simulation/ScatteringSimulation.h
+++ b/Sim/Simulation/ScatteringSimulation.h
@@ -38,6 +38,9 @@ public:
     Beam& beam() { return *m_beam; }
     IDetector& detector() { return *m_detector; }
 
+    void addParameterDistribution(ParameterDistribution::WhichParameter whichParameter,
+                                  const IDistribution1D& distribution);
+
 #ifndef SWIG
     const Beam& beam() const { return *m_beam; }
     const IDetector& detector() const { return *m_detector; }
@@ -57,7 +60,7 @@ private:
     //! Returns the number of elements this simulation needs to calculate
     size_t nElements() const override;
 
-    Datafield packResult() override;
+    Datafield packResult() const override;
 
     //... Model components:
     std::shared_ptr<Beam> m_beam;
diff --git a/Sim/Simulation/SpecularSimulation.cpp b/Sim/Simulation/SpecularSimulation.cpp
index 5d27ddc5d07242937f5214171bf2640870e4413b..5f74b839007ce5ab8e4804955dc617ad4a655f06 100644
--- a/Sim/Simulation/SpecularSimulation.cpp
+++ b/Sim/Simulation/SpecularSimulation.cpp
@@ -22,7 +22,7 @@
 #include "Device/Pol/PolFilter.h"
 #include "Param/Distrib/Distributions.h"
 #include "Param/Distrib/ParameterSample.h"
-#include "Resample/Element/SpecularElement.h"
+#include "Resample/Element/ScanElement.h"
 #include "Resample/Processed/ReSample.h"
 #include "Resample/Specular/ComputeFluxMagnetic.h"
 #include "Resample/Specular/ComputeFluxScalar.h"
@@ -43,14 +43,14 @@ SpecularSimulation::~SpecularSimulation() = default;
 
 //... Overridden executors:
 
-void SpecularSimulation::initElementVector()
+void SpecularSimulation::initScanElementVector()
 {
     m_eles = m_scan->generateElements();
 }
 
 void SpecularSimulation::runComputation(const ReSample& re_sample, size_t iElement, double weight)
 {
-    SpecularElement& ele = *(m_eles.begin() + static_cast<long>(iElement));
+    ScanElement& ele = *(m_eles.begin() + static_cast<long>(iElement));
 
     double refl = 0;
     if (ele.isCalculated()) {
@@ -71,7 +71,7 @@ void SpecularSimulation::runComputation(const ReSample& re_sample, size_t iEleme
         }
     }
 
-    m_cache[iElement] += refl * ele.footprint() * weight;
+    m_cache[iElement] += refl * ele.footprint() * weight * ele.beamIntensity() * ele.weight();
 
     progress().incrementDone(1);
 }
@@ -85,15 +85,15 @@ bool SpecularSimulation::force_polarized() const
 
 size_t SpecularSimulation::nElements() const
 {
-    return m_scan->nScan() * m_scan->nDistributionSamples();
+    return m_scan->nElements();
 }
 
-Datafield SpecularSimulation::packResult()
+Datafield SpecularSimulation::packResult() const
 {
     std::vector<double> vec(m_scan->nScan(), 0.0);
     for (size_t i = 0; i < nElements(); i++) {
-        const SpecularElement& ele = m_eles.at(i);
-        vec.at(ele.i_out()) += ele.beamIntensity() * m_cache.at(i) * ele.weight();
+        const ScanElement& ele = m_eles.at(i);
+        vec.at(ele.i_out()) += m_cache.at(i);
     }
 
     double common_intensity = m_scan->commonIntensity();
@@ -110,6 +110,5 @@ Datafield SpecularSimulation::packResult()
         else
             vec[i] = 0;
 
-    Datafield data(std::vector<const Scale*>{m_scan->coordinateAxis()->clone()}, vec);
-    return {data};
+    return {std::vector<const Scale*>{m_scan->coordinateAxis()->clone()}, vec};
 }
diff --git a/Sim/Simulation/SpecularSimulation.h b/Sim/Simulation/SpecularSimulation.h
index 308b295dace71b68c52494f86ed6aef7b8fe13be..99229d4a1e565f1b70d11e43fb01fe9885300f17 100644
--- a/Sim/Simulation/SpecularSimulation.h
+++ b/Sim/Simulation/SpecularSimulation.h
@@ -18,7 +18,7 @@
 #include "Sim/Simulation/ISimulation.h"
 
 class BeamScan;
-class SpecularElement;
+class ScanElement;
 
 //! Specular reflectometry simulation.
 //!
@@ -38,7 +38,7 @@ public:
 
 private:
     //... Overridden executors:
-    void initElementVector() override;
+    void initScanElementVector() override;
 
     void runComputation(const ReSample& re_sample, size_t iElement, double weight) override;
 
@@ -48,13 +48,13 @@ private:
     //! Returns the number of elements this simulation needs to calculate
     size_t nElements() const override;
 
-    Datafield packResult() override;
+    Datafield packResult() const override;
 
     //... Model components:
     std::unique_ptr<const BeamScan> m_scan;
 
     //... Caches:
-    std::vector<SpecularElement> m_eles;
+    std::vector<ScanElement> m_eles;
 #endif // SWIG
 };
 
diff --git a/auto/Examples/varia/Resonator.py b/auto/Examples/varia/Resonator.py
index 1b3118ad0a794d37ed1f8419996027cdecab8e63..cef29118ef66f1baa71186a2cfe9162fd1d0a6fb 100755
--- a/auto/Examples/varia/Resonator.py
+++ b/auto/Examples/varia/Resonator.py
@@ -72,13 +72,12 @@ def get_simulation(sample):
     scan = ba.AlphaScan(na, ai_min, ai_max)
     scan.setWavelength(wl)
 
+    alpha_distr = ba.DistributionGaussian(0, d_ang, 25, 3.)
+    scan.setAngleDistribution(alpha_distr)
+
     z_axis = ba.EquiDivision("z (nm)", nz, z_min, z_max)
     simulation = ba.DepthprobeSimulation(scan, sample, z_axis)
 
-    alpha_distr = ba.DistributionGaussian(0, d_ang, 25, 3.)
-    simulation.addParameterDistribution(
-        ba.ParameterDistribution.BeamInclinationAngle, alpha_distr)
-
     return simulation
 
 
diff --git a/auto/FigExamples/varia/Resonator.py b/auto/FigExamples/varia/Resonator.py
index c3ee54cad494e57e0cac01f9651b8104ae58494e..00e151eeb8a8b5f7bd00f396fd3a1a0198c9ecf8 100755
--- a/auto/FigExamples/varia/Resonator.py
+++ b/auto/FigExamples/varia/Resonator.py
@@ -72,13 +72,12 @@ def get_simulation(sample):
     scan = ba.AlphaScan(na, ai_min, ai_max)
     scan.setWavelength(wl)
 
+    alpha_distr = ba.DistributionGaussian(0, d_ang, 25, 3.)
+    scan.setAngleDistribution(alpha_distr)
+
     z_axis = ba.EquiDivision("z (nm)", nz, z_min, z_max)
     simulation = ba.DepthprobeSimulation(scan, sample, z_axis)
 
-    alpha_distr = ba.DistributionGaussian(0, d_ang, 25, 3.)
-    simulation.addParameterDistribution(
-        ba.ParameterDistribution.BeamInclinationAngle, alpha_distr)
-
     return simulation
 
 
diff --git a/auto/MiniExamples/varia/Resonator.py b/auto/MiniExamples/varia/Resonator.py
index 62069040fd3177c48bc50275065acd135367628d..148abe2f953e342bf6dd0a9ae60123dbeb554dff 100755
--- a/auto/MiniExamples/varia/Resonator.py
+++ b/auto/MiniExamples/varia/Resonator.py
@@ -72,13 +72,12 @@ def get_simulation(sample):
     scan = ba.AlphaScan(na, ai_min, ai_max)
     scan.setWavelength(wl)
 
+    alpha_distr = ba.DistributionGaussian(0, d_ang, 25, 3.)
+    scan.setAngleDistribution(alpha_distr)
+
     z_axis = ba.EquiDivision("z (nm)", nz, z_min, z_max)
     simulation = ba.DepthprobeSimulation(scan, sample, z_axis)
 
-    alpha_distr = ba.DistributionGaussian(0, d_ang, 25, 3.)
-    simulation.addParameterDistribution(
-        ba.ParameterDistribution.BeamInclinationAngle, alpha_distr)
-
     return simulation
 
 
diff --git a/auto/Wrap/libBornAgainSim.py b/auto/Wrap/libBornAgainSim.py
index efc33a2ee736f5c0af21bbfbae0c3f6182297e90..31d69e0af9921cec02980a3f802c3f96c024629d 100644
--- a/auto/Wrap/libBornAgainSim.py
+++ b/auto/Wrap/libBornAgainSim.py
@@ -2253,10 +2253,6 @@ class ISimulation(libBornAgainParam.INode):
         r"""setBackground(ISimulation self, IBackground bg)"""
         return _libBornAgainSim.ISimulation_setBackground(self, bg)
 
-    def addParameterDistribution(self, whichParameter, distribution):
-        r"""addParameterDistribution(ISimulation self, ParameterDistribution::WhichParameter whichParameter, IDistribution1D const & distribution)"""
-        return _libBornAgainSim.ISimulation_addParameterDistribution(self, whichParameter, distribution)
-
     def setTerminalProgressMonitor(self):
         r"""setTerminalProgressMonitor(ISimulation self)"""
         return _libBornAgainSim.ISimulation_setTerminalProgressMonitor(self)
@@ -2294,6 +2290,10 @@ class ScatteringSimulation(ISimulation):
         r"""detector(ScatteringSimulation self) -> IDetector &"""
         return _libBornAgainSim.ScatteringSimulation_detector(self)
 
+    def addParameterDistribution(self, whichParameter, distribution):
+        r"""addParameterDistribution(ScatteringSimulation self, ParameterDistribution::WhichParameter whichParameter, IDistribution1D const & distribution)"""
+        return _libBornAgainSim.ScatteringSimulation_addParameterDistribution(self, whichParameter, distribution)
+
 # Register ScatteringSimulation in _libBornAgainSim:
 _libBornAgainSim.ScatteringSimulation_swigregister(ScatteringSimulation)
 class DepthprobeSimulation(ISimulation):
diff --git a/auto/Wrap/libBornAgainSim_wrap.cpp b/auto/Wrap/libBornAgainSim_wrap.cpp
index 07dadf4ed0be3c4eb642dadf14e73e420ad03666..6f2cd75c3d02c21eb060eeda91d20228a822b199 100644
--- a/auto/Wrap/libBornAgainSim_wrap.cpp
+++ b/auto/Wrap/libBornAgainSim_wrap.cpp
@@ -30785,57 +30785,6 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_ISimulation_addParameterDistribution(PyObject *self, PyObject *args) {
-  PyObject *resultobj = 0;
-  ISimulation *arg1 = (ISimulation *) 0 ;
-  ParameterDistribution::WhichParameter arg2 ;
-  IDistribution1D *arg3 = 0 ;
-  void *argp1 = 0 ;
-  int res1 = 0 ;
-  int val2 ;
-  int ecode2 = 0 ;
-  void *argp3 = 0 ;
-  int res3 = 0 ;
-  PyObject *swig_obj[3] ;
-  
-  (void)self;
-  if (!SWIG_Python_UnpackTuple(args, "ISimulation_addParameterDistribution", 3, 3, swig_obj)) SWIG_fail;
-  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_ISimulation, 0 |  0 );
-  if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "ISimulation_addParameterDistribution" "', argument " "1"" of type '" "ISimulation *""'"); 
-  }
-  arg1 = reinterpret_cast< ISimulation * >(argp1);
-  ecode2 = SWIG_AsVal_int(swig_obj[1], &val2);
-  if (!SWIG_IsOK(ecode2)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "ISimulation_addParameterDistribution" "', argument " "2"" of type '" "ParameterDistribution::WhichParameter""'");
-  } 
-  arg2 = static_cast< ParameterDistribution::WhichParameter >(val2);
-  res3 = SWIG_ConvertPtr(swig_obj[2], &argp3, SWIGTYPE_p_IDistribution1D,  0  | 0);
-  if (!SWIG_IsOK(res3)) {
-    SWIG_exception_fail(SWIG_ArgError(res3), "in method '" "ISimulation_addParameterDistribution" "', argument " "3"" of type '" "IDistribution1D const &""'"); 
-  }
-  if (!argp3) {
-    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "ISimulation_addParameterDistribution" "', argument " "3"" of type '" "IDistribution1D const &""'"); 
-  }
-  arg3 = reinterpret_cast< IDistribution1D * >(argp3);
-  {
-    try {
-      (arg1)->addParameterDistribution(arg2,(IDistribution1D const &)*arg3);
-    } 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_Py_Void();
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
 SWIGINTERN PyObject *_wrap_ISimulation_setTerminalProgressMonitor(PyObject *self, PyObject *args) {
   PyObject *resultobj = 0;
   ISimulation *arg1 = (ISimulation *) 0 ;
@@ -31137,6 +31086,57 @@ fail:
 }
 
 
+SWIGINTERN PyObject *_wrap_ScatteringSimulation_addParameterDistribution(PyObject *self, PyObject *args) {
+  PyObject *resultobj = 0;
+  ScatteringSimulation *arg1 = (ScatteringSimulation *) 0 ;
+  ParameterDistribution::WhichParameter arg2 ;
+  IDistribution1D *arg3 = 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  int val2 ;
+  int ecode2 = 0 ;
+  void *argp3 = 0 ;
+  int res3 = 0 ;
+  PyObject *swig_obj[3] ;
+  
+  (void)self;
+  if (!SWIG_Python_UnpackTuple(args, "ScatteringSimulation_addParameterDistribution", 3, 3, swig_obj)) SWIG_fail;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_ScatteringSimulation, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "ScatteringSimulation_addParameterDistribution" "', argument " "1"" of type '" "ScatteringSimulation *""'"); 
+  }
+  arg1 = reinterpret_cast< ScatteringSimulation * >(argp1);
+  ecode2 = SWIG_AsVal_int(swig_obj[1], &val2);
+  if (!SWIG_IsOK(ecode2)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "ScatteringSimulation_addParameterDistribution" "', argument " "2"" of type '" "ParameterDistribution::WhichParameter""'");
+  } 
+  arg2 = static_cast< ParameterDistribution::WhichParameter >(val2);
+  res3 = SWIG_ConvertPtr(swig_obj[2], &argp3, SWIGTYPE_p_IDistribution1D,  0  | 0);
+  if (!SWIG_IsOK(res3)) {
+    SWIG_exception_fail(SWIG_ArgError(res3), "in method '" "ScatteringSimulation_addParameterDistribution" "', argument " "3"" of type '" "IDistribution1D const &""'"); 
+  }
+  if (!argp3) {
+    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "ScatteringSimulation_addParameterDistribution" "', argument " "3"" of type '" "IDistribution1D const &""'"); 
+  }
+  arg3 = reinterpret_cast< IDistribution1D * >(argp3);
+  {
+    try {
+      (arg1)->addParameterDistribution(arg2,(IDistribution1D const &)*arg3);
+    } 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_Py_Void();
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
 SWIGINTERN PyObject *ScatteringSimulation_swigregister(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *obj;
   if (!SWIG_Python_UnpackTuple(args, "swigregister", 1, 1, &obj)) return NULL;
@@ -35660,7 +35660,6 @@ static PyMethodDef SwigMethods[] = {
 	 { "QzScan_swiginit", QzScan_swiginit, METH_VARARGS, NULL},
 	 { "delete_ISimulation", _wrap_delete_ISimulation, METH_O, "delete_ISimulation(ISimulation self)"},
 	 { "ISimulation_setBackground", _wrap_ISimulation_setBackground, METH_VARARGS, "ISimulation_setBackground(ISimulation self, IBackground bg)"},
-	 { "ISimulation_addParameterDistribution", _wrap_ISimulation_addParameterDistribution, METH_VARARGS, "ISimulation_addParameterDistribution(ISimulation self, ParameterDistribution::WhichParameter whichParameter, IDistribution1D const & distribution)"},
 	 { "ISimulation_setTerminalProgressMonitor", _wrap_ISimulation_setTerminalProgressMonitor, METH_O, "ISimulation_setTerminalProgressMonitor(ISimulation self)"},
 	 { "ISimulation_options", _wrap_ISimulation_options, METH_O, "ISimulation_options(ISimulation self) -> SimulationOptions &"},
 	 { "ISimulation_simulate", _wrap_ISimulation_simulate, METH_O, "ISimulation_simulate(ISimulation self) -> Datafield"},
@@ -35670,6 +35669,7 @@ static PyMethodDef SwigMethods[] = {
 	 { "ScatteringSimulation_className", _wrap_ScatteringSimulation_className, METH_O, "ScatteringSimulation_className(ScatteringSimulation self) -> std::string"},
 	 { "ScatteringSimulation_beam", _wrap_ScatteringSimulation_beam, METH_O, "ScatteringSimulation_beam(ScatteringSimulation self) -> Beam &"},
 	 { "ScatteringSimulation_detector", _wrap_ScatteringSimulation_detector, METH_O, "ScatteringSimulation_detector(ScatteringSimulation self) -> IDetector &"},
+	 { "ScatteringSimulation_addParameterDistribution", _wrap_ScatteringSimulation_addParameterDistribution, METH_VARARGS, "ScatteringSimulation_addParameterDistribution(ScatteringSimulation self, ParameterDistribution::WhichParameter whichParameter, IDistribution1D const & distribution)"},
 	 { "ScatteringSimulation_swigregister", ScatteringSimulation_swigregister, METH_O, NULL},
 	 { "ScatteringSimulation_swiginit", ScatteringSimulation_swiginit, METH_VARARGS, NULL},
 	 { "new_DepthprobeSimulation", _wrap_new_DepthprobeSimulation, METH_VARARGS, "DepthprobeSimulation(BeamScan scan, Sample const & sample, Scale zaxis, int flags=0)"},
diff --git a/rawEx/varia/Resonator.py b/rawEx/varia/Resonator.py
index 239813eb81bd4e12204c847dfe5bc1d1e075a49b..5de1e42735dec9523a84e7ec694f23d95118d117 100755
--- a/rawEx/varia/Resonator.py
+++ b/rawEx/varia/Resonator.py
@@ -72,13 +72,12 @@ def get_simulation(sample):
     scan = ba.AlphaScan(na, ai_min, ai_max)
     scan.setWavelength(wl)
 
+    alpha_distr = ba.DistributionGaussian(0, d_ang, 25, 3.)
+    scan.setAngleDistribution(alpha_distr)
+
     z_axis = ba.EquiDivision("z (nm)", nz, z_min, z_max)
     simulation = ba.DepthprobeSimulation(scan, sample, z_axis)
 
-    alpha_distr = ba.DistributionGaussian(0, d_ang, 25, 3.)
-    simulation.addParameterDistribution(
-        ba.ParameterDistribution.BeamInclinationAngle, alpha_distr)
-
     return simulation