diff --git a/Examples/Python/sim01_Particles/CylindersAndPrisms.py b/Examples/Python/sim01_Particles/CylindersAndPrisms.py
index cbab1b3b6d47c381c6f2adc6dc66d69f608dbf28..f6fe8987c66abbde0a2f7265205b92ee7d3ebe84 100644
--- a/Examples/Python/sim01_Particles/CylindersAndPrisms.py
+++ b/Examples/Python/sim01_Particles/CylindersAndPrisms.py
@@ -1,8 +1,9 @@
 """
 Mixture of cylinders and prisms without interference
 """
+import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, nm, nm2, kvector_t
+from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
 
 
 def get_sample():
@@ -48,22 +49,16 @@ def get_sample():
 
 
 def get_simulation():
-    """
-    Returns a GISAXS simulation with beam and detector defined.
-    """
-    simulation = ba.GISASSimulation()
-    simulation.setDetectorParameters(100, -1.0*deg, 1.0*deg, 100, 0.0*deg,
-                                     2.0*deg)
-    simulation.setBeamParameters(1.0*angstrom, 0.2*deg, 0.0*deg)
+    beam = ba.Beam(1.0, 0.1*nm, ba.Direction(0.2*deg, 0.0*deg))
+    detector = ba.SphericalDetector(100, -1.0*deg, 1.0*deg, 100, 0.0*deg,
+                                    2.0*deg)
+
+    simulation = ba.GISASSimulation(beam, get_sample(), detector)
     return simulation
 
 
 def run_simulation():
-    """
-    Runs simulation and returns resulting intensity map.
-    """
     simulation = get_simulation()
-    simulation.setSample(get_sample())
     simulation.runSimulation()
     return simulation.result()
 
diff --git a/Examples/Python/sim01_Particles/CylindersInBA.py b/Examples/Python/sim01_Particles/CylindersInBA.py
index 26855e44a93ba0f63acf6bb3e266b453afc3050d..50c91c0e3a87b4414c39677c4d72b215c2b1641c 100644
--- a/Examples/Python/sim01_Particles/CylindersInBA.py
+++ b/Examples/Python/sim01_Particles/CylindersInBA.py
@@ -1,8 +1,9 @@
 """
 Cylinder form factor in Born approximation
 """
+import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, nm, nm2, kvector_t
+from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
 
 
 def get_sample():
@@ -39,22 +40,16 @@ def get_sample():
 
 
 def get_simulation():
-    """
-    Returns a GISAXS simulation with beam and detector defined
-    """
-    simulation = ba.GISASSimulation()
-    simulation.setDetectorParameters(200, -2.0*deg, 2.0*deg, 200, 0.0*deg,
-                                     2.0*deg)
-    simulation.setBeamParameters(1.0*angstrom, 0.2*deg, 0.0*deg)
+    beam = ba.Beam(1.0, 0.1*nm, ba.Direction(0.2*deg, 0.0*deg))
+    detector = ba.SphericalDetector(200, -2.0*deg, 2.0*deg, 200, 0.0*deg,
+                                    2.0*deg)
+
+    simulation = ba.GISASSimulation(beam, get_sample(), detector)
     return simulation
 
 
 def run_simulation():
-    """
-    Runs simulation and returns intensity map.
-    """
     simulation = get_simulation()
-    simulation.setSample(get_sample())
     simulation.runSimulation()
     return simulation.result()
 
diff --git a/Examples/Python/sim01_Particles/CylindersInDWBA.py b/Examples/Python/sim01_Particles/CylindersInDWBA.py
index 28d884bf62df1e399244fdd48aba86bc76b148d0..97f141dc8699fb130d762551a34435ae8f645a62 100644
--- a/Examples/Python/sim01_Particles/CylindersInDWBA.py
+++ b/Examples/Python/sim01_Particles/CylindersInDWBA.py
@@ -1,8 +1,9 @@
 """
 Cylinder form factor in DWBA
 """
+import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, nm, nm2, kvector_t
+from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
 
 
 def get_sample():
@@ -41,22 +42,16 @@ def get_sample():
 
 
 def get_simulation():
-    """
-    Returns a GISAXS simulation with beam and detector defined
-    """
-    simulation = ba.GISASSimulation()
-    simulation.setDetectorParameters(200, -2.0*deg, 2.0*deg, 200, 0.0*deg,
-                                     2.0*deg)
-    simulation.setBeamParameters(1.0*angstrom, 0.2*deg, 0.0*deg)
+    beam = ba.Beam(1.0, 0.1*nm, ba.Direction(0.2*deg, 0.0*deg))
+    detector = ba.SphericalDetector(200, -2.0*deg, 2.0*deg, 200, 0.0*deg,
+                                    2.0*deg)
+
+    simulation = ba.GISASSimulation(beam, get_sample(), detector)
     return simulation
 
 
 def run_simulation():
-    """
-    Runs simulation and returns intensity map.
-    """
     simulation = get_simulation()
-    simulation.setSample(get_sample())
     simulation.runSimulation()
     return simulation.result()
 
diff --git a/Examples/Python/sim01_Particles/CylindersWithSizeDistribution.py b/Examples/Python/sim01_Particles/CylindersWithSizeDistribution.py
index adc1fd1decbcae25800c687551080346903e08d2..1407410ce6f5febb1caf1c0ea5a541b4e993af76 100644
--- a/Examples/Python/sim01_Particles/CylindersWithSizeDistribution.py
+++ b/Examples/Python/sim01_Particles/CylindersWithSizeDistribution.py
@@ -1,8 +1,9 @@
 """
 Cylinders with size distribution
 """
+import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, nm, nm2, kvector_t
+from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
 
 
 def get_sample():
@@ -45,22 +46,16 @@ def get_sample():
 
 
 def get_simulation():
-    """
-    Create and return GISAXS simulation with beam and detector defined
-    """
-    simulation = ba.GISASSimulation()
-    simulation.setDetectorParameters(200, 0.0*deg, 2.0*deg, 200, 0.0*deg,
-                                     2.0*deg)
-    simulation.setBeamParameters(1.0*angstrom, 0.2*deg, 0.0*deg)
+    beam = ba.Beam(1.0, 0.1*nm, ba.Direction(0.2*deg, 0.0*deg))
+    detector = ba.SphericalDetector(200, 0.0*deg, 2.0*deg, 200, 0.0*deg,
+                                    2.0*deg)
+
+    simulation = ba.GISASSimulation(beam, get_sample(), detector)
     return simulation
 
 
 def run_simulation():
-    """
-    Runs simulation and returns intensity map.
-    """
     simulation = get_simulation()
-    simulation.setSample(get_sample())
     simulation.runSimulation()
     return simulation.result()
 
diff --git a/Examples/Python/sim01_Particles/RotatedPyramids.py b/Examples/Python/sim01_Particles/RotatedPyramids.py
index 0077becb7bfe0fbb2de1ba8ed717227810b81960..a758f45164b54dc48c78964a70635fa52e59eb1d 100644
--- a/Examples/Python/sim01_Particles/RotatedPyramids.py
+++ b/Examples/Python/sim01_Particles/RotatedPyramids.py
@@ -1,8 +1,9 @@
 """
 Rotated pyramids on top of substrate
 """
+import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, nm, nm2, kvector_t
+from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
 
 
 def get_sample():
@@ -43,22 +44,16 @@ def get_sample():
 
 
 def get_simulation():
-    """
-    Returns a GISAXS simulation with beam and detector defined.
-    """
-    simulation = ba.GISASSimulation()
-    simulation.setDetectorParameters(200, -2.0*deg, 2.0*deg, 200, 0.0*deg,
-                                     2.0*deg)
-    simulation.setBeamParameters(1.0*angstrom, 0.2*deg, 0.0*deg)
+    beam = ba.Beam(1.0, 0.1*nm, ba.Direction(0.2*deg, 0.0*deg))
+    detector = ba.SphericalDetector(200, -2.0*deg, 2.0*deg, 200, 0.0*deg,
+                                    2.0*deg)
+
+    simulation = ba.GISASSimulation(beam, get_sample(), detector)
     return simulation
 
 
 def run_simulation():
-    """
-    Runs simulation and returns intensity map.
-    """
     simulation = get_simulation()
-    simulation.setSample(get_sample())
     simulation.runSimulation()
     return simulation.result()
 
diff --git a/Examples/Python/sim01_Particles/TwoTypesOfCylindersWithSizeDistribution.py b/Examples/Python/sim01_Particles/TwoTypesOfCylindersWithSizeDistribution.py
index 4609d50824a4e035e5a8a960b09a00ec1564754b..3bb2f240db48102dc12704f35a55c1424e48c8ad 100644
--- a/Examples/Python/sim01_Particles/TwoTypesOfCylindersWithSizeDistribution.py
+++ b/Examples/Python/sim01_Particles/TwoTypesOfCylindersWithSizeDistribution.py
@@ -1,8 +1,9 @@
 """
 Mixture cylinder particles with different size distribution
 """
+import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, nm, nm2, kvector_t
+from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
 
 
 def get_sample():
@@ -54,22 +55,16 @@ def get_sample():
 
 
 def get_simulation():
-    """
-    Create and return GISAXS simulation with beam and detector defined
-    """
-    simulation = ba.GISASSimulation()
-    simulation.setDetectorParameters(200, 0.0*deg, 2.0*deg, 200, 0.0*deg,
-                                     2.0*deg)
-    simulation.setBeamParameters(1.0*angstrom, 0.2*deg, 0.0*deg)
+    beam = ba.Beam(1.0, 0.1*nm, ba.Direction(0.2*deg, 0.0*deg))
+    detector = ba.SphericalDetector(200, 0.0*deg, 2.0*deg, 200, 0.0*deg,
+                                    2.0*deg)
+
+    simulation = ba.GISASSimulation(beam, get_sample(), detector)
     return simulation
 
 
 def run_simulation():
-    """
-    Runs simulation and returns intensity map.
-    """
     simulation = get_simulation()
-    simulation.setSample(get_sample())
     simulation.runSimulation()
     return simulation.result()
 
diff --git a/Examples/Python/sim02_Complexes/BiMaterialCylinders.py b/Examples/Python/sim02_Complexes/BiMaterialCylinders.py
index c99b9de98ba2a6c4bfefebcc54c150e46627ea53..14b770a15d0e67c7366c32180fa6b94f12be1695 100644
--- a/Examples/Python/sim02_Complexes/BiMaterialCylinders.py
+++ b/Examples/Python/sim02_Complexes/BiMaterialCylinders.py
@@ -2,30 +2,9 @@
 Cylindrical particle made from two materials.
 Particle crosses air/substrate interface.
 """
+import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, nm, nm2, kvector_t
-
-
-def get_composition(top_material,
-                    bottom_material,
-                    top_height=4.0,
-                    bottom_height=10.0):
-    """
-    Returns cylindrical particle made of two different materials.
-    """
-
-    cylinder_radius = 10*nm
-
-    topPart = ba.Particle(top_material,
-                          ba.FormFactorCylinder(cylinder_radius, top_height))
-    bottomPart = ba.Particle(
-        bottom_material, ba.FormFactorCylinder(cylinder_radius, bottom_height))
-
-    result = ba.ParticleComposition()
-    result.addParticle(topPart, ba.kvector_t(0.0, 0.0, bottom_height))
-    result.addParticle(bottomPart)
-
-    return result
+from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
 
 
 def get_sample():
@@ -79,23 +58,16 @@ def get_sample():
 
 
 def get_simulation():
-    """
-    Returns a GISAXS simulation with beam and detector defined.
-    """
-    simulation = ba.GISASSimulation()
-    simulation.setDetectorParameters(100, -1.0*deg, 1.0*deg, 100, 0.0*deg,
-                                     2.0*deg)
-    simulation.setBeamParameters(1.0*angstrom, 0.2*deg, 0.0*deg)
-    simulation.beam().setIntensity(1.0e+08)
+    beam = ba.Beam(100000000.0, 0.1*nm, ba.Direction(0.2*deg, 0.0*deg))
+    detector = ba.SphericalDetector(100, -1.0*deg, 1.0*deg, 100, 0.0*deg,
+                                    2.0*deg)
+
+    simulation = ba.GISASSimulation(beam, get_sample(), detector)
     return simulation
 
 
 def run_simulation():
-    """
-    Runs simulation and returns resulting intensity map.
-    """
     simulation = get_simulation()
-    simulation.setSample(get_sample())
     simulation.runSimulation()
     return simulation.result()
 
diff --git a/Examples/Python/sim02_Complexes/CoreShellNanoparticles.py b/Examples/Python/sim02_Complexes/CoreShellNanoparticles.py
index 0ff5326cf1f6ba7d591fdfa20063ec2540d3fde2..8fab24ca4b6c592a6b33a55f883cb713a4852e55 100644
--- a/Examples/Python/sim02_Complexes/CoreShellNanoparticles.py
+++ b/Examples/Python/sim02_Complexes/CoreShellNanoparticles.py
@@ -1,8 +1,9 @@
 """
 Core shell nanoparticles
 """
+import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, nm, nm2, kvector_t
+from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
 
 
 def get_sample():
@@ -48,22 +49,16 @@ def get_sample():
 
 
 def get_simulation():
-    """
-    Returns a GISAXS simulation with beam and detector defined.
-    """
-    simulation = ba.GISASSimulation()
-    simulation.setDetectorParameters(200, -1.0*deg, 1.0*deg, 200, 0.0*deg,
-                                     2.0*deg)
-    simulation.setBeamParameters(1.0*angstrom, 0.2*deg, 0.0*deg)
+    beam = ba.Beam(1.0, 0.1*nm, ba.Direction(0.2*deg, 0.0*deg))
+    detector = ba.SphericalDetector(200, -1.0*deg, 1.0*deg, 200, 0.0*deg,
+                                    2.0*deg)
+
+    simulation = ba.GISASSimulation(beam, get_sample(), detector)
     return simulation
 
 
 def run_simulation():
-    """
-    Runs simulation and returns intensity map.
-    """
     simulation = get_simulation()
-    simulation.setSample(get_sample())
     simulation.runSimulation()
     return simulation.result()
 
diff --git a/Examples/Python/sim02_Complexes/HexagonalLatticesWithBasis.py b/Examples/Python/sim02_Complexes/HexagonalLatticesWithBasis.py
index 860732d97d96c21cdd8a05304edf2b4bc84b9855..25aeb82c29fb2d67705aa9c0878179ccde8d4711 100644
--- a/Examples/Python/sim02_Complexes/HexagonalLatticesWithBasis.py
+++ b/Examples/Python/sim02_Complexes/HexagonalLatticesWithBasis.py
@@ -1,9 +1,9 @@
 """
 Spheres on two hexagonal close packed layers
 """
-import numpy
+import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, nm, nm2, kvector_t
+from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
 
 
 def get_sample():
@@ -61,22 +61,16 @@ def get_sample():
 
 
 def get_simulation():
-    """
-    Returns a GISAXS simulation with beam and detector defined.
-    """
-    simulation = ba.GISASSimulation()
-    simulation.setDetectorParameters(200, -1.0*deg, 1.0*deg, 200, 0.0*deg,
-                                     1.0*deg)
-    simulation.setBeamParameters(1.0*angstrom, 0.2*deg, 0.0*deg)
+    beam = ba.Beam(1.0, 0.1*nm, ba.Direction(0.2*deg, 0.0*deg))
+    detector = ba.SphericalDetector(200, -1.0*deg, 1.0*deg, 200, 0.0*deg,
+                                    1.0*deg)
+
+    simulation = ba.GISASSimulation(beam, get_sample(), detector)
     return simulation
 
 
 def run_simulation():
-    """
-    Runs simulation and returns intensity map.
-    """
     simulation = get_simulation()
-    simulation.setSample(get_sample())
     simulation.runSimulation()
     return simulation.result()
 
diff --git a/Examples/Python/sim02_Complexes/MesoCrystal.py b/Examples/Python/sim02_Complexes/MesoCrystal.py
index 3bd203daadd04ab6d752e726b4b8dc616ddf2e31..de317479061d9f6afbcb586fbc51f214a381ced5 100644
--- a/Examples/Python/sim02_Complexes/MesoCrystal.py
+++ b/Examples/Python/sim02_Complexes/MesoCrystal.py
@@ -1,8 +1,9 @@
 """
 Cylindrical mesocrystal on a substrate
 """
+import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, nm, nm2, kvector_t
+from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
 
 
 def get_sample():
@@ -53,22 +54,16 @@ def get_sample():
 
 
 def get_simulation():
-    """
-    Returns a GISAXS simulation with beam and detector defined
-    """
-    simulation = ba.GISASSimulation()
-    simulation.setDetectorParameters(200, -2.0*deg, 2.0*deg, 200, 0.0*deg,
-                                     2.0*deg)
-    simulation.setBeamParameters(1.0*angstrom, 0.2*deg, 0.0*deg)
+    beam = ba.Beam(1.0, 0.1*nm, ba.Direction(0.2*deg, 0.0*deg))
+    detector = ba.SphericalDetector(200, -2.0*deg, 2.0*deg, 200, 0.0*deg,
+                                    2.0*deg)
+
+    simulation = ba.GISASSimulation(beam, get_sample(), detector)
     return simulation
 
 
 def run_simulation():
-    """
-    Runs simulation and returns intensity map.
-    """
     simulation = get_simulation()
-    simulation.setSample(get_sample())
     simulation.runSimulation()
     return simulation.result()
 
diff --git a/Examples/Python/sim02_Complexes/ParticlesCrossingInterface.py b/Examples/Python/sim02_Complexes/ParticlesCrossingInterface.py
index 4065dbd16b0bebe16f00d52bd210fcf489f372e0..43fe5da2928b0b85ff45e924e4a41b9c9c1f71c7 100644
--- a/Examples/Python/sim02_Complexes/ParticlesCrossingInterface.py
+++ b/Examples/Python/sim02_Complexes/ParticlesCrossingInterface.py
@@ -16,11 +16,9 @@ adjusts calculations accordingly.
 For example, X or Y rotated particles can not yet cross interfaces (exception
 will be thrown when trying to simulate such geometries).
 """
+import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, nm, nm2, kvector_t
-
-phi_min, phi_max = -1.0, 1.0
-alpha_min, alpha_max = 0.0, 2.0
+from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
 
 
 def get_sample():
@@ -70,23 +68,16 @@ def get_sample():
 
 
 def get_simulation():
-    """
-    Returns a GISAXS simulation with beam and detector defined.
-    """
-    simulation = ba.GISASSimulation()
-    simulation.setDetectorParameters(100, phi_min*deg, phi_max*deg, 100,
-                                     alpha_min*deg, alpha_max*deg)
-    simulation.setBeamParameters(1.0*angstrom, 0.2*deg, 0.0*deg)
+    beam = ba.Beam(1.0, 0.1*nm, ba.Direction(0.2*deg, 0.0*deg))
+    detector = ba.SphericalDetector(100, -1.0*deg, 1.0*deg, 100, 0.0*deg,
+                                    2.0*deg)
+
+    simulation = ba.GISASSimulation(beam, get_sample(), detector)
     return simulation
 
 
 def run_simulation():
-    """
-    Runs simulation and returns resulting intensity map.
-    """
-    sample = get_sample()
     simulation = get_simulation()
-    simulation.setSample(sample)
     simulation.runSimulation()
     return simulation.result()
 
diff --git a/Examples/Python/sim03_Structures/ApproximationDA.py b/Examples/Python/sim03_Structures/ApproximationDA.py
index c072968b804ba59f5effe416b447e620af14c038..2af0151cdde1678614eecb9f8e4e546a99a0d557 100644
--- a/Examples/Python/sim03_Structures/ApproximationDA.py
+++ b/Examples/Python/sim03_Structures/ApproximationDA.py
@@ -1,8 +1,9 @@
 """
 Cylinders of two different sizes in Decoupling Approximation
 """
+import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, nm, nm2, kvector_t
+from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
 
 
 def get_sample():
@@ -51,22 +52,16 @@ def get_sample():
 
 
 def get_simulation():
-    """
-    Create and return GISAXS simulation with beam and detector defined
-    """
-    simulation = ba.GISASSimulation()
-    simulation.setDetectorParameters(200, 0.0*deg, 2.0*deg, 200, 0.0*deg,
-                                     2.0*deg)
-    simulation.setBeamParameters(1.0*angstrom, 0.2*deg, 0.0*deg)
+    beam = ba.Beam(1.0, 0.1*nm, ba.Direction(0.2*deg, 0.0*deg))
+    detector = ba.SphericalDetector(200, 0.0*deg, 2.0*deg, 200, 0.0*deg,
+                                    2.0*deg)
+
+    simulation = ba.GISASSimulation(beam, get_sample(), detector)
     return simulation
 
 
 def run_simulation():
-    """
-    Runs simulation and returns intensity map.
-    """
     simulation = get_simulation()
-    simulation.setSample(get_sample())
     simulation.runSimulation()
     return simulation.result()
 
diff --git a/Examples/Python/sim03_Structures/ApproximationLMA.py b/Examples/Python/sim03_Structures/ApproximationLMA.py
index 73c8eebb9b3947f23b9da1de47a4608d8b17e872..b054a61ababdc69656c2de5c8ca62e773c413705 100644
--- a/Examples/Python/sim03_Structures/ApproximationLMA.py
+++ b/Examples/Python/sim03_Structures/ApproximationLMA.py
@@ -1,8 +1,9 @@
 """
 Cylinders of two different sizes in Local Monodisperse Approximation
 """
+import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, nm, nm2, kvector_t
+from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
 
 
 def get_sample():
@@ -59,22 +60,16 @@ def get_sample():
 
 
 def get_simulation():
-    """
-    Create and return GISAXS simulation with beam and detector defined
-    """
-    simulation = ba.GISASSimulation()
-    simulation.setDetectorParameters(200, 0.0*deg, 2.0*deg, 200, 0.0*deg,
-                                     2.0*deg)
-    simulation.setBeamParameters(1.0*angstrom, 0.2*deg, 0.0*deg)
+    beam = ba.Beam(1.0, 0.1*nm, ba.Direction(0.2*deg, 0.0*deg))
+    detector = ba.SphericalDetector(200, 0.0*deg, 2.0*deg, 200, 0.0*deg,
+                                    2.0*deg)
+
+    simulation = ba.GISASSimulation(beam, get_sample(), detector)
     return simulation
 
 
 def run_simulation():
-    """
-    Runs simulation and returns intensity map.
-    """
     simulation = get_simulation()
-    simulation.setSample(get_sample())
     simulation.runSimulation()
     return simulation.result()
 
diff --git a/Examples/Python/sim03_Structures/ApproximationSSCA.py b/Examples/Python/sim03_Structures/ApproximationSSCA.py
index abc4a8cb83ba6dd9afedf41a54b0aa6ac54e53cf..f8346585db469dab6d6a7aec867db33bdfc5398a 100644
--- a/Examples/Python/sim03_Structures/ApproximationSSCA.py
+++ b/Examples/Python/sim03_Structures/ApproximationSSCA.py
@@ -1,8 +1,9 @@
 """
 Cylinders of two different sizes in Size-Spacing Coupling Approximation
 """
+import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, nm, nm2, kvector_t
+from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
 
 
 def get_sample():
@@ -52,22 +53,16 @@ def get_sample():
 
 
 def get_simulation():
-    """
-    Create and return GISAXS simulation with beam and detector defined
-    """
-    simulation = ba.GISASSimulation()
-    simulation.setDetectorParameters(200, 0.0*deg, 2.0*deg, 200, 0.0*deg,
-                                     2.0*deg)
-    simulation.setBeamParameters(1.0*angstrom, 0.2*deg, 0.0*deg)
+    beam = ba.Beam(1.0, 0.1*nm, ba.Direction(0.2*deg, 0.0*deg))
+    detector = ba.SphericalDetector(200, 0.0*deg, 2.0*deg, 200, 0.0*deg,
+                                    2.0*deg)
+
+    simulation = ba.GISASSimulation(beam, get_sample(), detector)
     return simulation
 
 
 def run_simulation():
-    """
-    Runs simulation and returns intensity map.
-    """
     simulation = get_simulation()
-    simulation.setSample(get_sample())
     simulation.runSimulation()
     return simulation.result()
 
diff --git a/Examples/Python/sim03_Structures/CosineRipplesAtRectLattice.py b/Examples/Python/sim03_Structures/CosineRipplesAtRectLattice.py
index b4ae226890224b83dc62439efe2f945e434d0c20..b14ee49a5cc64696a7ca30691d08da38d5d2397f 100644
--- a/Examples/Python/sim03_Structures/CosineRipplesAtRectLattice.py
+++ b/Examples/Python/sim03_Structures/CosineRipplesAtRectLattice.py
@@ -1,9 +1,9 @@
 """
 Cosine ripple on a 2D lattice
 """
-import numpy
+import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, nm, nm2, kvector_t
+from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
 
 
 def get_sample():
@@ -52,22 +52,16 @@ def get_sample():
 
 
 def get_simulation():
-    """
-    characterizing the input beam and output detector
-    """
-    simulation = ba.GISASSimulation()
-    simulation.setDetectorParameters(100, -1.5*deg, 1.5*deg, 100, 0.0*deg,
-                                     2.5*deg)
-    simulation.setBeamParameters(1.6*angstrom, 0.3*deg, 0.0*deg)
+    beam = ba.Beam(1.0, 0.16*nm, ba.Direction(0.3*deg, 0.0*deg))
+    detector = ba.SphericalDetector(100, -1.5*deg, 1.5*deg, 100, 0.0*deg,
+                                    2.5*deg)
+
+    simulation = ba.GISASSimulation(beam, get_sample(), detector)
     return simulation
 
 
 def run_simulation():
-    """
-    Runs simulation and returns intensity map.
-    """
     simulation = get_simulation()
-    simulation.setSample(get_sample())
     simulation.runSimulation()
     return simulation.result()
 
diff --git a/Examples/Python/sim03_Structures/Interference1DRadialParaCrystal.py b/Examples/Python/sim03_Structures/Interference1DRadialParaCrystal.py
index 12e550a5a155149517eeca0327417769bcc95125..29281c48b85420952d7c40467be19a3f5baee079 100644
--- a/Examples/Python/sim03_Structures/Interference1DRadialParaCrystal.py
+++ b/Examples/Python/sim03_Structures/Interference1DRadialParaCrystal.py
@@ -1,11 +1,9 @@
 """
 radial paracrystal
 """
+import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, nm, nm2, kvector_t
-
-phi_min, phi_max = -2.0, 2.0
-alpha_min, alpha_max = 0.0, 2.0
+from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
 
 
 def get_sample():
@@ -50,22 +48,16 @@ def get_sample():
 
 
 def get_simulation():
-    """
-    Create and return GISAXS simulation with beam and detector defined
-    """
-    simulation = ba.GISASSimulation()
-    simulation.setDetectorParameters(200, -2.0*deg, 2.0*deg, 200, 0.0*deg,
-                                     2.0*deg)
-    simulation.setBeamParameters(1.0*angstrom, 0.2*deg, 0.0*deg)
+    beam = ba.Beam(1.0, 0.1*nm, ba.Direction(0.2*deg, 0.0*deg))
+    detector = ba.SphericalDetector(200, -2.0*deg, 2.0*deg, 200, 0.0*deg,
+                                    2.0*deg)
+
+    simulation = ba.GISASSimulation(beam, get_sample(), detector)
     return simulation
 
 
 def run_simulation():
-    """
-    Runs simulation and returns intensity map.
-    """
     simulation = get_simulation()
-    simulation.setSample(get_sample())
     simulation.runSimulation()
     return simulation.result()
 
diff --git a/Examples/Python/sim03_Structures/Interference2DCenteredSquareLattice.py b/Examples/Python/sim03_Structures/Interference2DCenteredSquareLattice.py
index 03049b3fb7b280223d910daf2450e1c90088194a..696953b7f7670a1c85e35cfbdcff1a251ab3d2cb 100644
--- a/Examples/Python/sim03_Structures/Interference2DCenteredSquareLattice.py
+++ b/Examples/Python/sim03_Structures/Interference2DCenteredSquareLattice.py
@@ -1,9 +1,9 @@
 """
 2D lattice with disorder, centered square lattice
 """
-import numpy
+import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, nm, nm2, kvector_t
+from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
 
 
 def get_sample():
@@ -61,23 +61,16 @@ def get_sample():
 
 
 def get_simulation():
-    """
-    Create and return GISAS simulation with beam and detector defined
-    """
-    simulation = ba.GISASSimulation()
-    simulation.setDetectorParameters(200, -2.0*deg, 2.0*deg, 200, 0.0*deg,
-                                     2.0*deg)
-    simulation.setBeamParameters(1.0*angstrom, 0.2*deg, 0.0*deg)
+    beam = ba.Beam(1.0, 0.1*nm, ba.Direction(0.2*deg, 0.0*deg))
+    detector = ba.SphericalDetector(200, -2.0*deg, 2.0*deg, 200, 0.0*deg,
+                                    2.0*deg)
 
+    simulation = ba.GISASSimulation(beam, get_sample(), detector)
     return simulation
 
 
 def run_simulation():
-    """
-    Runs simulation and returns intensity map.
-    """
     simulation = get_simulation()
-    simulation.setSample(get_sample())
     simulation.runSimulation()
     return simulation.result()
 
diff --git a/Examples/Python/sim03_Structures/Interference2DLatticeSumOfRotated.py b/Examples/Python/sim03_Structures/Interference2DLatticeSumOfRotated.py
index 1948f822ed174e40c629bbea026678b28806d920..7ef1948867e45ef0a95861c0c4237a8b1d0534c7 100644
--- a/Examples/Python/sim03_Structures/Interference2DLatticeSumOfRotated.py
+++ b/Examples/Python/sim03_Structures/Interference2DLatticeSumOfRotated.py
@@ -1,7 +1,6 @@
-# 2D lattice with different disorder (IsGISAXS example #6), sum of rotated lattices
-import numpy
+import numpy, sys
 import bornagain as ba
-from bornagain import deg, angstrom, nm
+from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
 
 
 def get_sample():
@@ -38,32 +37,17 @@ def get_sample():
 
 
 def get_simulation():
-    """
-    Returns a GISAXS simulation with beam and detector defined.
-    Assigns additional distribution to lattice rotation angle.
-    """
-    simulation = ba.GISASSimulation()
-    simulation.setDetectorParameters(100, 0.0*deg, 2.0*deg, 100, 0.0*deg,
-                                     2.0*deg)
-    simulation.setBeamParameters(1.0*angstrom, 0.2*deg, 0.0*deg)
-
-    simulation.setSample(get_sample())
-
-    xi_min = 0.0*deg
-    xi_max = 240.0*deg
-    xi_distr = ba.DistributionGate(xi_min, xi_max)
-
-    # assigns distribution with 3 equidistant points to lattice rotation angle
-    simulation.addParameterDistribution("*/SquareLattice2D/Xi", xi_distr, 3)
+    beam = ba.Beam(1.0, 0.1*nm, ba.Direction(0.2*deg, 0.0*deg))
+    detector = ba.SphericalDetector(100, 0.0*deg, 2.0*deg, 100, 0.0*deg,
+                                    2.0*deg)
 
+    simulation = ba.GISASSimulation(beam, get_sample(), detector)
+    distr_1 = ba.DistributionGate(0.0*deg, 240.0*deg)
+    simulation.addParameterDistribution("*/SquareLattice2D/Xi", distr_1, 3, 0.0)
     return simulation
 
 
 def run_simulation():
-    """
-    Runs simulation and returns intensity map.
-    """
-
     simulation = get_simulation()
     simulation.runSimulation()
     return simulation.result()
diff --git a/Examples/Python/sim03_Structures/Interference2DParaCrystal.py b/Examples/Python/sim03_Structures/Interference2DParaCrystal.py
index dfae0d5c546f0ba8c4e51376d145db064bdf33d3..c6106e79db38766efafd5dee52f485c07e4a5e47 100644
--- a/Examples/Python/sim03_Structures/Interference2DParaCrystal.py
+++ b/Examples/Python/sim03_Structures/Interference2DParaCrystal.py
@@ -1,8 +1,9 @@
 """
 2D paracrystal
 """
+import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, nm, nm2, kvector_t
+from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
 
 
 def get_sample():
@@ -53,25 +54,16 @@ def get_sample():
 
 
 def get_simulation():
-    """
-    Returns a GISAXS simulation with beam and detector defined.
-    """
-    simulation = ba.GISASSimulation()
-    # coarse grid because this simulation takes rather long
-    simulation.setDetectorParameters(200, -2.0*deg, 2.0*deg, 200, 0.0*deg,
-                                     2.0*deg)
-    simulation.setBeamParameters(1.0*angstrom, 0.2*deg, 0.0*deg)
+    beam = ba.Beam(1.0, 0.1*nm, ba.Direction(0.2*deg, 0.0*deg))
+    detector = ba.SphericalDetector(200, -2.0*deg, 2.0*deg, 200, 0.0*deg,
+                                    2.0*deg)
+
+    simulation = ba.GISASSimulation(beam, get_sample(), detector)
     return simulation
 
 
 def run_simulation():
-    """
-    Runs simulation and returns intensity map.
-    """
     simulation = get_simulation()
-    simulation.setSample(get_sample())
-    if not "__no_terminal__" in globals():
-        simulation.setTerminalProgressMonitor()
     simulation.runSimulation()
     return simulation.result()
 
diff --git a/Examples/Python/sim03_Structures/Interference2DRotatedSquareLattice.py b/Examples/Python/sim03_Structures/Interference2DRotatedSquareLattice.py
index 33ec3f79809db7ed9ef3666054966b6723863b70..bc6f76a91bf6e414b6b8230ecadba63dda5245bc 100644
--- a/Examples/Python/sim03_Structures/Interference2DRotatedSquareLattice.py
+++ b/Examples/Python/sim03_Structures/Interference2DRotatedSquareLattice.py
@@ -1,9 +1,9 @@
 """
 Cylinders on a rotated 2D lattice
 """
-import numpy
+import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, nm, nm2, kvector_t
+from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
 
 
 def get_sample():
@@ -52,23 +52,16 @@ def get_sample():
 
 
 def get_simulation():
-    """
-    Create and return GISAXS simulation with beam and detector defined
-    """
-    simulation = ba.GISASSimulation()
-    simulation.setDetectorParameters(200, -2.0*deg, 2.0*deg, 200, 0.0*deg,
-                                     2.0*deg)
-    simulation.setBeamParameters(1.0*angstrom, 0.2*deg, 0.0*deg)
+    beam = ba.Beam(1.0, 0.1*nm, ba.Direction(0.2*deg, 0.0*deg))
+    detector = ba.SphericalDetector(200, -2.0*deg, 2.0*deg, 200, 0.0*deg,
+                                    2.0*deg)
 
+    simulation = ba.GISASSimulation(beam, get_sample(), detector)
     return simulation
 
 
 def run_simulation():
-    """
-    Runs simulation and returns intensity map.
-    """
     simulation = get_simulation()
-    simulation.setSample(get_sample())
     simulation.runSimulation()
     return simulation.result()
 
diff --git a/Examples/Python/sim03_Structures/Interference2DSquareFiniteLattice.py b/Examples/Python/sim03_Structures/Interference2DSquareFiniteLattice.py
index e018b063c317bb9935640c8dce3f10b1fd4ffff5..5fee0b8c214dade6e3e00a74ab6f60fd0d58fec9 100644
--- a/Examples/Python/sim03_Structures/Interference2DSquareFiniteLattice.py
+++ b/Examples/Python/sim03_Structures/Interference2DSquareFiniteLattice.py
@@ -1,9 +1,9 @@
 """
 Cylinders on a 2D square lattice
 """
-import numpy
+import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, nm, nm2, kvector_t
+from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
 
 
 def get_sample():
@@ -50,22 +50,16 @@ def get_sample():
 
 
 def get_simulation():
-    """
-    Create and return GISAXS simulation with beam and detector defined
-    """
-    simulation = ba.GISASSimulation()
-    simulation.setDetectorParameters(200, -2.0*deg, 2.0*deg, 200, 0.0*deg,
-                                     2.0*deg)
-    simulation.setBeamParameters(1.0*angstrom, 0.2*deg, 0.0*deg)
+    beam = ba.Beam(1.0, 0.1*nm, ba.Direction(0.2*deg, 0.0*deg))
+    detector = ba.SphericalDetector(200, -2.0*deg, 2.0*deg, 200, 0.0*deg,
+                                    2.0*deg)
+
+    simulation = ba.GISASSimulation(beam, get_sample(), detector)
     return simulation
 
 
 def run_simulation():
-    """
-    Runs simulation and returns intensity map.
-    """
     simulation = get_simulation()
-    simulation.setSample(get_sample())
     simulation.runSimulation()
     return simulation.result()
 
diff --git a/Examples/Python/sim03_Structures/Interference2DSquareLattice.py b/Examples/Python/sim03_Structures/Interference2DSquareLattice.py
index affd7715ec59e66c594c4acf5a22332d0da8d990..8ff54b2eb00713be38cd6d49b2e1b550b8bfe815 100644
--- a/Examples/Python/sim03_Structures/Interference2DSquareLattice.py
+++ b/Examples/Python/sim03_Structures/Interference2DSquareLattice.py
@@ -1,9 +1,9 @@
 """
 Cylinders on a 2D square lattice
 """
-import numpy
+import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, nm, nm2, kvector_t
+from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
 
 
 def get_sample():
@@ -51,22 +51,16 @@ def get_sample():
 
 
 def get_simulation():
-    """
-    Create and return GISAXS simulation with beam and detector defined
-    """
-    simulation = ba.GISASSimulation()
-    simulation.setDetectorParameters(200, -2.0*deg, 2.0*deg, 200, 0.0*deg,
-                                     2.0*deg)
-    simulation.setBeamParameters(1.0*angstrom, 0.2*deg, 0.0*deg)
+    beam = ba.Beam(1.0, 0.1*nm, ba.Direction(0.2*deg, 0.0*deg))
+    detector = ba.SphericalDetector(200, -2.0*deg, 2.0*deg, 200, 0.0*deg,
+                                    2.0*deg)
+
+    simulation = ba.GISASSimulation(beam, get_sample(), detector)
     return simulation
 
 
 def run_simulation():
-    """
-    Runs simulation and returns intensity map.
-    """
     simulation = get_simulation()
-    simulation.setSample(get_sample())
     simulation.runSimulation()
     return simulation.result()
 
diff --git a/Examples/Python/sim03_Structures/RectangularGrating.py b/Examples/Python/sim03_Structures/RectangularGrating.py
index da5b67337499317af369d351c39d7c99a1c57af4..d1f6722bd3726d45785c56c9460f71842ee09d26 100644
--- a/Examples/Python/sim03_Structures/RectangularGrating.py
+++ b/Examples/Python/sim03_Structures/RectangularGrating.py
@@ -3,8 +3,9 @@ Simulation of grating using very long boxes and 1D lattice.
 Monte-carlo integration is used to get rid of
 large-particle form factor oscillations.
 """
+import numpy, sys
 import bornagain as ba
-from bornagain import deg, angstrom, nm, micrometer
+from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
 
 
 def get_sample(lattice_rotation_angle=0.0*deg):
@@ -52,26 +53,17 @@ def get_sample(lattice_rotation_angle=0.0*deg):
 
 
 def get_simulation():
-    """
-    Create and return GISAXS simulation with beam and detector defined
-    """
-    simulation = ba.GISASSimulation()
-    simulation.setDetectorParameters(200, -0.5*deg, 0.5*deg, 200, 0.0*deg,
-                                     0.6*deg)
-    simulation.setBeamParameters(1.34*angstrom, 0.4*deg, 0.0*deg)
-    simulation.beam().setIntensity(1e+08)
+    beam = ba.Beam(100000000.0, 0.134*nm, ba.Direction(0.4*deg, 0.0*deg))
+    detector = ba.SphericalDetector(200, -0.5*deg, 0.5*deg, 200, 0.0*deg,
+                                    0.6*deg)
+
+    simulation = ba.GISASSimulation(beam, get_sample(), detector)
     simulation.getOptions().setMonteCarloIntegration(True, 100)
     return simulation
 
 
 def run_simulation():
-    """
-    Runs simulation and returns intensity map.
-    """
     simulation = get_simulation()
-    simulation.setSample(get_sample())
-    if not "__no_terminal__" in globals():
-        simulation.setTerminalProgressMonitor()
     simulation.runSimulation()
     return simulation.result()
 
diff --git a/Examples/Python/sim03_Structures/SpheresAtHexLattice.py b/Examples/Python/sim03_Structures/SpheresAtHexLattice.py
index 7f3c49d042b4f7e05fa5b942c17900b3816cb40e..8318bd33c930ba2ea8fa30cd976e8cf79c9c166f 100644
--- a/Examples/Python/sim03_Structures/SpheresAtHexLattice.py
+++ b/Examples/Python/sim03_Structures/SpheresAtHexLattice.py
@@ -1,8 +1,9 @@
 """
 Spheres on a hexagonal lattice
 """
+import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, nm, nm2, kvector_t
+from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
 
 
 def get_sample():
@@ -51,22 +52,16 @@ def get_sample():
 
 
 def get_simulation():
-    """
-    Create and return GISAXS simulation with beam and detector defined
-    """
-    simulation = ba.GISASSimulation()
-    simulation.setDetectorParameters(200, -1.0*deg, 1.0*deg, 200, 0.0*deg,
-                                     1.0*deg)
-    simulation.setBeamParameters(1.0*angstrom, 0.2*deg, 0.0*deg)
+    beam = ba.Beam(1.0, 0.1*nm, ba.Direction(0.2*deg, 0.0*deg))
+    detector = ba.SphericalDetector(200, -1.0*deg, 1.0*deg, 200, 0.0*deg,
+                                    1.0*deg)
+
+    simulation = ba.GISASSimulation(beam, get_sample(), detector)
     return simulation
 
 
 def run_simulation():
-    """
-    Runs simulation and returns intensity map.
-    """
     simulation = get_simulation()
-    simulation.setSample(get_sample())
     simulation.runSimulation()
     return simulation.result()
 
diff --git a/Examples/Python/sim03_Structures/TriangularRipple.py b/Examples/Python/sim03_Structures/TriangularRipple.py
index 2815287042a10a22f4e7cb04849b1132650fd8b9..d89e7ecae0ffadf827dc37d391e557d314803fd6 100644
--- a/Examples/Python/sim03_Structures/TriangularRipple.py
+++ b/Examples/Python/sim03_Structures/TriangularRipple.py
@@ -1,9 +1,9 @@
 """
  Sample from the article D. Babonneau et. al., Phys. Rev. B 85, 235415, 2012 (Fig.3)
 """
-import numpy
+import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, nm, nm2, kvector_t
+from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
 
 
 def get_sample():
@@ -53,22 +53,16 @@ def get_sample():
 
 
 def get_simulation():
-    """
-    characterizing the input beam and output detector
-    """
-    simulation = ba.GISASSimulation()
-    simulation.setDetectorParameters(200, -1.5*deg, 1.5*deg, 200, 0.0*deg,
-                                     2.5*deg)
-    simulation.setBeamParameters(1.6*angstrom, 0.3*deg, 0.0*deg)
+    beam = ba.Beam(1.0, 0.16*nm, ba.Direction(0.3*deg, 0.0*deg))
+    detector = ba.SphericalDetector(200, -1.5*deg, 1.5*deg, 200, 0.0*deg,
+                                    2.5*deg)
+
+    simulation = ba.GISASSimulation(beam, get_sample(), detector)
     return simulation
 
 
 def run_simulation():
-    """
-    Runs simulation and returns intensity map.
-    """
     simulation = get_simulation()
-    simulation.setSample(get_sample())
     simulation.runSimulation()
     return simulation.result()
 
diff --git a/Examples/Python/sim04_Multilayers/BuriedParticles.py b/Examples/Python/sim04_Multilayers/BuriedParticles.py
index d543571da20b84d28b480c3c3dcc2cec353c3d3c..0d3dfe0e880832dc83a9591279991ccbffd2f938 100644
--- a/Examples/Python/sim04_Multilayers/BuriedParticles.py
+++ b/Examples/Python/sim04_Multilayers/BuriedParticles.py
@@ -1,8 +1,9 @@
 """
 Spherical particles embedded in the middle of the layer on top of substrate.
 """
+import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, nm, nm2, kvector_t
+from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
 
 
 def get_sample():
@@ -48,22 +49,16 @@ def get_sample():
 
 
 def get_simulation():
-    """
-    Returns a GISAXS simulation.
-    """
-    simulation = ba.GISASSimulation()
-    simulation.setSample(get_sample())
-    simulation.setDetectorParameters(200, -1*deg, +1*deg, 200, 0*deg, +2*deg)
-    simulation.setBeamParameters(1.5*angstrom, 0.15*deg, 0.0*deg)
+    beam = ba.Beam(1.0, 0.15*nm, ba.Direction(0.15*deg, 0.0*deg))
+    detector = ba.SphericalDetector(200, -1.0*deg, 1.0*deg, 200, 0.0*deg,
+                                    2.0*deg)
+
+    simulation = ba.GISASSimulation(beam, get_sample(), detector)
     return simulation
 
 
 def run_simulation():
-    """
-    Runs simulation and returns intensity map.
-    """
     simulation = get_simulation()
-    simulation.setSample(get_sample())
     simulation.runSimulation()
     return simulation.result()
 
diff --git a/Examples/Python/sim04_Multilayers/CorrelatedRoughness.py b/Examples/Python/sim04_Multilayers/CorrelatedRoughness.py
index a32410384e7c997082411d3a6dbf060feee46284..32b9d86210e6ea9d3c765ea02c8544099f0c42ac 100644
--- a/Examples/Python/sim04_Multilayers/CorrelatedRoughness.py
+++ b/Examples/Python/sim04_Multilayers/CorrelatedRoughness.py
@@ -1,8 +1,9 @@
 """
 MultiLayer with correlated roughness
 """
+import numpy, sys
 import bornagain as ba
-from bornagain import deg, angstrom, nm
+from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
 
 
 def get_sample():
@@ -45,23 +46,16 @@ def get_sample():
 
 
 def get_simulation():
-    """
-    Characterizing the input beam and output detector
-    """
-    simulation = ba.GISASSimulation()
-    simulation.setDetectorParameters(200, -0.5*deg, 0.5*deg, 200, 0.0*deg,
-                                     1.0*deg)
-    simulation.setBeamParameters(1.0*angstrom, 0.2*deg, 0.0*deg)
-    simulation.beam().setIntensity(5e11)
+    beam = ba.Beam(500000000000.0, 0.1*nm, ba.Direction(0.2*deg, 0.0*deg))
+    detector = ba.SphericalDetector(200, -0.5*deg, 0.5*deg, 200, 0.0*deg,
+                                    1.0*deg)
+
+    simulation = ba.GISASSimulation(beam, get_sample(), detector)
     return simulation
 
 
 def run_simulation():
-    """
-    Runs simulation and returns intensity map.
-    """
     simulation = get_simulation()
-    simulation.setSample(get_sample())
     simulation.runSimulation()
     return simulation.result()
 
diff --git a/Examples/Python/sim04_Multilayers/CylindersInAverageLayer.py b/Examples/Python/sim04_Multilayers/CylindersInAverageLayer.py
index bfbe4d269e5e858c40f11980d172608c49d41ac1..e8c31fe3efc91a4d7e35c56af67a9b8c92a20346 100644
--- a/Examples/Python/sim04_Multilayers/CylindersInAverageLayer.py
+++ b/Examples/Python/sim04_Multilayers/CylindersInAverageLayer.py
@@ -1,8 +1,9 @@
 """
 Square lattice of cylinders inside finite layer with usage of average material
 """
+import numpy, sys
 import bornagain as ba
-from bornagain import deg, angstrom, nm
+from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
 
 
 def get_sample(cyl_height=5*nm):
@@ -42,23 +43,17 @@ def get_sample(cyl_height=5*nm):
 
 
 def get_simulation():
-    """
-    Returns a GISAXS simulation with beam and detector defined
-    """
-    simulation = ba.GISASSimulation()
-    simulation.setDetectorParameters(100, -2.0*deg, 2.0*deg, 100, 0.0*deg,
-                                     2.0*deg)
-    simulation.setBeamParameters(1.0*angstrom, 0.2*deg, 0.0*deg)
+    beam = ba.Beam(1.0, 0.1*nm, ba.Direction(0.2*deg, 0.0*deg))
+    detector = ba.SphericalDetector(100, -2.0*deg, 2.0*deg, 100, 0.0*deg,
+                                    2.0*deg)
+
+    simulation = ba.GISASSimulation(beam, get_sample(), detector)
     simulation.getOptions().setUseAvgMaterials(True)
     return simulation
 
 
 def run_simulation():
-    """
-    Runs simulation and returns intensity map.
-    """
     simulation = get_simulation()
-    simulation.setSample(get_sample())
     simulation.runSimulation()
     return simulation.result()
 
diff --git a/Examples/Python/sim04_Multilayers/HalfSpheresInAverageTopLayer.py b/Examples/Python/sim04_Multilayers/HalfSpheresInAverageTopLayer.py
index a691fbab52d14a78ea9c97b1a90129aade4781e1..eb9e7a9624d4d79e3c5c675fda5bfa09dd21365b 100644
--- a/Examples/Python/sim04_Multilayers/HalfSpheresInAverageTopLayer.py
+++ b/Examples/Python/sim04_Multilayers/HalfSpheresInAverageTopLayer.py
@@ -2,11 +2,9 @@
 Square lattice of half spheres on substrate with usage of average material
 and slicing
 """
+import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, nm, nm2, kvector_t
-
-sphere_radius = 5*nm
-n_slices = 10
+from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
 
 
 def get_sample():
@@ -55,24 +53,17 @@ def get_sample():
 
 
 def get_simulation():
-    """
-    Returns a GISAXS simulation with beam and detector defined
-    """
-    simulation = ba.GISASSimulation()
-    simulation.setDetectorParameters(100, -2.0*deg, 2.0*deg, 100, 0.0*deg,
-                                     2.0*deg)
-    simulation.setBeamParameters(1.0*angstrom, 0.2*deg, 0.0*deg)
+    beam = ba.Beam(1.0, 0.1*nm, ba.Direction(0.2*deg, 0.0*deg))
+    detector = ba.SphericalDetector(100, -2.0*deg, 2.0*deg, 100, 0.0*deg,
+                                    2.0*deg)
+
+    simulation = ba.GISASSimulation(beam, get_sample(), detector)
     simulation.getOptions().setUseAvgMaterials(True)
     return simulation
 
 
 def run_simulation():
-    """
-    Runs simulation and returns intensity map.
-    """
-    sample = get_sample()
     simulation = get_simulation()
-    simulation.setSample(sample)
     simulation.runSimulation()
     return simulation.result()
 
diff --git a/Examples/Python/sim05_Magnetism/MagneticSpheres.py b/Examples/Python/sim05_Magnetism/MagneticSpheres.py
index 171baed9ac7d00e256b8d063e3e2bed67e150fad..fc5b64a85c85ed06c553c2367ccfbd87de5f94af 100644
--- a/Examples/Python/sim05_Magnetism/MagneticSpheres.py
+++ b/Examples/Python/sim05_Magnetism/MagneticSpheres.py
@@ -1,11 +1,9 @@
 """
 Simulation demo: magnetic spheres in substrate
 """
+import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, nm, nm2, kvector_t
-
-# Magnetization of the particle's material (A/m)
-magnetization_particle = ba.kvector_t(0.0, 0.0, 1e7)
+from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
 
 
 def get_sample():
@@ -48,29 +46,20 @@ def get_sample():
 
 
 def get_simulation():
-    """
-    Returns a GISAXS simulation with beam and detector defined
-    """
-    simulation = ba.GISASSimulation()
-    simulation.setDetectorParameters(200, -3.0*deg, 3.0*deg, 200, 0.0*deg,
-                                     6.0*deg)
-    simulation.setBeamParameters(1.*angstrom, 0.5*deg, 0.0*deg)
-    simulation.beam().setIntensity(1e12)
-
-    analyzer_dir = ba.kvector_t(0.0, 0.0, -1.0)
-    beampol = ba.kvector_t(0.0, 0.0, 1.0)
-    simulation.setBeamPolarization(beampol)
-    simulation.setAnalyzerProperties(analyzer_dir, 1.0, 0.5)
-
+    beam = ba.Beam(1e+12, 0.1*nm, ba.Direction(0.5*deg, 0.0*deg))
+    beam_polarization = kvector_t(0.0, 0.0, 1.0)
+    beam.setPolarization(beam_polarization)
+    detector = ba.SphericalDetector(200, -3.0*deg, 3.0*deg, 200, 0.0*deg,
+                                    6.0*deg)
+
+    simulation = ba.GISASSimulation(beam, get_sample(), detector)
+    analyzer_direction = kvector_t(0.0, 0.0, -1.0)
+    simulation.setAnalyzerProperties(analyzer_direction, 1.0, 0.5)
     return simulation
 
 
 def run_simulation():
-    """
-    Runs simulation and returns intensity map.
-    """
     simulation = get_simulation()
-    simulation.setSample(get_sample())
     simulation.runSimulation()
     return simulation.result()
 
diff --git a/Examples/Python/sim11_Device/BeamDivergence.py b/Examples/Python/sim11_Device/BeamDivergence.py
index 9fdc00c06f875932142dd7e433719193725f30e0..dbcf9cd66a0bf18e382923bc810bb564977dbe70 100644
--- a/Examples/Python/sim11_Device/BeamDivergence.py
+++ b/Examples/Python/sim11_Device/BeamDivergence.py
@@ -1,8 +1,9 @@
 """
 Cylinder form factor in DWBA with beam divergence
 """
+import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, nm, nm2, kvector_t
+from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
 
 
 def get_sample():
@@ -41,32 +42,24 @@ def get_sample():
 
 
 def get_simulation():
-    """
-    Returns a GISAXS simulation with beam (+ divergence) and detector defined.
-    """
-    simulation = ba.GISASSimulation()
-    simulation.setDetectorParameters(100, 0.0*deg, 2.0*deg, 100, 0.0*deg,
-                                     2.0*deg)
-    simulation.setBeamParameters(1.0*angstrom, 0.2*deg, 0.0*deg)
-    wavelength_distr = ba.DistributionLogNormal(1.0*angstrom, 0.1)
-    alpha_distr = ba.DistributionGaussian(0.2*deg, 0.1*deg)
-    phi_distr = ba.DistributionGaussian(0.0*deg, 0.1*deg)
-    simulation.addParameterDistribution("*/Beam/Wavelength", wavelength_distr,
-                                        5)
-    simulation.addParameterDistribution("*/Beam/InclinationAngle", alpha_distr,
-                                        5)
-    simulation.addParameterDistribution("*/Beam/AzimuthalAngle", phi_distr, 5)
+    beam = ba.Beam(1.0, 0.1*nm, ba.Direction(0.2*deg, 0.0*deg))
+    detector = ba.SphericalDetector(100, 0.0*deg, 2.0*deg, 100, 0.0*deg,
+                                    2.0*deg)
+
+    simulation = ba.GISASSimulation(beam, get_sample(), detector)
+    distr_1 = ba.DistributionLogNormal(0.1*nm, 0.1)
+    simulation.addParameterDistribution("*/Beam/Wavelength", distr_1, 5, 0.0)
+    distr_2 = ba.DistributionGaussian(0.2*deg, 0.1*deg)
+    simulation.addParameterDistribution("*/Beam/InclinationAngle", distr_2, 5,
+                                        0.0)
+    distr_3 = ba.DistributionGaussian(0.0*deg, 0.1*deg)
+    simulation.addParameterDistribution("*/Beam/AzimuthalAngle", distr_3, 5,
+                                        0.0)
     return simulation
 
 
 def run_simulation():
-    """
-    Runs simulation and returns intensity map.
-    """
     simulation = get_simulation()
-    simulation.setSample(get_sample())
-    print(simulation.treeToString())
-    print(simulation.parametersToString())
     simulation.runSimulation()
     return simulation.result()
 
diff --git a/Examples/Python/sim11_Device/ConstantBackground.py b/Examples/Python/sim11_Device/ConstantBackground.py
index a1e5e980907b73953ac9a142f39ef471663e25ad..b3bb577106ea15d8d60e9359571dd1a737c2ea13 100644
--- a/Examples/Python/sim11_Device/ConstantBackground.py
+++ b/Examples/Python/sim11_Device/ConstantBackground.py
@@ -1,8 +1,9 @@
 """
 Cylinder form factor in DWBA with constant background
 """
+import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, nm, nm2, kvector_t
+from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
 
 
 def get_sample():
@@ -41,26 +42,18 @@ def get_sample():
 
 
 def get_simulation():
-    """
-    Returns a GISAXS simulation with a constant backround.
-    """
-    simulation = ba.GISASSimulation()
-    simulation.setDetectorParameters(100, 0.0*deg, 2.0*deg, 100, 0.0*deg,
-                                     2.0*deg)
-    simulation.setBeamParameters(1.0*angstrom, 0.2*deg, 0.0*deg)
-    simulation.beam().setIntensity(1e6)
-    bg = ba.ConstantBackground(1e3)
-    simulation.setBackground(bg)
+    beam = ba.Beam(1000000.0, 0.1*nm, ba.Direction(0.2*deg, 0.0*deg))
+    detector = ba.SphericalDetector(100, 0.0*deg, 2.0*deg, 100, 0.0*deg,
+                                    2.0*deg)
+
+    simulation = ba.GISASSimulation(beam, get_sample(), detector)
+    background = ba.ConstantBackground(1.0e+03)
+    simulation.setBackground(background)
     return simulation
 
 
 def run_simulation():
-    """
-    Runs simulation and returns intensity map.
-    """
     simulation = get_simulation()
-    simulation.setSample(get_sample())
-    print(simulation.parametersToString())
     simulation.runSimulation()
     return simulation.result()
 
diff --git a/Examples/Python/sim11_Device/DetectorResolutionFunction.py b/Examples/Python/sim11_Device/DetectorResolutionFunction.py
index 5add835f7998d4aad0f4c544aee6734061adfb37..1142a1e9c8a8dd11837218477e07b1da5a8967a8 100644
--- a/Examples/Python/sim11_Device/DetectorResolutionFunction.py
+++ b/Examples/Python/sim11_Device/DetectorResolutionFunction.py
@@ -1,8 +1,9 @@
 """
 Cylinder form factor in DWBA with detector resolution function applied
 """
+import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, nm, nm2, kvector_t
+from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
 
 
 def get_sample():
@@ -41,24 +42,18 @@ def get_sample():
 
 
 def get_simulation():
-    """
-    Returns a GISAXS simulation with detector resolution function defined.
-    """
-    simulation = ba.GISASSimulation()
-    simulation.setDetectorParameters(100, 0.0*deg, 2.0*deg, 100, 0.0*deg,
-                                     2.0*deg)
-    simulation.setBeamParameters(1.0*angstrom, 0.2*deg, 0.0*deg)
+    beam = ba.Beam(1.0, 0.1*nm, ba.Direction(0.2*deg, 0.0*deg))
+    detector = ba.SphericalDetector(100, 0.0*deg, 2.0*deg, 100, 0.0*deg,
+                                    2.0*deg)
+
+    simulation = ba.GISASSimulation(beam, get_sample(), detector)
     simulation.setDetectorResolutionFunction(
         ba.ResolutionFunction2DGaussian(0.02*deg, 0.02*deg))
     return simulation
 
 
 def run_simulation():
-    """
-    Runs simulation and returns intensity map.
-    """
     simulation = get_simulation()
-    simulation.setSample(get_sample())
     simulation.runSimulation()
     return simulation.result()
 
diff --git a/Examples/Python/sim22_OffSpecular/BoxesWithSpecularPeak.py b/Examples/Python/sim22_OffSpecular/BoxesWithSpecularPeak.py
index e5d58ff02734e888aacc75f51f14792d2d29278d..ca4bda70e1cec9492d6ad1dc5956fc64665f76ea 100644
--- a/Examples/Python/sim22_OffSpecular/BoxesWithSpecularPeak.py
+++ b/Examples/Python/sim22_OffSpecular/BoxesWithSpecularPeak.py
@@ -1,8 +1,9 @@
 """
 Square lattice of boxes on substrate (including the specular peak)
 """
+import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, nm, nm2, kvector_t
+from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
 
 
 def get_sample():
@@ -50,25 +51,18 @@ def get_sample():
 
 
 def get_simulation():
-    """
-    Returns a GISAXS simulation with beam and detector defined
-    """
-    simulation = ba.GISASSimulation()
-    simulation.setDetectorParameters(101, -2.0*deg, 2.0*deg, 101, 0.0*deg,
-                                     2.0*deg)
-    simulation.setBeamParameters(1.0*angstrom, 0.2*deg, 0.0*deg)
+    beam = ba.Beam(1.0, 0.1*nm, ba.Direction(0.2*deg, 0.0*deg))
+    detector = ba.SphericalDetector(101, -2.0*deg, 2.0*deg, 101, 0.0*deg,
+                                    2.0*deg)
+
+    simulation = ba.GISASSimulation(beam, get_sample(), detector)
     simulation.getOptions().setUseAvgMaterials(True)
     simulation.getOptions().setIncludeSpecular(True)
     return simulation
 
 
 def run_simulation():
-    """
-    Runs simulation and returns intensity map.
-    """
-    sample = get_sample()
     simulation = get_simulation()
-    simulation.setSample(sample)
     simulation.runSimulation()
     return simulation.result()