diff --git a/Core/Computation/ConstantBackground.cpp b/Core/Computation/ConstantBackground.cpp
index 22b34d0dcbefe714e313c15455465c662a426fa8..2e08f90a16184259884309cdd87bd17e09b0c1b3 100644
--- a/Core/Computation/ConstantBackground.cpp
+++ b/Core/Computation/ConstantBackground.cpp
@@ -14,7 +14,6 @@
 
 #include "ConstantBackground.h"
 #include "RealParameter.h"
-#include "SimulationElement.h"
 
 
 ConstantBackground::ConstantBackground(double background_value)
diff --git a/Core/Computation/IBackground.cpp b/Core/Computation/IBackground.cpp
index df3af9661cfb763c4e2734e43eacc0ed18c846e5..c11a562c678cce1c7b16db131f1184138126f5f3 100644
--- a/Core/Computation/IBackground.cpp
+++ b/Core/Computation/IBackground.cpp
@@ -13,6 +13,5 @@
 // ************************************************************************** //
 
 #include "IBackground.h"
-#include "SimulationElement.h"
 
 IBackground::~IBackground() = default;
diff --git a/Core/Computation/IBackground.h b/Core/Computation/IBackground.h
index dfe43cbd9a90f33a8a37f737a58fd7452a0b5676..4b784c5655ea33eb9e47e12e3eb6c16c64e1f80b 100644
--- a/Core/Computation/IBackground.h
+++ b/Core/Computation/IBackground.h
@@ -17,7 +17,6 @@
 
 #include "ICloneable.h"
 #include "INode.h"
-#include <vector>
 
 class SimulationElement;
 
diff --git a/Core/Computation/PoissonNoiseBackground.cpp b/Core/Computation/PoissonNoiseBackground.cpp
index 4dd3a9f531826045b6ab212d7d8a0619af308745..89b1adf2f36ecba393e859c8102925a9bc69162f 100644
--- a/Core/Computation/PoissonNoiseBackground.cpp
+++ b/Core/Computation/PoissonNoiseBackground.cpp
@@ -14,7 +14,6 @@
 
 #include "PoissonNoiseBackground.h"
 #include "MathFunctions.h"
-#include "SimulationElement.h"
 
 PoissonNoiseBackground::PoissonNoiseBackground()
 {
diff --git a/Core/Computation/SpecularComputation.cpp b/Core/Computation/SpecularComputation.cpp
index 0bff01d9911e07f57cd046e9e072070b1c001a20..e0a0f862d6b468094f64b047bf49d2807027207a 100644
--- a/Core/Computation/SpecularComputation.cpp
+++ b/Core/Computation/SpecularComputation.cpp
@@ -19,7 +19,6 @@
 #include "MultiLayer.h"
 #include "ScalarFresnelMap.h"
 #include "ProgressHandler.h"
-#include "SimulationElement.h"
 #include "SpecularComputationTerm.h"
 
 static_assert(std::is_copy_constructible<SpecularComputation>::value == false,
diff --git a/Core/Computation/SpecularComputationTerm.cpp b/Core/Computation/SpecularComputationTerm.cpp
index 7879085b61dffaa70f0ba8d3a02b5e6fc9121ab4..769309d198e740355ac4333c03180494542007cb 100644
--- a/Core/Computation/SpecularComputationTerm.cpp
+++ b/Core/Computation/SpecularComputationTerm.cpp
@@ -1,7 +1,6 @@
 #include "IFresnelMap.h"
 #include "ILayerRTCoefficients.h"
 #include "MultiLayer.h"
-#include "SimulationElement.h"
 #include "SpecularComputationTerm.h"
 #include "SpecularData.h"
 #include "SpecularSimulationElement.h"
@@ -12,18 +11,6 @@ SpecularComputationTerm::SpecularComputationTerm(const MultiLayer* p_multi_layer
     , mp_fresnel_map(p_fresnel_map)
 {}
 
-void SpecularComputationTerm::eval(ProgressHandler*,
-    const std::vector<SimulationElement>::iterator& begin_it,
-    const std::vector<SimulationElement>::iterator& end_it) const
-{
-    if (mp_multilayer->requiresMatrixRTCoefficients())
-        return;
-
-    for (auto it = begin_it; it != end_it; ++it)
-        if (it->specularData())
-            evalSingle(it);
-}
-
 void SpecularComputationTerm::eval(ProgressHandler*, const SpecularElementIter& begin_it,
                                    const SpecularElementIter& end_it) const
 {
@@ -34,13 +21,6 @@ void SpecularComputationTerm::eval(ProgressHandler*, const SpecularElementIter&
         evalSingle(it);
 }
 
-void SpecularComputationTerm::evalSingle(const std::vector<SimulationElement>::iterator& iter) const
-{
-    mp_fresnel_map->fillSpecularData(*iter);
-    const ILayerRTCoefficients& layer_data = (*iter->specularData())[0];
-    iter->setIntensity(std::norm(layer_data.getScalarR()));
-}
-
 void SpecularComputationTerm::evalSingle(const SpecularElementIter& iter) const
 {
     mp_fresnel_map->fillSpecularData(*iter);
diff --git a/Core/Computation/SpecularComputationTerm.h b/Core/Computation/SpecularComputationTerm.h
index 3f20a80be5548a3f6c453a8d56e46115db9cd349..e314868554eef5e84fbf3670a08eb51ed6726bf5 100644
--- a/Core/Computation/SpecularComputationTerm.h
+++ b/Core/Computation/SpecularComputationTerm.h
@@ -33,14 +33,10 @@ class SpecularComputationTerm
 public:
     SpecularComputationTerm(const MultiLayer* p_multi_layer, const IFresnelMap* p_fresnel_map);
 
-    void eval(ProgressHandler* progress, const std::vector<SimulationElement>::iterator& begin_it,
-              const std::vector<SimulationElement>::iterator& end_it) const;
-
     void eval(ProgressHandler* progress, const SpecularElementIter& begin_it,
               const SpecularElementIter& end_it) const;
 
 private:
-    void evalSingle(const std::vector<SimulationElement>::iterator& iter) const;
     void evalSingle(const SpecularElementIter& iter) const;
 
     const MultiLayer* mp_multilayer;
diff --git a/Core/Instrument/SpecularDetector1D.cpp b/Core/Instrument/SpecularDetector1D.cpp
index 0d822ed696954453051a0c6cf175a7d5794fded5..b7dbd5a14ba0929a59da0e95ead091a41e2a916c 100644
--- a/Core/Instrument/SpecularDetector1D.cpp
+++ b/Core/Instrument/SpecularDetector1D.cpp
@@ -17,7 +17,6 @@
 #include "IPixel.h"
 #include "OutputData.h"
 #include "SimulationArea.h"
-#include "SimulationElement.h"
 #include "SpecularDetector1D.h"
 #include "SpecularSimulationElement.h"
 #include "Units.h"
@@ -80,33 +79,8 @@ double SpecularDetector1D::alphaI(size_t index) const
     return axis.getBin(axis_index).getMidPoint();
 }
 
-std::vector<SimulationElement> SpecularDetector1D::createSimulationElements(const Beam& beam)
-{
-    std::vector<SimulationElement> result;
-
-    const double wavelength = beam.getWavelength();
-    const Eigen::Matrix2cd beam_polarization = beam.getPolarization();
-    const Eigen::Matrix2cd analyzer_operator = detectionProperties().analyzerOperator();
-    const double phi_i = 0; // Assuming that beam is always parallel to x-axis of the sample
-
-    SimulationArea area(this);
-    result.reserve(area.totalSize());
-    for (SimulationArea::iterator it = area.begin(); it != area.end(); ++it) {
-        // opposite sign for alpha_i since e_{z_beam} = - e_{z_detector}
-        const double alpha_i = -alphaI(it.detectorIndex());
-        result.emplace_back(wavelength, alpha_i, phi_i,
-                            std::make_unique<SpecularPixel>(alpha_i));
-        auto& sim_element = result.back();
-        sim_element.setPolarization(beam_polarization);
-        sim_element.setAnalyzerOperator(analyzer_operator);
-        sim_element.setSpecular();
-    }
-
-    return result;
-}
-
 std::vector<SpecularSimulationElement>
-SpecularDetector1D::createSpecularSimulationElements(const Beam& beam)
+SpecularDetector1D::createSimulationElements(const Beam& beam)
 {
     std::vector<SpecularSimulationElement> result;
 
diff --git a/Core/Instrument/SpecularDetector1D.h b/Core/Instrument/SpecularDetector1D.h
index 0c0fd34e87d3bf9b19c2e32000f720b1200f19db..6d7a04a148ff4e6d5c8514a35d188f14a95cea53 100644
--- a/Core/Instrument/SpecularDetector1D.h
+++ b/Core/Instrument/SpecularDetector1D.h
@@ -35,17 +35,13 @@ public:
     const DetectorMask* detectorMask() const override { return nullptr; }
 
 #ifndef SWIG
-    //! Create a vector of SimulationElement objects according to the detector
-    std::vector<SimulationElement> createSimulationElements(const Beam& beam);
-
     //! Create a vector of SpecularSimulationElement objects according to the detector
-    std::vector<SpecularSimulationElement> createSpecularSimulationElements(const Beam& beam);
+    std::vector<SpecularSimulationElement> createSimulationElements(const Beam& beam);
 
     //! Create a vector of DetectorElement objects according to the detector
     std::vector<DetectorElement> createDetectorElements(const Beam& beam) override;
 #endif // SWIG
 
-    using IDetector::createDetectorIntensity;
     //! Returns new intensity map with detector resolution applied and axes in requested units
     OutputData<double>*
     createDetectorIntensity(const std::vector<SpecularSimulationElement>& elements,
diff --git a/Core/Multilayer/IFresnelMap.h b/Core/Multilayer/IFresnelMap.h
index 46642dfe53c587022732e4c7e87e60b100c079c2..95c932276b0f70b8a3e74271aa6bd871d069fa37 100644
--- a/Core/Multilayer/IFresnelMap.h
+++ b/Core/Multilayer/IFresnelMap.h
@@ -44,7 +44,6 @@ public:
             const SimulationElement& sim_element, size_t layer_index) const =0;
 
     //! Fills simulation element specular data
-    virtual void fillSpecularData(SimulationElement& sim_element) const = 0;
     virtual void fillSpecularData(SpecularSimulationElement& sim_element) const = 0;
 
     //! Sets the multilayer to be used for the Fresnel calculations.
diff --git a/Core/Multilayer/MatrixFresnelMap.cpp b/Core/Multilayer/MatrixFresnelMap.cpp
index 37fc2181bb81ad3274bf2eead1b72e5d0f21b8da..c3bebcd308e85e0dc788dfda4c483291d0d4811f 100644
--- a/Core/Multilayer/MatrixFresnelMap.cpp
+++ b/Core/Multilayer/MatrixFresnelMap.cpp
@@ -47,17 +47,6 @@ MatrixFresnelMap::getInCoefficients(const SimulationElement& sim_element, size_t
     return getCoefficients(sim_element.getKi(), layer_index, *mP_multilayer, m_hash_table_in);
 }
 
-void MatrixFresnelMap::fillSpecularData(SimulationElement& sim_element) const
-{
-    const auto& kvec = sim_element.getKi();
-    std::vector<MatrixRTCoefficients> coef_vector;
-    if (m_use_cache)
-        coef_vector = getCoefficientsFromCache(kvec, *mP_multilayer, m_hash_table_in);
-    else
-        coef_vector = calculateCoefficients(*mP_multilayer, kvec);
-    sim_element.setSpecular(std::make_unique<SpecularData>(std::move(coef_vector)));
-}
-
 void MatrixFresnelMap::fillSpecularData(SpecularSimulationElement& sim_element) const
 {
     const auto& kvec = sim_element.getKi();
diff --git a/Core/Multilayer/MatrixFresnelMap.h b/Core/Multilayer/MatrixFresnelMap.h
index 32ee7fc63262f5e5f68f1348e6625d063572e95e..b02c461588c5e07bfc49648cf6357513907f861e 100644
--- a/Core/Multilayer/MatrixFresnelMap.h
+++ b/Core/Multilayer/MatrixFresnelMap.h
@@ -44,7 +44,6 @@ public:
     void setMultilayer(const MultiLayer& multilayer) final override;
 
     //! Fills simulation element specular data
-    void fillSpecularData(SimulationElement& sim_element) const override;
     void fillSpecularData(SpecularSimulationElement& sim_element) const override;
 
     typedef std::unordered_map<kvector_t, std::vector<MatrixRTCoefficients>, HashKVector>
diff --git a/Core/Multilayer/ScalarFresnelMap.cpp b/Core/Multilayer/ScalarFresnelMap.cpp
index f3e217cc44a7f53a66a3b664191caf0680dfd487..e15740d61a34c8edebe69d9fa372a26c9750bdf9 100644
--- a/Core/Multilayer/ScalarFresnelMap.cpp
+++ b/Core/Multilayer/ScalarFresnelMap.cpp
@@ -43,17 +43,6 @@ const ILayerRTCoefficients* ScalarFresnelMap::getInCoefficients(
     return getCoefficients(sim_element.getKi(), layer_index);
 }
 
-void ScalarFresnelMap::fillSpecularData(SimulationElement& sim_element) const
-{
-    const auto& kvec = sim_element.getKi();
-    std::vector<ScalarRTCoefficients> coef_vector;
-    if (m_use_cache)
-        coef_vector = getCoefficientsFromCache(kvec);
-    else
-        coef_vector = calculateCoefficients(*mP_multilayer, kvec);
-    sim_element.setSpecular(std::make_unique<SpecularData>(std::move(coef_vector)));
-}
-
 void ScalarFresnelMap::fillSpecularData(SpecularSimulationElement& sim_element) const
 {
     const auto& kvec = sim_element.getKi();
diff --git a/Core/Multilayer/ScalarFresnelMap.h b/Core/Multilayer/ScalarFresnelMap.h
index bd0b7fcabcab81fb3642a6d1577b4943a8d0f987..7ebd52cc264e2d928078899bb04232b5310d073a 100644
--- a/Core/Multilayer/ScalarFresnelMap.h
+++ b/Core/Multilayer/ScalarFresnelMap.h
@@ -44,7 +44,6 @@ public:
                                                   size_t layer_index) const final override;
 
     //! Fills simulation element specular data
-    void fillSpecularData(SimulationElement& sim_element) const override;
     void fillSpecularData(SpecularSimulationElement& sim_element) const override;
 
 private:
diff --git a/Core/Simulation/SpecularSimulation.cpp b/Core/Simulation/SpecularSimulation.cpp
index 94295a02ca6f5929ca2c246466a2ba0b9764934b..f2d438f88bbbe24036453cc242fb27b39fd8431b 100644
--- a/Core/Simulation/SpecularSimulation.cpp
+++ b/Core/Simulation/SpecularSimulation.cpp
@@ -17,10 +17,8 @@
 #include "IFootprintFactor.h"
 #include "IMultiLayerBuilder.h"
 #include "MultiLayer.h"
-#include "SpecularMatrix.h"
 #include "MaterialUtils.h"
 #include "Histogram1D.h"
-#include "SimulationElement.h"
 #include "SpecularComputation.h"
 #include "SpecularData.h"
 #include "SpecularDetector1D.h"
@@ -112,7 +110,7 @@ void SpecularSimulation::initSimulationElementVector()
 std::vector<SpecularSimulationElement> SpecularSimulation::generateSimulationElements(const Beam& beam)
 {
     auto p_detector = SpecDetector(m_instrument);
-    return p_detector->createSpecularSimulationElements(beam);
+    return p_detector->createSimulationElements(beam);
 }
 
 std::vector<complex_t> SpecularSimulation::getData(size_t i_layer, DataGetter fn_ptr) const
diff --git a/Tests/UnitTests/Core/Detector/SpecularDetector1DTest.h b/Tests/UnitTests/Core/Detector/SpecularDetector1DTest.h
index 117ae5da7747ae2e83dfabfe725c8daf1c43dd46..d992346b8c9227d6f610d693667879099e6fbfda 100644
--- a/Tests/UnitTests/Core/Detector/SpecularDetector1DTest.h
+++ b/Tests/UnitTests/Core/Detector/SpecularDetector1DTest.h
@@ -2,9 +2,9 @@
 #include "BornAgainNamespace.h"
 #include "FixedBinAxis.h"
 #include "OutputData.h"
-#include "SimulationElement.h"
 #include "SimulationArea.h"
 #include "SpecularDetector1D.h"
+#include "SpecularSimulationElement.h"
 #include "Units.h"
 #include "google_test.h"
 #include <memory>
@@ -116,20 +116,10 @@ TEST_F(SpecularDetectorTest, SimulationElements)
     EXPECT_EQ(5u, sim_elements.size());
 
     EXPECT_NEAR(axis.getBinCenter(0), -sim_elements[0].getAlphaI(), 1e-10);
-    EXPECT_NEAR(0.0, sim_elements[0].getPhiI(), 1e-10);
     EXPECT_NEAR(beam.getWavelength(), sim_elements[0].getWavelength(), 1e-10);
-    EXPECT_NEAR(axis.getBinCenter(0), sim_elements[0].getAlphaMean(), 1e-10);
-    EXPECT_NEAR(0.0, sim_elements[0].getPhiMean(), 1e-10);
-    EXPECT_EQ(sim_elements[0].getAlphaMean(), sim_elements[0].getAlpha(0, 0));
-    EXPECT_EQ(sim_elements[0].getAlphaMean(), sim_elements[0].getAlpha(1, 1));
 
     EXPECT_NEAR(axis.getBinCenter(4), -sim_elements[4].getAlphaI(), 1e-10);
-    EXPECT_NEAR(0.0, sim_elements[4].getPhiI(), 1e-10);
     EXPECT_NEAR(beam.getWavelength(), sim_elements[4].getWavelength(), 1e-10);
-    EXPECT_NEAR(axis.getBinCenter(4), sim_elements[4].getAlphaMean(), 1e-10);
-    EXPECT_NEAR(0.0, sim_elements[4].getPhiMean(), 1e-10);
-    EXPECT_EQ(sim_elements[4].getAlphaMean(), sim_elements[4].getAlpha(0, 0));
-    EXPECT_EQ(sim_elements[4].getAlphaMean(), sim_elements[4].getAlpha(1, 1));
 }
 
 TEST_F(SpecularDetectorTest, Clone)