diff --git a/Core/Examples/StandardSimulations.cpp b/Core/Examples/StandardSimulations.cpp
index 78e32e93fa9cceff0191f462fc5f4e1a1ff541db..ea69b50a350950d8f88d34fe8abdcfc9b184b691 100644
--- a/Core/Examples/StandardSimulations.cpp
+++ b/Core/Examples/StandardSimulations.cpp
@@ -160,9 +160,9 @@ GISASSimulation* StandardSimulations::MiniGISASPolarizationPP()
     GISASSimulation* result = MiniGISAS();
 
     kvector_t analyzer_dir(0.0, 0.0, 1.0);
-    kvector_t beampol(0.0, 0.0, 1.0);
+    kvector_t polarizer_dir(0.0, 0.0, 1.0);
 
-    result->beam().setPolarization(beampol);
+    result->beam().setPolarization(polarizer_dir);
     result->detector().setAnalyzer(analyzer_dir, 1.0, 0.5);
     return result;
 }
@@ -172,9 +172,9 @@ GISASSimulation* StandardSimulations::MiniGISASPolarizationPM()
     GISASSimulation* result = MiniGISAS();
 
     kvector_t analyzer_dir(0.0, 0.0, -1.0);
-    kvector_t beampol(0.0, 0.0, 1.0);
+    kvector_t polarizer_dir(0.0, 0.0, 1.0);
 
-    result->beam().setPolarization(beampol);
+    result->beam().setPolarization(polarizer_dir);
     result->detector().setAnalyzer(analyzer_dir, 1.0, 0.5);
     return result;
 }
@@ -184,9 +184,9 @@ GISASSimulation* StandardSimulations::MiniGISASPolarizationMP()
     GISASSimulation* result = MiniGISAS();
 
     kvector_t analyzer_dir(0.0, 0.0, 1.0);
-    kvector_t beampol(0.0, 0.0, -1.0);
+    kvector_t polarizer_dir(0.0, 0.0, -1.0);
 
-    result->beam().setPolarization(beampol);
+    result->beam().setPolarization(polarizer_dir);
     result->detector().setAnalyzer(analyzer_dir, 1.0, 0.5);
     return result;
 }
@@ -196,9 +196,9 @@ GISASSimulation* StandardSimulations::MiniGISASPolarizationMM()
     GISASSimulation* result = MiniGISAS();
 
     kvector_t analyzer_dir(0.0, 0.0, -1.0);
-    kvector_t beampol(0.0, 0.0, -1.0);
+    kvector_t polarizer_dir(0.0, 0.0, -1.0);
 
-    result->beam().setPolarization(beampol);
+    result->beam().setPolarization(polarizer_dir);
     result->detector().setAnalyzer(analyzer_dir, 1.0, 0.5);
     return result;
 }
diff --git a/Examples/fit55_SpecularIntro/PolarizedSpinAsymmetryFit.py b/Examples/fit55_SpecularIntro/PolarizedSpinAsymmetryFit.py
index 93b0805fccf80d65d69a0a500bd2b59b6980de14..7a97e3dd79d086b33aeacc9c114c8818550b29f3 100755
--- a/Examples/fit55_SpecularIntro/PolarizedSpinAsymmetryFit.py
+++ b/Examples/fit55_SpecularIntro/PolarizedSpinAsymmetryFit.py
@@ -64,7 +64,7 @@ def get_sample(params):
     return multi_layer
 
 
-def get_simulation(q_axis, parameters, polarization, analyzer):
+def get_simulation(q_axis, parameters, polarizer_dir, analyzer_dir):
     """
     Returns a simulation object.
     Polarization, analyzer and resolution are set
@@ -81,14 +81,14 @@ def get_simulation(q_axis, parameters, polarization, analyzer):
     distr = ba.RangedDistributionGaussian(n_samples, n_sig)
     scan.setAbsoluteQResolution(distr, parameters["q_res"])
 
-    simulation.beam().setPolarization(polarization)
-    simulation.detector().setAnalyzer(analyzer, 1, 0.5)
+    simulation.beam().setPolarization(polarizer_dir)
+    simulation.detector().setAnalyzer(analyzer_dir, 1, 0.5)
 
     simulation.setScan(scan)
     return simulation
 
 
-def run_simulation(q_axis, fitParams, *, polarization, analyzer):
+def run_simulation(q_axis, fitParams, *, polarizer_dir, analyzer_dir):
     """
     Run a simulation on the given q-axis, where the sample is constructed
     with the given parameters.
@@ -97,7 +97,7 @@ def run_simulation(q_axis, fitParams, *, polarization, analyzer):
     parameters = dict(fitParams, **fixedParams)
 
     sample = get_sample(parameters)
-    simulation = get_simulation(q_axis, parameters, polarization, analyzer)
+    simulation = get_simulation(q_axis, parameters, polarizer_dir, analyzer_dir)
 
     simulation.setSample(sample)
     simulation.runSimulation()
diff --git a/Examples/scatter2d/MagneticSpheres.py b/Examples/scatter2d/MagneticSpheres.py
index bbb544b11112e75fdfec976cf0a0fd9a4de51b49..5d43278af3538c8ec22c03159c8dfad4c41ba72e 100755
--- a/Examples/scatter2d/MagneticSpheres.py
+++ b/Examples/scatter2d/MagneticSpheres.py
@@ -48,12 +48,12 @@ def get_sample():
 
 def get_simulation(sample):
     beam = ba.Beam(1e12, 0.1*nm, ba.Direction(0.5*deg, 0))
-    beam_polpair = kvector_t(0, 0, 1)
-    beam.setPolarization(beam_polpair)
+    polarizer_dir = kvector_t(0, 0, 1)
+    beam.setPolarization(polarizer_dir)
     detector = ba.SphericalDetector(200, 6*deg, 0, 3*deg)
     simulation = ba.GISASSimulation(beam, sample, detector)
-    analyzer_direction = kvector_t(0, 0, -1)
-    simulation.detector().setAnalyzer(analyzer_direction, 1, 0.5)
+    analyzer_dir = kvector_t(0, 0, -1)
+    simulation.detector().setAnalyzer(analyzer_dir, 1, 0.5)
     return simulation
 
 
diff --git a/Examples/varia/BasicPolarizedReflectometry.py b/Examples/varia/BasicPolarizedReflectometry.py
index 4e5af3c774e4541fc90a64715e6ef2334f225ef1..6e9c278727245e7a43e4f50d62293e98083ec0ae 100755
--- a/Examples/varia/BasicPolarizedReflectometry.py
+++ b/Examples/varia/BasicPolarizedReflectometry.py
@@ -47,8 +47,8 @@ def get_simulation(sample, scan_size=500):
     return simulation
 
 
-def run_simulation(polarization=ba.kvector_t(0, 1, 0),
-                   analyzer=ba.kvector_t(0, 1, 0)):
+def run_simulation(polarizer_dir=ba.kvector_t(0, 1, 0),
+                   analyzer_dir=ba.kvector_t(0, 1, 0)):
     """
     Runs simulation and returns its result.
     """
@@ -56,8 +56,8 @@ def run_simulation(polarization=ba.kvector_t(0, 1, 0),
     simulation = get_simulation(sample)
 
     # adding polarization and analyzer operator
-    simulation.beam().setPolarization(polarization)
-    simulation.detector().setAnalyzer(analyzer, 1, 0.5)
+    simulation.beam().setPolarization(polarizer_dir)
+    simulation.detector().setAnalyzer(analyzer_dir, 1, 0.5)
 
     simulation.runSimulation()
     return simulation.result()
diff --git a/Examples/varia/PolarizedNoAnalyzer.py b/Examples/varia/PolarizedNoAnalyzer.py
index 95d25769a07d865ec20813a2dd25b100799f864c..3d69f71350b3a70ca0aa74287a21d4326f8c374f 100755
--- a/Examples/varia/PolarizedNoAnalyzer.py
+++ b/Examples/varia/PolarizedNoAnalyzer.py
@@ -48,17 +48,17 @@ def get_simulation(sample, scan_size=500):
     return simulation
 
 
-def run_simulation(polarization=ba.kvector_t(0, 1, 0), analyzer=None):
+def run_simulation(polarizer_dir=ba.kvector_t(0, 1, 0), analyzer_dir=None):
     """
     Runs simulation and returns its result.
     """
     sample = get_sample()
     simulation = get_simulation(sample)
 
-    # adding polarization and analyzer operator
-    simulation.beam().setPolarization(polarization)
-    if analyzer:
-        simulation.detector().setAnalyzer(analyzer, 1, 0.5)
+    # adding polarizer and analyzer operator
+    simulation.beam().setPolarization(polarizer_dir)
+    if analyzer_dir:
+        simulation.detector().setAnalyzer(analyzer_dir, 1, 0.5)
 
     simulation.runSimulation()
     result = simulation.result()
diff --git a/Examples/varia/PolarizedNonperfectAnalyzerPolarizer.py b/Examples/varia/PolarizedNonperfectAnalyzerPolarizer.py
index 023e8f9c826041cacc4c646984fbd9c099399a60..be01dcda29d92056b7b8457c7c2ce285ab97e8b7 100755
--- a/Examples/varia/PolarizedNonperfectAnalyzerPolarizer.py
+++ b/Examples/varia/PolarizedNonperfectAnalyzerPolarizer.py
@@ -71,9 +71,9 @@ def get_simulation(sample, scan_size=1500):
 
 
 def run_simulation(*,
-                   polarization=ba.kvector_t(0, 1, 0),
+                   polarizer_dir=ba.kvector_t(0, 1, 0),
                    polarizer_efficiency=1,
-                   analyzer=ba.kvector_t(0, 1, 0),
+                   analyzer_dir=ba.kvector_t(0, 1, 0),
                    analyzer_efficiency=1):
     """
     Runs simulation and returns its result.
@@ -81,8 +81,8 @@ def run_simulation(*,
     sample = get_sample()
     simulation = get_simulation(sample)
 
-    simulation.beam().setPolarization(polarization*polarizer_efficiency)
-    simulation.detector().setAnalyzer(analyzer,
+    simulation.beam().setPolarization(polarizer_dir*polarizer_efficiency)
+    simulation.detector().setAnalyzer(analyzer_dir,
                                                 analyzer_efficiency, 0.5)
 
     simulation.setBackground(ba.ConstantBackground(1e-7))
@@ -116,21 +116,21 @@ if __name__ == '__main__':
     polarizer_efficiency = 0.986
     analyzer_efficiency = 0.970
 
-    results_pp = run_simulation(polarization=ba.kvector_t(0, 1, 0),
-                                analyzer=ba.kvector_t(0, 1, 0),
+    results_pp = run_simulation(polarizer_dir=ba.kvector_t(0, 1, 0),
+                                analyzer_dir=ba.kvector_t(0, 1, 0),
                                 polarizer_efficiency=polarizer_efficiency,
                                 analyzer_efficiency=analyzer_efficiency)
-    results_mm = run_simulation(polarization=ba.kvector_t(0, -1, 0),
-                                analyzer=ba.kvector_t(0, -1, 0),
+    results_mm = run_simulation(polarizer_dir=ba.kvector_t(0, -1, 0),
+                                analyzer_dir=ba.kvector_t(0, -1, 0),
                                 polarizer_efficiency=polarizer_efficiency,
                                 analyzer_efficiency=analyzer_efficiency)
 
-    results_pm = run_simulation(polarization=ba.kvector_t(0, 1, 0),
-                                analyzer=ba.kvector_t(0, -1, 0),
+    results_pm = run_simulation(polarizer_dir=ba.kvector_t(0, 1, 0),
+                                analyzer_dir=ba.kvector_t(0, -1, 0),
                                 polarizer_efficiency=polarizer_efficiency,
                                 analyzer_efficiency=analyzer_efficiency)
-    results_mp = run_simulation(polarization=ba.kvector_t(0, -1, 0),
-                                analyzer=ba.kvector_t(0, 1, 0),
+    results_mp = run_simulation(polarizer_dir=ba.kvector_t(0, -1, 0),
+                                analyzer_dir=ba.kvector_t(0, 1, 0),
                                 polarizer_efficiency=polarizer_efficiency,
                                 analyzer_efficiency=analyzer_efficiency)
 
diff --git a/Examples/varia/PolarizedSpinAsymmetry.py b/Examples/varia/PolarizedSpinAsymmetry.py
index 307e2490099608f9c10c06641815d834a6c5effd..1ad373cc870c8090d3e2757a064d92eb1c6a34fc 100755
--- a/Examples/varia/PolarizedSpinAsymmetry.py
+++ b/Examples/varia/PolarizedSpinAsymmetry.py
@@ -63,7 +63,7 @@ def get_sample(params):
     return multi_layer
 
 
-def get_simulation(q_axis, parameters, polarization, analyzer):
+def get_simulation(q_axis, parameters, polarizer_dir, analyzer_dir):
     """
     Returns a simulation object.
     Polarization, analyzer and resolution are set
@@ -80,14 +80,14 @@ def get_simulation(q_axis, parameters, polarization, analyzer):
     distr = ba.RangedDistributionGaussian(n_samples, n_sig)
     scan.setAbsoluteQResolution(distr, parameters["q_res"])
 
-    simulation.beam().setPolarization(polarization)
-    simulation.detector().setAnalyzer(analyzer, 1, 0.5)
+    simulation.beam().setPolarization(polarizer_dir)
+    simulation.detector().setAnalyzer(analyzer_dir, 1, 0.5)
 
     simulation.setScan(scan)
     return simulation
 
 
-def run_simulation(q_axis, fitParams, *, polarization, analyzer):
+def run_simulation(q_axis, fitParams, *, polarizer_dir, analyzer_dir):
     """
     Run a simulation on the given q-axis, where the sample is
     constructed with the given parameters.
@@ -96,7 +96,7 @@ def run_simulation(q_axis, fitParams, *, polarization, analyzer):
     parameters = dict(fitParams, **fixedParams)
 
     sample = get_sample(parameters)
-    simulation = get_simulation(q_axis, parameters, polarization, analyzer)
+    simulation = get_simulation(q_axis, parameters, polarizer_dir, analyzer_dir)
 
     simulation.setSample(sample)
     simulation.runSimulation()
@@ -224,14 +224,14 @@ if __name__ == '__main__':
     def run_Simulation_pp(qzs, params):
         return run_simulation(qzs,
                               params,
-                              polarization=ba.kvector_t(0, 1, 0),
-                              analyzer=ba.kvector_t(0, 1, 0))
+                              polarizer_dir=ba.kvector_t(0, 1, 0),
+                              analyzer_dir=ba.kvector_t(0, 1, 0))
 
     def run_Simulation_mm(qzs, params):
         return run_simulation(qzs,
                               params,
-                              polarization=ba.kvector_t(0, -1, 0),
-                              analyzer=ba.kvector_t(0, -1, 0))
+                              polarizer_dir=ba.kvector_t(0, -1, 0),
+                              analyzer_dir=ba.kvector_t(0, -1, 0))
 
     qzs = numpy.linspace(qmin, qmax, scan_size)
     q_pp, r_pp = qr(run_Simulation_pp(qzs, fixedParams))
diff --git a/Examples/varia/PolarizedSpinFlip.py b/Examples/varia/PolarizedSpinFlip.py
index 5ed7f6f75a9ab3b382d07936dd41507297b81332..32e58dfbdb2a596fee414ecccd2aa3059e1667df 100755
--- a/Examples/varia/PolarizedSpinFlip.py
+++ b/Examples/varia/PolarizedSpinFlip.py
@@ -47,8 +47,8 @@ def get_simulation(sample, scan_size=500):
     return simulation
 
 
-def run_simulation(polarization=ba.kvector_t(0, 1, 0),
-                   analyzer=ba.kvector_t(0, 1, 0)):
+def run_simulation(polarizer_dir=ba.kvector_t(0, 1, 0),
+                   analyzer_dir=ba.kvector_t(0, 1, 0)):
     """
     Runs simulation and returns its result.
     """
@@ -56,8 +56,8 @@ def run_simulation(polarization=ba.kvector_t(0, 1, 0),
     simulation = get_simulation(sample)
 
     # adding polarization and analyzer operator
-    simulation.beam().setPolarization(polarization)
-    simulation.detector().setAnalyzer(analyzer, 1, 0.5)
+    simulation.beam().setPolarization(polarizer_dir)
+    simulation.detector().setAnalyzer(analyzer_dir, 1, 0.5)
 
     simulation.runSimulation()
     return simulation.result()
diff --git a/Tests/Functional/PyCore/polmagcylinders2.py b/Tests/Functional/PyCore/polmagcylinders2.py
index bf34dc910114ea00fab19f488fcc4dea088c7cdc..eded466d8b4d415b4b42f70b63becc2550a4b427 100644
--- a/Tests/Functional/PyCore/polmagcylinders2.py
+++ b/Tests/Functional/PyCore/polmagcylinders2.py
@@ -12,7 +12,7 @@ from bornagain import *
 # ----------------------------------
 # describe sample and run simulation
 # ----------------------------------
-def getSimulationIntensity(rho_beam, efficiency):
+def getSimulationIntensity(polarizer_dir, efficiency):
     print("- simulate", flush=True)
     # defining materials
     mAmbience = HomogeneousMaterial("Vacuum", 0, 0)
@@ -45,7 +45,7 @@ def getSimulationIntensity(rho_beam, efficiency):
     zplus = kvector_t(0, 0, 1)
     simulation.detector().setAnalyzer(zplus, efficiency, 0.5)
     simulation.setBeamParameters(1*angstrom, 0.2*deg, 0)
-    simulation.beam().setPolarization(rho_beam)
+    simulation.beam().setPolarization(polarizer_dir)
     simulation.setSample(multi_layer)
     simulation.beam().setIntensity(1e9)
     simulation.runSimulation()
diff --git a/Tests/Unit/Core/SpecularSimulationTest.cpp b/Tests/Unit/Core/SpecularSimulationTest.cpp
index 2a1e4ad767e9a27ebc153bf2a90898a4c1a671e9..e0678e83d09d51a72ce649425d4554d161959e6e 100644
--- a/Tests/Unit/Core/SpecularSimulationTest.cpp
+++ b/Tests/Unit/Core/SpecularSimulationTest.cpp
@@ -98,10 +98,10 @@ TEST_F(SpecularSimulationTest, SetAngularScan)
 
     SpecularSimulation sim4;
     AlphaScan scan4(1.0, 10, .0 * Units::deg, 2.0 * Units::deg);
-    const auto polarization = kvector_t({0., 0., 0.876});
-    const auto analyzer = kvector_t({0., 0., 1.});
-    sim4.beam().setPolarization(polarization);
-    sim4.detector().setAnalyzer(analyzer, 0.33, 0.22);
+    const auto polarizer_dir = kvector_t({0., 0., 0.876});
+    const auto analyzer_dir = kvector_t({0., 0., 1.});
+    sim4.beam().setPolarization(polarizer_dir);
+    sim4.detector().setAnalyzer(analyzer_dir, 0.33, 0.22);
     sim4.setScan(scan4);
     EXPECT_THROW(sim4.setScan(scan4), std::runtime_error);
 
@@ -113,8 +113,8 @@ TEST_F(SpecularSimulationTest, SetAngularScan)
     EXPECT_EQ(0.0, sim4.beam().direction().alpha());
     EXPECT_EQ(0.0, sim4.beam().direction().phi());
 
-    EXPECT_EQ(sim4.beam().getBlochVector(), polarization);
-    EXPECT_EQ(sim4.detector().analyzer().polDirection(), analyzer);
+    EXPECT_EQ(sim4.beam().getBlochVector(), polarizer_dir);
+    EXPECT_EQ(sim4.detector().analyzer().polDirection(), analyzer_dir);
     EXPECT_EQ(sim4.detector().analyzer().polEfficiency(), 0.33);
     EXPECT_EQ(sim4.detector().analyzer().totalTransmission(), 0.22);
 }
@@ -152,10 +152,10 @@ TEST_F(SpecularSimulationTest, SetQScan)
 
     SpecularSimulation sim3;
     QzScan scan3(10, 1.0, 10.0);
-    const auto polarization = kvector_t({0., 0., 0.876});
-    const auto analyzer = kvector_t({0., 0., 1.});
-    sim3.beam().setPolarization(polarization);
-    sim3.detector().setAnalyzer(analyzer, 0.33, 0.22);
+    const auto polarizer_dir = kvector_t({0., 0., 0.876});
+    const auto analyzer_dir = kvector_t({0., 0., 1.});
+    sim3.beam().setPolarization(polarizer_dir);
+    sim3.detector().setAnalyzer(analyzer_dir, 0.33, 0.22);
     sim3.setScan(scan3);
 
     EXPECT_EQ(1.0, sim3.coordinateAxis()->lowerBound());
@@ -166,8 +166,8 @@ TEST_F(SpecularSimulationTest, SetQScan)
     EXPECT_EQ(0.0, sim3.beam().direction().alpha());
     EXPECT_EQ(0.0, sim3.beam().direction().phi());
 
-    EXPECT_EQ(sim3.beam().getBlochVector(), polarization);
-    EXPECT_EQ(sim3.detector().analyzer().polDirection(), analyzer);
+    EXPECT_EQ(sim3.beam().getBlochVector(), polarizer_dir);
+    EXPECT_EQ(sim3.detector().analyzer().polDirection(), analyzer_dir);
     EXPECT_EQ(sim3.detector().analyzer().polEfficiency(), 0.33);
     EXPECT_EQ(sim3.detector().analyzer().totalTransmission(), 0.22);
 }
diff --git a/Tests/Unit/Device/BeamTest.cpp b/Tests/Unit/Device/BeamTest.cpp
index 420e70e6fdc0800089d69756f710a335089046e5..39158e26172dcb98908f9df79b4d9dfa7cd8f8e7 100644
--- a/Tests/Unit/Device/BeamTest.cpp
+++ b/Tests/Unit/Device/BeamTest.cpp
@@ -12,8 +12,8 @@ class BeamTest : public ::testing::Test {
 TEST_F(BeamTest, BeamPolarization)
 {
     Beam beam = Beam::horizontalBeam();
-    kvector_t polarization(0.1, -0.2, 0.4);
-    beam.setPolarization(polarization);
+    kvector_t polarizer_dir(0.1, -0.2, 0.4);
+    beam.setPolarization(polarizer_dir);
 
     kvector_t bloch_vector = beam.getBlochVector();
     EXPECT_NEAR(0.1, bloch_vector.x(), 1e-8);