diff --git a/Examples/anaklasis/ba_anaklasis.py b/Examples/anaklasis/ba_anaklasis.py
index 0257cb237100799bfdd5e2d364975dfbe0dcf1c7..0822d7671f718465571ed987f821a352706e9b5b 100755
--- a/Examples/anaklasis/ba_anaklasis.py
+++ b/Examples/anaklasis/ba_anaklasis.py
@@ -449,7 +449,6 @@ def simulate(sample,
     if qmax:
         qzs = np.linspace(0., qmax[0]/angstrom, scan_size)
 
-    simulation = ba.SpecularSimulation()
 
     n_sig = 4.0
     n_samples = 25
@@ -462,11 +461,9 @@ def simulate(sample,
     else:
         scan.setRelativeQResolution(distr, 0.5*resolution)
 
-    simulation.setScan(scan)
-
+    simulation = ba.SpecularSimulation(scan, sample)
     simulation.setBackground(ba.ConstantBackground(background/scale))
 
-    simulation.setSample(sample)
     result = simulation.simulate()
 
     return np.array(result.axis(
diff --git a/Examples/bayesian/likelihood_sampling.py b/Examples/bayesian/likelihood_sampling.py
index 760a909cb8c06a87fba3c9189e083584bf43d2a4..be455ef29aa87bff79ce372758d5aa04fc05c4d1 100755
--- a/Examples/bayesian/likelihood_sampling.py
+++ b/Examples/bayesian/likelihood_sampling.py
@@ -76,15 +76,13 @@ def get_real_data():
 
 
 # Define the simulation
-def get_simulation(alpha):
+def get_simulation(sample, alpha):
     """
     Defines and returns a specular simulation.
     """
     wavelength = 0.154  #nm
     scan = ba.AlphaScan(wavelength, alpha)
-    simulation = ba.SpecularSimulation()
-    simulation.setScan(scan)
-    return simulation
+    return ba.SpecularSimulation(scan, sample)
 
 
 # Run the simulation
@@ -97,8 +95,8 @@ def run_simulation(alpha, ni_thickness, ti_thickness):
     :return: simulated reflected intensity
     """
     sample = get_sample(ni_thickness, ti_thickness)
-    simulation = get_simulation(alpha)
-    simulation.setSample(sample)
+    simulation = get_simulation(sample, alpha)
+
     result = simulation.simulate()
     return result.array()
 
diff --git a/Examples/fit/specular/FitSpecularBasics.py b/Examples/fit/specular/FitSpecularBasics.py
index b086df09c3c40bbd25c47004057ab09134050b71..022344a0440cc353cbf359b836b851f61093f711 100755
--- a/Examples/fit/specular/FitSpecularBasics.py
+++ b/Examples/fit/specular/FitSpecularBasics.py
@@ -92,12 +92,10 @@ def get_simulation(params):
     Create and return specular simulation with its instrument defined
     """
     wavelength = 1.54*ba.angstrom  # beam wavelength
-
-    simulation = ba.SpecularSimulation()
     scan = ba.AlphaScan(wavelength, get_real_data_axis())
-    simulation.setScan(scan)
-    simulation.setSample(get_sample(params))
-    return simulation
+    sample = get_sample(params)
+
+    return ba.SpecularSimulation(scan, sample)
 
 
 def run_fitting():
diff --git a/Examples/fit/specular/Honeycomb_fit.py b/Examples/fit/specular/Honeycomb_fit.py
index 3f01274449a186a73a51e7949eb33847bcfd36fc..f8755973fab99b6ea405dc403347b4156097f5bb 100755
--- a/Examples/fit/specular/Honeycomb_fit.py
+++ b/Examples/fit/specular/Honeycomb_fit.py
@@ -88,9 +88,6 @@ def get_simulation(q_axis, fitParams, sign, ms150=False):
     scan = ba.QzScan(q_axis)
     scan.setAbsoluteQResolution(q_distr, dq)
 
-    simulation = ba.SpecularSimulation()
-    simulation.setScan(scan)
-    simulation.beam().setIntensity(parameters["intensity"])
 
     if ms150:
         sample = get_sample(parameters=parameters,
@@ -99,7 +96,8 @@ def get_simulation(q_axis, fitParams, sign, ms150=False):
     else:
         sample = get_sample(parameters=parameters, sign=sign, ms150=1)
 
-    simulation.setSample(sample)
+    simulation = ba.SpecularSimulation(scan, sample)
+    simulation.beam().setIntensity(parameters["intensity"])
     simulation.setBackground(ba.ConstantBackground(5e-7))
 
     return simulation
diff --git a/Examples/fit/specular/Pt_layer_fit.py b/Examples/fit/specular/Pt_layer_fit.py
index cb0b31afb9f40e9446286fe353171ef440336daa..55ee70c557e33b7b0b6a54ded6c0d48f14fff675 100755
--- a/Examples/fit/specular/Pt_layer_fit.py
+++ b/Examples/fit/specular/Pt_layer_fit.py
@@ -56,6 +56,8 @@ def get_sample(params):
 
 def get_simulation(q_axis, parameters):
 
+    sample = get_sample(parameters)
+
     scan = ba.QzScan(q_axis)
     scan.setOffset(parameters["q_offset"])
 
@@ -65,21 +67,18 @@ def get_simulation(q_axis, parameters):
     distr = ba.RangedDistributionGaussian(n_samples, n_sig)
     scan.setAbsoluteQResolution(distr, parameters["q_res/q"])
 
-    simulation = ba.SpecularSimulation()
+    simulation = ba.SpecularSimulation(scan, sample)
+
     simulation.beam().setIntensity(parameters["intensity"])
 
-    simulation.setScan(scan)
     return simulation
 
 
 def run_simulation(q_axis, fitParams):
     parameters = dict(fitParams, **fixedParams)
 
-    sample = get_sample(parameters)
     simulation = get_simulation(q_axis, parameters)
 
-    simulation.setSample(sample)
-
     return simulation.simulate()
 
 
diff --git a/Examples/fit/specular/RealLifeReflectometryFitting.py b/Examples/fit/specular/RealLifeReflectometryFitting.py
index 787d8123cd6581fdeb8bc294a7cf959cf8c8ef80..348d6aaa34a266d62fe675246610b38da5796d53 100755
--- a/Examples/fit/specular/RealLifeReflectometryFitting.py
+++ b/Examples/fit/specular/RealLifeReflectometryFitting.py
@@ -88,7 +88,7 @@ def get_weights(start, end):
     return expdata[start:end, 2]
 
 
-def create_simulation(arg_dict, bin_start, bin_end):
+def create_simulation(sample, arg_dict, bin_start, bin_end):
     """
     Creates and returns specular simulation
     """
@@ -100,8 +100,7 @@ def create_simulation(arg_dict, bin_start, bin_end):
     scan.setAbsoluteAngularResolution(alpha_distr, arg_dict["divergence"])
     scan.setFootprintFactor(footprint)
 
-    simulation = ba.SpecularSimulation()
-    simulation.setScan(scan)
+    simulation = ba.SpecularSimulation(scan, sample)
     simulation.beam().setIntensity(arg_dict["intensity"])
     return simulation
 
@@ -139,8 +138,8 @@ def run_simulation(arg_dict, bin_start=0, bin_end=-1):
     Runs simulation and returns its result
     """
 
-    simulation = create_simulation(arg_dict, bin_start, bin_end)
-    simulation.setSample(buildSample(arg_dict))
+    sample = buildSample(arg_dict)
+    simulation = create_simulation(sample, arg_dict, bin_start, bin_end)
 
     return simulation.simulate()
 
diff --git a/Examples/specular/AlternatingLayers1.py b/Examples/specular/AlternatingLayers1.py
index 536e3fc78142ea0e45b98df6ab1413a878c0661c..b25da0888d52673db1689fe9add428a7ebc01171 100755
--- a/Examples/specular/AlternatingLayers1.py
+++ b/Examples/specular/AlternatingLayers1.py
@@ -42,11 +42,8 @@ def get_simulation(sample):
     A standard specular simulation setup.
     """
     n = bp.simargs['n']
-    simulation = ba.SpecularSimulation()
     scan = ba.AlphaScan(1.54*angstrom, n, 0, 2*deg)
-    simulation.setScan(scan)
-    simulation.setSample(sample)
-    return simulation
+    return ba.SpecularSimulation(scan, sample)
 
 
 if __name__ == '__main__':
diff --git a/Examples/specular/AlternatingLayers2.py b/Examples/specular/AlternatingLayers2.py
index 47023a901a3004406127fc456c0adce294dd54cc..acac1dfc64d237698f274876cfa04c6efefa9b0c 100755
--- a/Examples/specular/AlternatingLayers2.py
+++ b/Examples/specular/AlternatingLayers2.py
@@ -18,11 +18,8 @@ def get_simulation(sample):
     A standard specular simulation setup.
     """
     n = bp.simargs['n']
-    simulation = ba.SpecularSimulation()
     scan = ba.AlphaScan(1.54*angstrom, n, 0, 2*deg)
-    simulation.setScan(scan)
-    simulation.setSample(sample)
-    return simulation
+    return ba.SpecularSimulation(scan, sample)
 
 
 if __name__ == '__main__':
diff --git a/Examples/specular/BasicPolarizedReflectometry.py b/Examples/specular/BasicPolarizedReflectometry.py
index a20d9685c81a045038d78f7e240f5c49edc7c3d1..af7248250a9322cd75df652b154a0ca2c4a55155 100755
--- a/Examples/specular/BasicPolarizedReflectometry.py
+++ b/Examples/specular/BasicPolarizedReflectometry.py
@@ -36,18 +36,10 @@ def get_sample():
     return sample
 
 
-def get_simulation(sample, scan_size=500):
-    """
-    Defines and returns a specular simulation.
-    """
-
-
-def simulate(polarizer_dir, analyzer_dir, title):
+def simulate(sample, polarizer_dir, analyzer_dir, title):
     n = bp.simargs['n']
-    simulation = ba.SpecularSimulation()
     scan = ba.AlphaScan(1.54*angstrom, n, 0, 5*deg)
-    simulation.setScan(scan)
-    simulation.setSample(get_sample())
+    simulation = ba.SpecularSimulation(scan, sample)
 
     simulation.setPolarization(polarizer_dir)
     simulation.setAnalyzer(analyzer_dir, 1, 0.5)
@@ -57,16 +49,13 @@ def simulate(polarizer_dir, analyzer_dir, title):
     return result
 
 
-def run_simulations():
-    return ret
-
-
 if __name__ == '__main__':
     bp.parse_args(sim_n=500)
 
+    sample = get_sample()
     results = [
-        simulate(ba.R3(0, +1, 0), ba.R3(0, +1, 0), "$++$"),
-        simulate(ba.R3(0, -1, 0), ba.R3(0, -1, 0), "$--$"),
+        simulate(sample, ba.R3(0, +1, 0), ba.R3(0, +1, 0), "$++$"),
+        simulate(sample, ba.R3(0, -1, 0), ba.R3(0, -1, 0), "$--$"),
     ]
 
     bp.plot_multicurve_specular(results)
diff --git a/Examples/specular/BeamAngularDivergence.py b/Examples/specular/BeamAngularDivergence.py
index b8e3f6a38a3cb9dd312477fd901caa22709afaff..557fcc62be80823b0ba95ee22cfb5637755beb9f 100755
--- a/Examples/specular/BeamAngularDivergence.py
+++ b/Examples/specular/BeamAngularDivergence.py
@@ -51,11 +51,7 @@ def get_simulation(sample, **kwargs):
     scan.setFootprintFactor(footprint)
     scan.setAbsoluteAngularResolution(alpha_distr, d_ang)
 
-    simulation = ba.SpecularSimulation()
-    simulation.setScan(scan)
-    simulation.setSample(sample)
-
-    return simulation
+    return ba.SpecularSimulation(scan, sample)
 
 
 if __name__ == '__main__':
diff --git a/Examples/specular/BeamFullDivergence.py b/Examples/specular/BeamFullDivergence.py
index bb3a9b3f6ba177bce998d215e4c680d3c4dca34a..cd09f06fd02d2914f16b156ca1eb145971db8857 100755
--- a/Examples/specular/BeamFullDivergence.py
+++ b/Examples/specular/BeamFullDivergence.py
@@ -35,9 +35,7 @@ def get_simulation(sample):
     scan.setAbsoluteAngularResolution(alpha_distr, d_ang)
     scan.setAbsoluteWavelengthResolution(wavelength_distr, d_wl)
 
-    simulation = ba.SpecularSimulation()
-    simulation.setScan(scan)
-    simulation.setSample(sample)
+    return ba.SpecularSimulation(scan, sample)
 
     return simulation
 
diff --git a/Examples/specular/FootprintCorrection.py b/Examples/specular/FootprintCorrection.py
index 793d5cdcfb295fce01d2ca54016177aa74526bf0..c9889c4fda2e1d8fee110b8f20a5128702490672 100755
--- a/Examples/specular/FootprintCorrection.py
+++ b/Examples/specular/FootprintCorrection.py
@@ -12,12 +12,11 @@ sample = std_samples.alternating_layers()
 
 
 def simulate(footprint, title):
-    simulation = ba.SpecularSimulation()
     n = bp.simargs['n']
     scan = ba.AlphaScan(1.54*angstrom, n, 0, 0.6*deg)
     scan.setFootprintFactor(footprint)
-    simulation.setScan(scan)
-    simulation.setSample(sample)
+    simulation = ba.SpecularSimulation(scan, sample)
+
     result = simulation.simulate()
     result.setTitle(title)
     return result
diff --git a/Examples/specular/PolarizedNoAnalyzer.py b/Examples/specular/PolarizedNoAnalyzer.py
index 78a46aa29f05d4206e4e664af483848f3d97cd55..f2642fd5441375095a84ac11cf5c6a5cc74280e6 100755
--- a/Examples/specular/PolarizedNoAnalyzer.py
+++ b/Examples/specular/PolarizedNoAnalyzer.py
@@ -40,12 +40,9 @@ def get_simulation(sample):
     """
     Defines and returns a specular simulation.
     """
-    simulation = ba.SpecularSimulation()
     n = bp.simargs['n']
     scan = ba.AlphaScan(1.54*angstrom, n, 0, 5*deg)
-    simulation.setScan(scan)
-    simulation.setSample(sample)
-    return simulation
+    return ba.SpecularSimulation(scan, sample)
 
 
 def run_simulation(polarizer_dir=R3(0, 1, 0), analyzer_dir=None):
diff --git a/Examples/specular/PolarizedNonperfectAnalyzerPolarizer.py b/Examples/specular/PolarizedNonperfectAnalyzerPolarizer.py
index 03e0f386754917f583c68232b577f5e00824dd88..b7e8f34b56255fd361cb5365eefde291cae2a719 100755
--- a/Examples/specular/PolarizedNonperfectAnalyzerPolarizer.py
+++ b/Examples/specular/PolarizedNonperfectAnalyzerPolarizer.py
@@ -54,7 +54,6 @@ def get_simulation(sample, scan_size=1500):
     """
     Defines and returns a specular simulation.
     """
-    simulation = ba.SpecularSimulation()
     qzs = numpy.linspace(0.1, 1.5, scan_size)
 
     n_sig = 4.0
@@ -64,9 +63,7 @@ def get_simulation(sample, scan_size=1500):
     scan = ba.QzScan(qzs)
     scan.setAbsoluteQResolution(distr, 0.008)
 
-    simulation.setScan(scan)
-    simulation.setSample(sample)
-    return simulation
+    return ba.SpecularSimulation(scan, sample)
 
 
 def run_simulation(*,
diff --git a/Examples/specular/PolarizedSpinAsymmetry.py b/Examples/specular/PolarizedSpinAsymmetry.py
index b1c0fdaeac7952493d437cf4a79d1119a677faaa..a6e2bdeb1b4c100251c6db6dfa8c668c23845a05 100755
--- a/Examples/specular/PolarizedSpinAsymmetry.py
+++ b/Examples/specular/PolarizedSpinAsymmetry.py
@@ -63,13 +63,12 @@ def get_sample(params):
     return multi_layer
 
 
-def get_simulation(q_axis, parameters, polarizer_dir, analyzer_dir):
+def get_simulation(sample, q_axis, parameters, polarizer_dir, analyzer_dir):
     """
     Returns a simulation object.
     Polarization, analyzer and resolution are set
     from given parameters
     """
-    simulation = ba.SpecularSimulation()
     q_axis = q_axis + parameters["q_offset"]
     scan = ba.QzScan(q_axis)
 
@@ -80,10 +79,10 @@ def get_simulation(q_axis, parameters, polarizer_dir, analyzer_dir):
     distr = ba.RangedDistributionGaussian(n_samples, n_sig)
     scan.setAbsoluteQResolution(distr, parameters["q_res"])
 
+    simulation = ba.SpecularSimulation(scan, sample)
     simulation.setPolarization(polarizer_dir)
     simulation.setAnalyzer(analyzer_dir, 1, 0.5)
 
-    simulation.setScan(scan)
     return simulation
 
 
@@ -96,10 +95,9 @@ def run_simulation(q_axis, fitParams, *, polarizer_dir, analyzer_dir):
     parameters = dict(fitParams, **fixedParams)
 
     sample = get_sample(parameters)
-    simulation = get_simulation(q_axis, parameters, polarizer_dir,
+    simulation = get_simulation(sample, q_axis, parameters, polarizer_dir,
                                 analyzer_dir)
 
-    simulation.setSample(sample)
     return simulation.simulate()
 
 
diff --git a/Examples/specular/PolarizedSpinFlip.py b/Examples/specular/PolarizedSpinFlip.py
index 05dee45de5d95306ca0c881dd5c2c68909faccd2..96cf227411f319e15d0e4fd8a3b6f5ebb490bc48 100755
--- a/Examples/specular/PolarizedSpinFlip.py
+++ b/Examples/specular/PolarizedSpinFlip.py
@@ -36,25 +36,14 @@ def get_sample():
     return sample
 
 
-def get_simulation(sample, scan_size=500):
-    """
-    Defines and returns a specular simulation.
-    """
-    simulation = ba.SpecularSimulation()
-    scan = ba.AlphaScan(1.54*angstrom, scan_size, 0, 5*deg)
-    simulation.setScan(scan)
-    simulation.setSample(sample)
-    return simulation
-
-
-def run_simulation(polarizer_dir=ba.R3(0, 1, 0),
+def run_simulation(sample, polarizer_dir=ba.R3(0, 1, 0),
                    analyzer_dir=ba.R3(0, 1, 0)):
     """
     Runs simulation and returns its result.
     """
-    sample = get_sample()
-    simulation = get_simulation(sample)
+    scan = ba.AlphaScan(1.54*angstrom, 500, 0, 5*deg)
 
+    simulation = ba.SpecularSimulation(scan, sample)
     simulation.setPolarization(polarizer_dir)
     simulation.setAnalyzer(analyzer_dir, 1, 0.5)
 
@@ -81,11 +70,12 @@ def plot(data, labels):
 if __name__ == '__main__':
     bp.parse_args()
 
-    results_pp = run_simulation(ba.R3(0, 1, 0), ba.R3(0, 1, 0))
-    results_mm = run_simulation(ba.R3(0, -1, 0), ba.R3(0, -1, 0))
+    sample = get_sample()
 
-    results_pm = run_simulation(ba.R3(0, 1, 0), ba.R3(0, -1, 0))
-    results_mp = run_simulation(ba.R3(0, -1, 0), ba.R3(0, 1, 0))
+    results_pp = run_simulation(sample, ba.R3(0, 1, 0), ba.R3(0, 1, 0))
+    results_mm = run_simulation(sample, ba.R3(0, -1, 0), ba.R3(0, -1, 0))
+    results_pm = run_simulation(sample, ba.R3(0, 1, 0), ba.R3(0, -1, 0))
+    results_mp = run_simulation(sample, ba.R3(0, -1, 0), ba.R3(0, 1, 0))
 
     plot([results_pp, results_mm, results_pm, results_mp],
          ["$++$", "$--$", "$+-$", "$-+$"])
diff --git a/Examples/specular/RoughnessModel.py b/Examples/specular/RoughnessModel.py
index 1a9a208a6a4615198aaedd87cd5762818c5befd2..0291e12efe9993e0c4dbb4d9e36e3f21adcecdf4 100755
--- a/Examples/specular/RoughnessModel.py
+++ b/Examples/specular/RoughnessModel.py
@@ -47,12 +47,10 @@ def get_simulation(sample):
     """
     Defines and returns a specular simulation.
     """
-    simulation = ba.SpecularSimulation()
     n = bp.simargs['n']
     scan = ba.AlphaScan(1.54*angstrom, n, 0, 2*deg)
-    simulation.setScan(scan)
-    simulation.setSample(sample)
-    return simulation
+
+    return ba.SpecularSimulation(scan, sample)
 
 
 def simulate(roughness_model, title):
diff --git a/Examples/specular/SpecularSimulationWithRoughness.py b/Examples/specular/SpecularSimulationWithRoughness.py
index 57cb74e83f6ee6147a9e4f815bedb1b26a16f670..304598c0e3d552e78f7e1090a635d2f905bbcbfc 100755
--- a/Examples/specular/SpecularSimulationWithRoughness.py
+++ b/Examples/specular/SpecularSimulationWithRoughness.py
@@ -35,12 +35,9 @@ def get_sample():
 
 
 def get_simulation(sample):
-    simulation = ba.SpecularSimulation()
     n = bp.simargs['n']
     scan = ba.AlphaScan(1.54*angstrom, n, 0, 2*deg)
-    simulation.setScan(scan)
-    simulation.setSample(sample)
-    return simulation
+    return ba.SpecularSimulation(scan, sample)
 
 
 if __name__ == '__main__':
diff --git a/Examples/specular/TOFRWithResolution.py b/Examples/specular/TOFRWithResolution.py
index 558cadeed1f694ba3c5535a199781dd59d72055d..ddcb36870904102d1cafad9e19016607aca62864 100755
--- a/Examples/specular/TOFRWithResolution.py
+++ b/Examples/specular/TOFRWithResolution.py
@@ -36,11 +36,7 @@ def get_simulation(sample):
     scan = ba.QzScan(qzs)
     scan.setAbsoluteQResolution(distr, dq)
 
-    simulation = ba.SpecularSimulation()
-    simulation.setScan(scan)
-    simulation.setSample(sample)
-
-    return simulation
+    return ba.SpecularSimulation(scan, sample)
 
 
 if __name__ == '__main__':
diff --git a/Examples/specular/TimeOfFlightReflectometry.py b/Examples/specular/TimeOfFlightReflectometry.py
index 4926047765132d38ba4f398775463d8c5518a360..f656a4dac7c9455de3c12397ef97287a64f62caa 100755
--- a/Examples/specular/TimeOfFlightReflectometry.py
+++ b/Examples/specular/TimeOfFlightReflectometry.py
@@ -25,10 +25,7 @@ def get_simulation(sample):
     n = bp.simargs['n']
     qzs = np.linspace(0.01, 1, n)  # qz-values
     scan = ba.QzScan(qzs)
-    simulation = ba.SpecularSimulation()
-    simulation.setScan(scan)
-    simulation.setSample(sample)
-    return simulation
+    return ba.SpecularSimulation(scan, sample)
 
 
 if __name__ == '__main__':
diff --git a/Tests/Unit/Sim/SpecularSimulationTest.cpp b/Tests/Unit/Sim/SpecularSimulationTest.cpp
index 6f6e08228547f6b1dfa362233a91bb15d557bb09..370f643de33cac451f427669e47d6270a64ac34a 100644
--- a/Tests/Unit/Sim/SpecularSimulationTest.cpp
+++ b/Tests/Unit/Sim/SpecularSimulationTest.cpp
@@ -23,7 +23,7 @@ protected:
 
     std::unique_ptr<SpecularSimulation> defaultSimulation();
 
-    MultiLayer multilayer;
+    MultiLayer sample;
 };
 
 SpecularSimulationTest::SpecularSimulationTest()
@@ -36,25 +36,22 @@ SpecularSimulationTest::SpecularSimulationTest()
     Layer layer1(mat1, 10 * Units::nm);
     Layer layer2(mat2);
 
-    multilayer.addLayer(layer0);
-    multilayer.addLayer(layer1);
-    multilayer.addLayer(layer2);
+    sample.addLayer(layer0);
+    sample.addLayer(layer1);
+    sample.addLayer(layer2);
 }
 
 std::unique_ptr<SpecularSimulation> SpecularSimulationTest::defaultSimulation()
 {
-    auto result = std::make_unique<SpecularSimulation>();
     AlphaScan scan(1.0, FixedBinAxis("axis", 10, 0.0 * Units::deg, 2.0 * Units::deg));
-    result->setScan(scan);
-    result->setSample(multilayer);
-    return result;
+    return std::make_unique<SpecularSimulation>(scan, sample);
 }
 
 TEST_F(SpecularSimulationTest, SetAngularScan)
 {
-    SpecularSimulation sim;
     AlphaScan scan(1.0, std::vector<double>{1.0 * Units::deg, 3.0 * Units::deg});
-    sim.setScan(scan);
+    SpecularSimulation sim(scan, sample);
+
     const auto& beam = sim.beam();
 
     EXPECT_EQ(2u, sim.coordinateAxis()->size());
@@ -65,21 +62,8 @@ TEST_F(SpecularSimulationTest, SetAngularScan)
     EXPECT_EQ(0.0, beam.direction().alpha());
     EXPECT_EQ(0.0, beam.direction().phi());
 
-    SpecularSimulation sim2;
     AlphaScan scan2(1.0, 10, 1.0 * Units::deg, 10.0 * Units::deg);
-    sim2.setScan(scan2);
-    EXPECT_EQ(10u, sim2.coordinateAxis()->size());
-    EXPECT_EQ(1.0 * Units::deg, sim2.coordinateAxis()->min());
-    EXPECT_EQ(10.0 * Units::deg, sim2.coordinateAxis()->max());
-    EXPECT_EQ(1.0, beam.wavelength());
-    EXPECT_EQ(0.0, beam.direction().alpha());
-    EXPECT_EQ(0.0, beam.direction().phi());
-
-    SpecularSimulation sim3;
-    AlphaScan scan3(1.0, 10, -1.0 * Units::deg, 2.0 * Units::deg);
-    EXPECT_FAILED_ASSERT(sim3.setScan(scan3));
-
-    EXPECT_FAILED_ASSERT(sim2.setScan(scan2));
+    SpecularSimulation sim2(scan2, sample);
     EXPECT_EQ(10u, sim2.coordinateAxis()->size());
     EXPECT_EQ(1.0 * Units::deg, sim2.coordinateAxis()->min());
     EXPECT_EQ(10.0 * Units::deg, sim2.coordinateAxis()->max());
@@ -87,14 +71,12 @@ TEST_F(SpecularSimulationTest, SetAngularScan)
     EXPECT_EQ(0.0, beam.direction().alpha());
     EXPECT_EQ(0.0, beam.direction().phi());
 
-    SpecularSimulation sim4;
     AlphaScan scan4(1.0, 10, .0 * Units::deg, 2.0 * Units::deg);
+    SpecularSimulation sim4(scan4, sample);
     const auto polarizer_dir = R3({0., 0., 0.876});
     const auto analyzer_dir = R3({0., 0., 1.});
     sim4.setPolarization(polarizer_dir);
     sim4.setAnalyzer(analyzer_dir, 0.33, 0.22);
-    sim4.setScan(scan4);
-    EXPECT_FAILED_ASSERT(sim4.setScan(scan4));
 
     EXPECT_EQ(.0 * Units::deg, sim4.coordinateAxis()->min());
     EXPECT_EQ(2.0 * Units::deg, sim4.coordinateAxis()->max());
@@ -112,10 +94,8 @@ TEST_F(SpecularSimulationTest, SetAngularScan)
 
 TEST_F(SpecularSimulationTest, SetQScan)
 {
-    SpecularSimulation sim;
-
     QzScan scan(std::vector<double>{1.0, 3.0});
-    sim.setScan(scan);
+    SpecularSimulation sim(scan, sample);
 
     const auto& beam = sim.beam();
 
@@ -130,9 +110,8 @@ TEST_F(SpecularSimulationTest, SetQScan)
     sim.beam().setIntensity(2.0);
     EXPECT_EQ(2.0, beam.intensity());
 
-    SpecularSimulation sim2;
     QzScan scan2(10, 1.0, 10.0);
-    sim2.setScan(scan2);
+    SpecularSimulation sim2(scan2, sample);
     EXPECT_EQ(10u, sim2.coordinateAxis()->size());
     EXPECT_EQ(1.0, sim2.coordinateAxis()->min());
     EXPECT_EQ(10.0, sim2.coordinateAxis()->max());
@@ -141,13 +120,11 @@ TEST_F(SpecularSimulationTest, SetQScan)
     EXPECT_EQ(0.0, beam.direction().alpha());
     EXPECT_EQ(0.0, beam.direction().phi());
 
-    SpecularSimulation sim3;
-    QzScan scan3(10, 1.0, 10.0);
+    SpecularSimulation sim3(scan2, sample);
     const auto polarizer_dir = R3({0., 0., 0.876});
     const auto analyzer_dir = R3({0., 0., 1.});
     sim3.setPolarization(polarizer_dir);
     sim3.setAnalyzer(analyzer_dir, 0.33, 0.22);
-    sim3.setScan(scan3);
 
     EXPECT_EQ(1.0, sim3.coordinateAxis()->min());
     EXPECT_EQ(10.0, sim3.coordinateAxis()->max());
diff --git a/Wrap/Python/std_simulations.py b/Wrap/Python/std_simulations.py
index 4d92918388ac42e7d6b895d145cc582184fb6d48..425868ed269dfc2f68e8853a4588e32a8e52bdcb 100644
--- a/Wrap/Python/std_simulations.py
+++ b/Wrap/Python/std_simulations.py
@@ -9,11 +9,8 @@ def specular(sample, scan_size):
     """
     Returns a standard specular simulation.
     """
-    simulation = ba.SpecularSimulation()
     scan = ba.AlphaScan(1.54 * angstrom, scan_size, 0, 2 * deg)
-    simulation.setScan(scan)
-    simulation.setSample(sample)
-    return simulation
+    return ba.SpecularSimulation(scan, sample)
 
 
 def sas(sample, npix):
diff --git a/auto/Wrap/doxygenSim.i b/auto/Wrap/doxygenSim.i
index 0a233d2465e8bc80339e61545b98d64132c5e662..b73c4c1c192ffb6abd0bc2d2a0f1dd36f9dfa03d 100644
--- a/auto/Wrap/doxygenSim.i
+++ b/auto/Wrap/doxygenSim.i
@@ -1894,6 +1894,9 @@ C++ includes: SpecularSimulation.h
 %feature("docstring")  SpecularSimulation::SpecularSimulation "SpecularSimulation::SpecularSimulation()
 ";
 
+%feature("docstring")  SpecularSimulation::SpecularSimulation "SpecularSimulation::SpecularSimulation(const ISpecularScan &scan, const MultiLayer &sample)
+";
+
 %feature("docstring")  SpecularSimulation::~SpecularSimulation "SpecularSimulation::~SpecularSimulation() override
 ";
 
diff --git a/auto/Wrap/libBornAgainSim.py b/auto/Wrap/libBornAgainSim.py
index 8865dfa8f1af2ab775cd7fdb2414b2eea5bba112..1357e91497246acb3288bd2de1c97085c9f5c659 100644
--- a/auto/Wrap/libBornAgainSim.py
+++ b/auto/Wrap/libBornAgainSim.py
@@ -3688,6 +3688,7 @@ class SpecularSimulation(ISimulation):
     def __init__(self, *args):
         r"""
         __init__(SpecularSimulation self) -> SpecularSimulation
+        __init__(SpecularSimulation self, ISpecularScan const & scan, MultiLayer const & sample) -> SpecularSimulation
         __init__(SpecularSimulation self, SpecularSimulation arg2) -> SpecularSimulation
         SpecularSimulation::SpecularSimulation(SpecularSimulation &&)
 
diff --git a/auto/Wrap/libBornAgainSim_wrap.cpp b/auto/Wrap/libBornAgainSim_wrap.cpp
index c0ca40167831ed3931c09218cb68f5398520b807..45a09be043ca6342853e5b17c8786e3ba095e0fb 100644
--- a/auto/Wrap/libBornAgainSim_wrap.cpp
+++ b/auto/Wrap/libBornAgainSim_wrap.cpp
@@ -38620,6 +38620,41 @@ fail:
 }
 
 
+SWIGINTERN PyObject *_wrap_new_SpecularSimulation__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
+  PyObject *resultobj = 0;
+  ISpecularScan *arg1 = 0 ;
+  MultiLayer *arg2 = 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  void *argp2 = 0 ;
+  int res2 = 0 ;
+  SpecularSimulation *result = 0 ;
+  
+  if ((nobjs < 2) || (nobjs > 2)) SWIG_fail;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1, SWIGTYPE_p_ISpecularScan,  0  | 0);
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "new_SpecularSimulation" "', argument " "1"" of type '" "ISpecularScan const &""'"); 
+  }
+  if (!argp1) {
+    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "new_SpecularSimulation" "', argument " "1"" of type '" "ISpecularScan const &""'"); 
+  }
+  arg1 = reinterpret_cast< ISpecularScan * >(argp1);
+  res2 = SWIG_ConvertPtr(swig_obj[1], &argp2, SWIGTYPE_p_MultiLayer,  0  | 0);
+  if (!SWIG_IsOK(res2)) {
+    SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "new_SpecularSimulation" "', argument " "2"" of type '" "MultiLayer const &""'"); 
+  }
+  if (!argp2) {
+    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "new_SpecularSimulation" "', argument " "2"" of type '" "MultiLayer const &""'"); 
+  }
+  arg2 = reinterpret_cast< MultiLayer * >(argp2);
+  result = (SpecularSimulation *)new SpecularSimulation((ISpecularScan const &)*arg1,(MultiLayer const &)*arg2);
+  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_SpecularSimulation, SWIG_POINTER_NEW |  0 );
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
 SWIGINTERN PyObject *_wrap_delete_SpecularSimulation(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
   SpecularSimulation *arg1 = (SpecularSimulation *) 0 ;
@@ -38642,7 +38677,7 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_new_SpecularSimulation__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
+SWIGINTERN PyObject *_wrap_new_SpecularSimulation__SWIG_2(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
   SpecularSimulation *arg1 = 0 ;
   void *argp1 = 0 ;
@@ -38668,11 +38703,11 @@ fail:
 
 SWIGINTERN PyObject *_wrap_new_SpecularSimulation(PyObject *self, PyObject *args) {
   Py_ssize_t argc;
-  PyObject *argv[2] = {
+  PyObject *argv[3] = {
     0
   };
   
-  if (!(argc = SWIG_Python_UnpackTuple(args, "new_SpecularSimulation", 0, 1, argv))) SWIG_fail;
+  if (!(argc = SWIG_Python_UnpackTuple(args, "new_SpecularSimulation", 0, 2, argv))) SWIG_fail;
   --argc;
   if (argc == 0) {
     return _wrap_new_SpecularSimulation__SWIG_0(self, argc, argv);
@@ -38683,7 +38718,19 @@ SWIGINTERN PyObject *_wrap_new_SpecularSimulation(PyObject *self, PyObject *args
     int res = SWIG_ConvertPtr(argv[0], &vptr, SWIGTYPE_p_SpecularSimulation, SWIG_POINTER_NO_NULL);
     _v = SWIG_CheckState(res);
     if (_v) {
-      return _wrap_new_SpecularSimulation__SWIG_1(self, argc, argv);
+      return _wrap_new_SpecularSimulation__SWIG_2(self, argc, argv);
+    }
+  }
+  if (argc == 2) {
+    int _v;
+    int res = SWIG_ConvertPtr(argv[0], 0, SWIGTYPE_p_ISpecularScan, SWIG_POINTER_NO_NULL | 0);
+    _v = SWIG_CheckState(res);
+    if (_v) {
+      int res = SWIG_ConvertPtr(argv[1], 0, SWIGTYPE_p_MultiLayer, SWIG_POINTER_NO_NULL | 0);
+      _v = SWIG_CheckState(res);
+      if (_v) {
+        return _wrap_new_SpecularSimulation__SWIG_1(self, argc, argv);
+      }
     }
   }
   
@@ -38691,6 +38738,7 @@ fail:
   SWIG_Python_RaiseOrModifyTypeError("Wrong number or type of arguments for overloaded function 'new_SpecularSimulation'.\n"
     "  Possible C/C++ prototypes are:\n"
     "    SpecularSimulation::SpecularSimulation()\n"
+    "    SpecularSimulation::SpecularSimulation(ISpecularScan const &,MultiLayer const &)\n"
     "    SpecularSimulation::SpecularSimulation(SpecularSimulation &&)\n");
   return 0;
 }
@@ -42102,6 +42150,7 @@ static PyMethodDef SwigMethods[] = {
 		""},
 	 { "new_SpecularSimulation", _wrap_new_SpecularSimulation, METH_VARARGS, "\n"
 		"SpecularSimulation()\n"
+		"SpecularSimulation(ISpecularScan const & scan, MultiLayer const & sample)\n"
 		"new_SpecularSimulation(SpecularSimulation arg1) -> SpecularSimulation\n"
 		"SpecularSimulation::SpecularSimulation(SpecularSimulation &&)\n"
 		"\n"