diff --git a/Examples/scatter2d/MagneticSpheres.py b/Examples/scatter2d/MagneticSpheres.py
index c1483742f7158a1f97c9f3227a6fb004244c838a..04a0ad72f272ddd72f3cd1985bd2acfdea0ab5ce 100755
--- a/Examples/scatter2d/MagneticSpheres.py
+++ b/Examples/scatter2d/MagneticSpheres.py
@@ -53,7 +53,7 @@ def get_simulation(sample):
 
     polarizer_dir = R3(0, 0, 1)
     analyzer_dir = R3(0, 0, -1)
-    simulation.beam().setPolarization(polarizer_dir)
+    simulation.setPolarization(polarizer_dir)
     simulation.detector().setAnalyzer(analyzer_dir, 1, 0.5)
 
     return simulation
diff --git a/Examples/specular/BasicPolarizedReflectometry.py b/Examples/specular/BasicPolarizedReflectometry.py
index c45f2dae4e46ee68cd82b0caf1dc836556ce9cc1..d733db722008cf994e309b717cd2d427fb79bddc 100755
--- a/Examples/specular/BasicPolarizedReflectometry.py
+++ b/Examples/specular/BasicPolarizedReflectometry.py
@@ -49,7 +49,7 @@ def simulate(polarizer_dir, analyzer_dir, title):
     simulation.setScan(scan)
     simulation.setSample(get_sample())
 
-    simulation.beam().setPolarization(polarizer_dir)
+    simulation.setPolarization(polarizer_dir)
     simulation.detector().setAnalyzer(analyzer_dir, 1, 0.5)
 
     result = simulation.simulate()
diff --git a/Examples/specular/PolarizedNoAnalyzer.py b/Examples/specular/PolarizedNoAnalyzer.py
index a3310d890a7a6a9e145d5bd05f45d1767d126da6..b1d51f7bdf5e420650d99e25ede793abf12b79e3 100755
--- a/Examples/specular/PolarizedNoAnalyzer.py
+++ b/Examples/specular/PolarizedNoAnalyzer.py
@@ -56,7 +56,7 @@ def run_simulation(polarizer_dir=R3(0, 1, 0), analyzer_dir=None):
     simulation = get_simulation(sample)
 
     # adding polarizer and analyzer operator
-    simulation.beam().setPolarization(polarizer_dir)
+    simulation.setPolarization(polarizer_dir)
     if analyzer_dir:
         simulation.detector().setAnalyzer(analyzer_dir, 1, 0.5)
 
diff --git a/Examples/specular/PolarizedNonperfectAnalyzerPolarizer.py b/Examples/specular/PolarizedNonperfectAnalyzerPolarizer.py
index 86ae7c0fa68e120c4934be55c9e8cff75da670d2..15fbf3791a0ec394236730172ca3d06716be9ce6 100755
--- a/Examples/specular/PolarizedNonperfectAnalyzerPolarizer.py
+++ b/Examples/specular/PolarizedNonperfectAnalyzerPolarizer.py
@@ -80,7 +80,7 @@ def run_simulation(*,
     sample = get_sample()
     simulation = get_simulation(sample)
 
-    simulation.beam().setPolarization(polarizer_dir*polarizer_efficiency)
+    simulation.setPolarization(polarizer_dir*polarizer_efficiency)
     simulation.detector().setAnalyzer(analyzer_dir, analyzer_efficiency, 0.5)
     simulation.setBackground(ba.ConstantBackground(1e-7))
 
diff --git a/Examples/specular/PolarizedSpinAsymmetry.py b/Examples/specular/PolarizedSpinAsymmetry.py
index aa40ca18d42c665f4b14fc17f29945433cda54a5..10f6cb539570e6818001bbe6028ef7796347583c 100755
--- a/Examples/specular/PolarizedSpinAsymmetry.py
+++ b/Examples/specular/PolarizedSpinAsymmetry.py
@@ -80,7 +80,7 @@ def get_simulation(q_axis, parameters, polarizer_dir, analyzer_dir):
     distr = ba.RangedDistributionGaussian(n_samples, n_sig)
     scan.setAbsoluteQResolution(distr, parameters["q_res"])
 
-    simulation.beam().setPolarization(polarizer_dir)
+    simulation.setPolarization(polarizer_dir)
     simulation.detector().setAnalyzer(analyzer_dir, 1, 0.5)
 
     simulation.setScan(scan)
diff --git a/Examples/specular/PolarizedSpinFlip.py b/Examples/specular/PolarizedSpinFlip.py
index 3bc5cdb4704e744080f0b455efe5a62b42b8f1a9..109088f228f4aca394c2a2e4a083b7ddd641f4d4 100755
--- a/Examples/specular/PolarizedSpinFlip.py
+++ b/Examples/specular/PolarizedSpinFlip.py
@@ -55,7 +55,7 @@ def run_simulation(polarizer_dir=ba.R3(0, 1, 0),
     sample = get_sample()
     simulation = get_simulation(sample)
 
-    simulation.beam().setPolarization(polarizer_dir)
+    simulation.setPolarization(polarizer_dir)
     simulation.detector().setAnalyzer(analyzer_dir, 1, 0.5)
 
     return simulation.simulate()
diff --git a/Sim/Simulation/ISimulation.cpp b/Sim/Simulation/ISimulation.cpp
index aa171fce490e0f94fa16c4ce2746facceb71954d..d85c8b41692143040601f2ed9e2e649c79c3dae2 100644
--- a/Sim/Simulation/ISimulation.cpp
+++ b/Sim/Simulation/ISimulation.cpp
@@ -174,6 +174,11 @@ void ISimulation::setTerminalProgressMonitor()
     });
 }
 
+void ISimulation::setPolarization(R3 bloch_vector)
+{
+    beam().setPolarization(bloch_vector);
+}
+
 void ISimulation::prepareSimulation()
 {
     if (!m_sample)
diff --git a/Sim/Simulation/ISimulation.h b/Sim/Simulation/ISimulation.h
index 582309d26141bd5d6940c7119ebf84cae9c65640..ce037f41f81f08a57dc07e0e591d6769aa3c064c 100644
--- a/Sim/Simulation/ISimulation.h
+++ b/Sim/Simulation/ISimulation.h
@@ -17,6 +17,7 @@
 
 #include "Param/Distrib/DistributionHandler.h"
 #include "Param/Node/INode.h"
+#include <heinz/Vectors3D.h>
 
 template <class T>
 class Powerfield;
@@ -71,6 +72,9 @@ public:
     Beam& beam() { return *m_beam; }
     IDetector& detector() { return *m_detector; }
 
+    //! Sets the polarization density matrix according to the given Bloch vector
+    void setPolarization(R3 bloch_vector);
+
 #ifndef SWIG
     void subscribe(const std::function<bool(size_t)>& inform);
 
diff --git a/Tests/Functional/Suite/MakeSimulations.cpp b/Tests/Functional/Suite/MakeSimulations.cpp
index 3222c21f8a33cc67ccade05462558a90efbb9bfb..592f7374693d959b4af70b47b13a27ef03717ff2 100644
--- a/Tests/Functional/Suite/MakeSimulations.cpp
+++ b/Tests/Functional/Suite/MakeSimulations.cpp
@@ -111,7 +111,7 @@ std::unique_ptr<ScatteringSimulation> test::makeSimulation::BasicGISAS(const Mul
 std::unique_ptr<ScatteringSimulation> test::makeSimulation::SpinflipGISAS(const MultiLayer& sample)
 {
     std::unique_ptr<ScatteringSimulation> result = BasicGISAS(sample);
-    result->beam().setPolarization(zplus);
+    result->setPolarization(zplus);
     result->detector().setAnalyzer(zplus, -1.0, 0.5);
     return result;
 }
@@ -204,7 +204,7 @@ std::unique_ptr<ScatteringSimulation> test::makeSimulation::MaxiGISAS(const Mult
 std::unique_ptr<ScatteringSimulation> test::makeSimulation::MaxiGISAS00(const MultiLayer& sample)
 {
     std::unique_ptr<ScatteringSimulation> result = MaxiGISAS(sample);
-    result->beam().setPolarization(zplus);
+    result->setPolarization(zplus);
     result->detector().setAnalyzer(zplus, 1.0, 0.5);
     return result;
 }
@@ -391,7 +391,7 @@ test::makeSimulation::BasicYPolarizedSpecular(const MultiLayer& sample, const st
 {
     const auto yCase = YPolarizationCases.at(polCase);
     auto simulation = BasicSpecular(sample, vsQ);
-    simulation->beam().setPolarization(yCase.first);
+    simulation->setPolarization(yCase.first);
     simulation->detector().setAnalyzer(yCase.second, 1.0, 0.5);
     simulation->setSample(sample);
     return simulation;
diff --git a/Tests/Unit/Sim/SpecularSimulationTest.cpp b/Tests/Unit/Sim/SpecularSimulationTest.cpp
index 4e007b1c2c26c4ec7807f93286dce9768aaf3f78..0fdaad9a2905edd49fc4181fbb84f4628412f2b5 100644
--- a/Tests/Unit/Sim/SpecularSimulationTest.cpp
+++ b/Tests/Unit/Sim/SpecularSimulationTest.cpp
@@ -91,7 +91,7 @@ TEST_F(SpecularSimulationTest, SetAngularScan)
     AlphaScan scan4(1.0, 10, .0 * Units::deg, 2.0 * Units::deg);
     const auto polarizer_dir = R3({0., 0., 0.876});
     const auto analyzer_dir = R3({0., 0., 1.});
-    sim4.beam().setPolarization(polarizer_dir);
+    sim4.setPolarization(polarizer_dir);
     sim4.detector().setAnalyzer(analyzer_dir, 0.33, 0.22);
     sim4.setScan(scan4);
     EXPECT_FAILED_ASSERT(sim4.setScan(scan4));
@@ -145,7 +145,7 @@ TEST_F(SpecularSimulationTest, SetQScan)
     QzScan scan3(10, 1.0, 10.0);
     const auto polarizer_dir = R3({0., 0., 0.876});
     const auto analyzer_dir = R3({0., 0., 1.});
-    sim3.beam().setPolarization(polarizer_dir);
+    sim3.setPolarization(polarizer_dir);
     sim3.detector().setAnalyzer(analyzer_dir, 0.33, 0.22);
     sim3.setScan(scan3);
 
diff --git a/auto/Wrap/doxygenDevice.i b/auto/Wrap/doxygenDevice.i
index 5a1f2119b5171bad4b28a9f38ebd7d54a992b10a..d62d9b57edaa2f71f787b82e04c4d0207deee6d5 100644
--- a/auto/Wrap/doxygenDevice.i
+++ b/auto/Wrap/doxygenDevice.i
@@ -1499,56 +1499,6 @@ Returns true if area defined by two bins is inside or on border of polygon (more
 ";
 
 
-// File: classInstrument.xml
-%feature("docstring") Instrument "
-
-Assembles beam, detector and their relative positions with respect to the sample.
-
-C++ includes: Instrument.h
-";
-
-%feature("docstring")  Instrument::Instrument "Instrument::Instrument()
-";
-
-%feature("docstring")  Instrument::Instrument "Instrument::Instrument(const Beam &beam, const IDetector &detector)
-";
-
-%feature("docstring")  Instrument::Instrument "Instrument::Instrument(const Instrument &other)
-";
-
-%feature("docstring")  Instrument::~Instrument "Instrument::~Instrument() override
-";
-
-%feature("docstring")  Instrument::className "std::string Instrument::className() const final
-";
-
-%feature("docstring")  Instrument::beam "Beam& Instrument::beam()
-";
-
-%feature("docstring")  Instrument::beam "const Beam& Instrument::beam() const
-";
-
-%feature("docstring")  Instrument::getDetector "IDetector * Instrument::getDetector()
-";
-
-%feature("docstring")  Instrument::getDetector "const IDetector * Instrument::getDetector() const
-";
-
-%feature("docstring")  Instrument::detector "IDetector & Instrument::detector()
-";
-
-%feature("docstring")  Instrument::detector "const IDetector & Instrument::detector() const
-";
-
-%feature("docstring")  Instrument::setDetector "void Instrument::setDetector(const IDetector &detector)
-
-Sets the detector (axes can be overwritten later) 
-";
-
-%feature("docstring")  Instrument::nodeChildren "std::vector< const INode * > Instrument::nodeChildren() const override
-";
-
-
 // File: classIOFactory.xml
 %feature("docstring") IOFactory "
 
@@ -2989,12 +2939,6 @@ make Swappable
 // File: ReadWriteTiff_8h.xml
 
 
-// File: Instrument_8cpp.xml
-
-
-// File: Instrument_8h.xml
-
-
 // File: DetectorMask_8cpp.xml
 
 
@@ -3085,9 +3029,6 @@ make Swappable
 // File: dir_9f013251ba980bff6504d6613b69183d.xml
 
 
-// File: dir_550e786a97bd4c801929243ea9773c04.xml
-
-
 // File: dir_4866552d576e04b61ad8ade47c8db877.xml
 
 
diff --git a/auto/Wrap/doxygenSim.i b/auto/Wrap/doxygenSim.i
index 0820f5320a9e14552ea39fe6c6fc8a3e987f617d..309a91d1d51e41983fa263478b21213d89ba7425 100644
--- a/auto/Wrap/doxygenSim.i
+++ b/auto/Wrap/doxygenSim.i
@@ -910,6 +910,11 @@ Initializes a progress monitor that prints to stdout.
 %feature("docstring")  ISimulation::detector "IDetector& ISimulation::detector()
 ";
 
+%feature("docstring")  ISimulation::setPolarization "void ISimulation::setPolarization(R3 bloch_vector)
+
+Sets the polarization density matrix according to the given Bloch vector. 
+";
+
 %feature("docstring")  ISimulation::subscribe "void ISimulation::subscribe(const std::function< bool(size_t)> &inform)
 ";
 
diff --git a/auto/Wrap/libBornAgainSim.py b/auto/Wrap/libBornAgainSim.py
index ce499f33a2cebe0c16fb9eb6e596af425cb4c3db..9050961aeb285441c5e971d7b09e0c51efc5d488 100644
--- a/auto/Wrap/libBornAgainSim.py
+++ b/auto/Wrap/libBornAgainSim.py
@@ -3422,6 +3422,16 @@ class ISimulation(libBornAgainParam.INode):
         """
         return _libBornAgainSim.ISimulation_detector(self)
 
+    def setPolarization(self, bloch_vector):
+        r"""
+        setPolarization(ISimulation self, R3 bloch_vector)
+        void ISimulation::setPolarization(R3 bloch_vector)
+
+        Sets the polarization density matrix according to the given Bloch vector. 
+
+        """
+        return _libBornAgainSim.ISimulation_setPolarization(self, bloch_vector)
+
 # Register ISimulation in _libBornAgainSim:
 _libBornAgainSim.ISimulation_swigregister(ISimulation)
 
diff --git a/auto/Wrap/libBornAgainSim_wrap.cpp b/auto/Wrap/libBornAgainSim_wrap.cpp
index 66ab5bbcf56b908761f515ba007c23f1b5631cb5..cb2aa8e42160cd261fc97b62e8ee49dd58ce9d49 100644
--- a/auto/Wrap/libBornAgainSim_wrap.cpp
+++ b/auto/Wrap/libBornAgainSim_wrap.cpp
@@ -37707,6 +37707,43 @@ fail:
 }
 
 
+SWIGINTERN PyObject *_wrap_ISimulation_setPolarization(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+  PyObject *resultobj = 0;
+  ISimulation *arg1 = (ISimulation *) 0 ;
+  R3 arg2 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  void *argp2 ;
+  int res2 = 0 ;
+  PyObject *swig_obj[2] ;
+  
+  if (!SWIG_Python_UnpackTuple(args, "ISimulation_setPolarization", 2, 2, 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_setPolarization" "', argument " "1"" of type '" "ISimulation *""'"); 
+  }
+  arg1 = reinterpret_cast< ISimulation * >(argp1);
+  {
+    res2 = SWIG_ConvertPtr(swig_obj[1], &argp2, SWIGTYPE_p_Vec3T_double_t,  0  | 0);
+    if (!SWIG_IsOK(res2)) {
+      SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "ISimulation_setPolarization" "', argument " "2"" of type '" "R3""'"); 
+    }  
+    if (!argp2) {
+      SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "ISimulation_setPolarization" "', argument " "2"" of type '" "R3""'");
+    } else {
+      R3 * temp = reinterpret_cast< R3 * >(argp2);
+      arg2 = *temp;
+      if (SWIG_IsNewObj(res2)) delete temp;
+    }
+  }
+  (arg1)->setPolarization(arg2);
+  resultobj = SWIG_Py_Void();
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
 SWIGINTERN PyObject *ISimulation_swigregister(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *obj;
   if (!SWIG_Python_UnpackTuple(args, "swigregister", 1, 1, &obj)) return NULL;
@@ -41850,6 +41887,13 @@ static PyMethodDef SwigMethods[] = {
 		"const IDetector& ISimulation::detector() const\n"
 		"\n"
 		""},
+	 { "ISimulation_setPolarization", _wrap_ISimulation_setPolarization, METH_VARARGS, "\n"
+		"ISimulation_setPolarization(ISimulation self, R3 bloch_vector)\n"
+		"void ISimulation::setPolarization(R3 bloch_vector)\n"
+		"\n"
+		"Sets the polarization density matrix according to the given Bloch vector. \n"
+		"\n"
+		""},
 	 { "ISimulation_swigregister", ISimulation_swigregister, METH_O, NULL},
 	 { "delete_ISimulation2D", _wrap_delete_ISimulation2D, METH_O, "\n"
 		"delete_ISimulation2D(ISimulation2D self)\n"