diff --git a/CMakeLists.txt b/CMakeLists.txt
index 7903e34465e01ae0f9d21cd9b2479bf0d436f7bc..6624dfe3a08621938f65a9f644068142c6ccbf6b 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -156,6 +156,7 @@ endif()
 add_subdirectory(Tests/UnitTests/Core)
 add_subdirectory(Tests/UnitTests/Numeric)
 add_subdirectory(Tests/Performance/Core)
+add_subdirectory(Tests/Examples)
 
 # GUI
 if(BORNAGAIN_GUI)
diff --git a/Core/Export/PyFmt.cpp b/Core/Export/PyFmt.cpp
index 410a88b08cc209c1e8a480a9f2f161ed548a1dc2..4a8c9faa7459a9c60a7828af9f282b9829ec3289 100644
--- a/Core/Export/PyFmt.cpp
+++ b/Core/Export/PyFmt.cpp
@@ -16,14 +16,23 @@
 #include "Base/Const/Units.h" // printDegrees
 #include "Base/Math/Constants.h"
 #include "Base/Utils/Algorithms.h"
+#include "Base/Utils/StringUtils.h"
 #include <iomanip>
 
 namespace pyfmt {
 
-std::string scriptPreamble() {
+std::string preambled(const std::string& code) {
+    std::vector<std::string> to_declare;
+    for (const std::string& key : {"angstrom", "deg", "nm", "nm2", "micrometer"})
+        if (code.find("*" + key) != std::string::npos)
+            to_declare.push_back(key);
+    for (const std::string& key : {"kvector_t"})
+        if (code.find(key) != std::string::npos)
+            to_declare.push_back(key);
     return "import numpy, sys\n"
            "import bornagain as ba\n"
-           "from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t\n\n\n";
+           "from bornagain import "
+           + StringUtils::join(to_declare, ", ") + "\n\n\n" + code;
 }
 
 std::string printBool(double value) {
diff --git a/Core/Export/PyFmt.h b/Core/Export/PyFmt.h
index d123014c5aaf557bdf3f75609fa4a470184f8f04..1b45fcdeb43eb50f43d69ee87c388222982d37a2 100644
--- a/Core/Export/PyFmt.h
+++ b/Core/Export/PyFmt.h
@@ -27,7 +27,7 @@
 
 namespace pyfmt {
 
-std::string scriptPreamble();
+std::string preambled(const std::string& code);
 
 std::string printInt(int value);
 std::string printBool(double value);
diff --git a/Core/Export/SimulationToPython.cpp b/Core/Export/SimulationToPython.cpp
index 70fe7e8dfa3f60ea94ef2d24c1541168fad00c6d..25494c0a27039f206f5e039cc916b93cb80d2dca 100644
--- a/Core/Export/SimulationToPython.cpp
+++ b/Core/Export/SimulationToPython.cpp
@@ -438,8 +438,9 @@ std::string defineSimulate(const ISimulation* simulation) {
 std::string simulationCode(const ISimulation& simulation) {
     if (simulation.sample() == nullptr)
         throw std::runtime_error("Cannot export: Simulation has no sample");
-    return pyfmt::scriptPreamble() + SampleToPython().sampleCode(*simulation.sample())
-           + defineSimulate(&simulation);
+    std::string code =
+        SampleToPython().sampleCode(*simulation.sample()) + defineSimulate(&simulation);
+    return pyfmt::preambled(code);
 }
 
 } // namespace
diff --git a/Examples/scatter2d/ApproximationDA.py b/Examples/scatter2d/ApproximationDA.py
index fb08f9e29b6f2f5c3e4ce267b25d846950140e7d..f6d60aa1180586eeeab8b939f12b8ceb6d32bb24 100644
--- a/Examples/scatter2d/ApproximationDA.py
+++ b/Examples/scatter2d/ApproximationDA.py
@@ -3,7 +3,7 @@ Cylinders of two different sizes in Decoupling Approximation
 """
 import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
+from bornagain import deg, nm
 
 
 def get_sample():
diff --git a/Examples/scatter2d/ApproximationLMA.py b/Examples/scatter2d/ApproximationLMA.py
index 6b0b30b9b039f88cdfa7f07e42220c05821ce2ff..1d69886e5bec82e76e1fd209f4b692dbd72abe92 100644
--- a/Examples/scatter2d/ApproximationLMA.py
+++ b/Examples/scatter2d/ApproximationLMA.py
@@ -3,7 +3,7 @@ Cylinders of two different sizes in Local Monodisperse Approximation
 """
 import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
+from bornagain import deg, nm
 
 
 def get_sample():
diff --git a/Examples/scatter2d/ApproximationSSCA.py b/Examples/scatter2d/ApproximationSSCA.py
index c2cfe56158b6b49afa396b0468fc74288280e41d..287dc47c68c67264dc589b8afcd11b9d6c380857 100644
--- a/Examples/scatter2d/ApproximationSSCA.py
+++ b/Examples/scatter2d/ApproximationSSCA.py
@@ -3,7 +3,7 @@ Cylinders of two different sizes in Size-Spacing Coupling Approximation
 """
 import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
+from bornagain import deg, nm
 
 
 def get_sample():
diff --git a/Examples/scatter2d/BeamDivergence.py b/Examples/scatter2d/BeamDivergence.py
index 302ae4f0cc2773f95b68c5412f1c59e6a31daa3d..e8d4a04fc849a9c83b6db60e7a0d6dae763f73f2 100644
--- a/Examples/scatter2d/BeamDivergence.py
+++ b/Examples/scatter2d/BeamDivergence.py
@@ -3,7 +3,7 @@ Cylinder form factor in DWBA with beam divergence
 """
 import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
+from bornagain import deg, nm
 
 
 def get_sample():
diff --git a/Examples/scatter2d/BiMaterialCylinders.py b/Examples/scatter2d/BiMaterialCylinders.py
index 978f0d4ab630d413ba62c29f11d41715bf8abd7b..ac991ce7e7984f89c7c0947be71ff589998533cd 100644
--- a/Examples/scatter2d/BiMaterialCylinders.py
+++ b/Examples/scatter2d/BiMaterialCylinders.py
@@ -4,7 +4,7 @@ Particle crosses air/substrate interface.
 """
 import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
+from bornagain import deg, nm, kvector_t
 
 
 def get_sample():
diff --git a/Examples/scatter2d/BoxesWithSpecularPeak.py b/Examples/scatter2d/BoxesWithSpecularPeak.py
index 7ee7bb4ffe0b89ec7bdcc6d5b12dfa5a48e0316f..ca87675492c734fa9eba542365b52daf66401d9f 100644
--- a/Examples/scatter2d/BoxesWithSpecularPeak.py
+++ b/Examples/scatter2d/BoxesWithSpecularPeak.py
@@ -3,7 +3,7 @@ Square lattice of boxes on substrate (including the specular peak)
 """
 import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
+from bornagain import deg, nm
 
 
 def get_sample():
diff --git a/Examples/scatter2d/BuriedParticles.py b/Examples/scatter2d/BuriedParticles.py
index 34b4180d02973812f183dcd138dbf3f18054eb38..6fb9cd0efa19b48c9610a586817e22e4e5bdec8d 100644
--- a/Examples/scatter2d/BuriedParticles.py
+++ b/Examples/scatter2d/BuriedParticles.py
@@ -3,7 +3,7 @@ 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, micrometer, nm, nm2, kvector_t
+from bornagain import deg, nm, kvector_t
 
 
 def get_sample():
diff --git a/Examples/scatter2d/ConstantBackground.py b/Examples/scatter2d/ConstantBackground.py
index a7e03143688a85f15ea437f78f7aab9c59d45fe6..5244d0ceb44ecb4422243b28553e06a4547a048e 100644
--- a/Examples/scatter2d/ConstantBackground.py
+++ b/Examples/scatter2d/ConstantBackground.py
@@ -3,7 +3,7 @@ Cylinder form factor in DWBA with constant background
 """
 import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
+from bornagain import deg, nm
 
 
 def get_sample():
diff --git a/Examples/scatter2d/CoreShellNanoparticles.py b/Examples/scatter2d/CoreShellNanoparticles.py
index f3b89ce66ef85fd8cdeb2d50084f3ae7a775f774..19162de038f9f1adab3e4f0dc3cf408886ad6a95 100644
--- a/Examples/scatter2d/CoreShellNanoparticles.py
+++ b/Examples/scatter2d/CoreShellNanoparticles.py
@@ -3,7 +3,7 @@ Core shell nanoparticles
 """
 import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
+from bornagain import deg, nm
 
 
 def get_sample():
diff --git a/Examples/scatter2d/CorrelatedRoughness.py b/Examples/scatter2d/CorrelatedRoughness.py
index de047c5913e137db3d70f47ceb4619dccc054a7a..916215b5f3d6495c0c40acd3e3af733d49bd805c 100644
--- a/Examples/scatter2d/CorrelatedRoughness.py
+++ b/Examples/scatter2d/CorrelatedRoughness.py
@@ -3,7 +3,7 @@ MultiLayer with correlated roughness
 """
 import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
+from bornagain import deg, nm
 
 
 def get_sample():
diff --git a/Examples/scatter2d/CosineRipplesAtRectLattice.py b/Examples/scatter2d/CosineRipplesAtRectLattice.py
index 253aca9bd09f814e0063d7c44a63e510d6277fd8..694f055b18d42239bd4c27d9755f64ea177dffd2 100644
--- a/Examples/scatter2d/CosineRipplesAtRectLattice.py
+++ b/Examples/scatter2d/CosineRipplesAtRectLattice.py
@@ -3,7 +3,7 @@ Cosine ripple on a 2D lattice
 """
 import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
+from bornagain import deg, nm
 
 
 def get_sample():
diff --git a/Examples/scatter2d/CylindersAndPrisms.py b/Examples/scatter2d/CylindersAndPrisms.py
index 6fd88b86c0f415a8016c873912fdc3b2c70b9a93..216d9afa21fc7e527523825747332f02a575f617 100644
--- a/Examples/scatter2d/CylindersAndPrisms.py
+++ b/Examples/scatter2d/CylindersAndPrisms.py
@@ -3,7 +3,7 @@ Mixture of cylinders and prisms without interference
 """
 import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
+from bornagain import deg, nm
 
 
 def get_sample():
diff --git a/Examples/scatter2d/CylindersInAverageLayer.py b/Examples/scatter2d/CylindersInAverageLayer.py
index 8f1b9206cd3a04479e4a311f6ab1e53e744d5ac9..e13db805276b53b4f7ef36268695407b952b2719 100644
--- a/Examples/scatter2d/CylindersInAverageLayer.py
+++ b/Examples/scatter2d/CylindersInAverageLayer.py
@@ -3,7 +3,7 @@ Square lattice of cylinders inside finite layer with usage of average material
 """
 import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
+from bornagain import deg, nm, kvector_t
 
 
 def get_sample(cyl_height=5*nm):
diff --git a/Examples/scatter2d/CylindersInBA.py b/Examples/scatter2d/CylindersInBA.py
index b458b02c73e3ba0e818301041151bd45628123f2..ec8d50a4d19c6726c6712150220df56ddfc9f6f4 100644
--- a/Examples/scatter2d/CylindersInBA.py
+++ b/Examples/scatter2d/CylindersInBA.py
@@ -3,7 +3,7 @@ Cylinder form factor in Born approximation
 """
 import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
+from bornagain import deg, nm
 
 
 def get_sample():
diff --git a/Examples/scatter2d/CylindersInDWBA.py b/Examples/scatter2d/CylindersInDWBA.py
index 0627bc21b0a9ff0d685733a7c5f575f364c063c8..076dd5960ea6f8c30f582d65bec6fcfe29418433 100644
--- a/Examples/scatter2d/CylindersInDWBA.py
+++ b/Examples/scatter2d/CylindersInDWBA.py
@@ -3,7 +3,7 @@ Cylinder form factor in DWBA
 """
 import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
+from bornagain import deg, nm
 
 
 def get_sample():
diff --git a/Examples/scatter2d/CylindersWithSizeDistribution.py b/Examples/scatter2d/CylindersWithSizeDistribution.py
index c81b9ee26e1f98ad9573b00dd2ec7151b6894696..65439c35962bcc4d4dc761d58ec248073f05f2ec 100644
--- a/Examples/scatter2d/CylindersWithSizeDistribution.py
+++ b/Examples/scatter2d/CylindersWithSizeDistribution.py
@@ -3,7 +3,7 @@ Cylinders with size distribution
 """
 import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
+from bornagain import deg, nm
 
 
 def get_sample():
diff --git a/Examples/scatter2d/DetectorResolutionFunction.py b/Examples/scatter2d/DetectorResolutionFunction.py
index 29850a777e2078dc0d93489d2156202086cce1a6..d6bde346d70437cfb77b3464cce1d7995823d38c 100644
--- a/Examples/scatter2d/DetectorResolutionFunction.py
+++ b/Examples/scatter2d/DetectorResolutionFunction.py
@@ -3,7 +3,7 @@ Cylinder form factor in DWBA with detector resolution function applied
 """
 import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
+from bornagain import deg, nm
 
 
 def get_sample():
diff --git a/Examples/scatter2d/HalfSpheresInAverageTopLayer.py b/Examples/scatter2d/HalfSpheresInAverageTopLayer.py
index f03fb8bd2e330c195f5e7e69954857f9d9110119..d8933487c15ee857e3fe4d076279a20aa13c5d1e 100644
--- a/Examples/scatter2d/HalfSpheresInAverageTopLayer.py
+++ b/Examples/scatter2d/HalfSpheresInAverageTopLayer.py
@@ -4,7 +4,7 @@ and slicing
 """
 import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
+from bornagain import deg, nm
 
 
 def get_sample():
diff --git a/Examples/scatter2d/HexagonalLatticesWithBasis.py b/Examples/scatter2d/HexagonalLatticesWithBasis.py
index 032dc1c7682685f3945d3a5f704ac617f4bc9f0e..dd99d3f2c12d94698c014dbb3744acd052252230 100644
--- a/Examples/scatter2d/HexagonalLatticesWithBasis.py
+++ b/Examples/scatter2d/HexagonalLatticesWithBasis.py
@@ -3,7 +3,7 @@ Spheres on two hexagonal close packed layers
 """
 import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
+from bornagain import deg, nm, kvector_t
 
 
 def get_sample():
diff --git a/Examples/scatter2d/Interference1DRadialParaCrystal.py b/Examples/scatter2d/Interference1DRadialParaCrystal.py
index 487b49d4766732a0a143016a4187345679c5bd94..b17da17637cf6aeaee65c27d0425fb0ad0e890be 100644
--- a/Examples/scatter2d/Interference1DRadialParaCrystal.py
+++ b/Examples/scatter2d/Interference1DRadialParaCrystal.py
@@ -3,7 +3,7 @@ radial paracrystal
 """
 import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
+from bornagain import deg, nm
 
 
 def get_sample():
diff --git a/Examples/scatter2d/Interference2DCenteredSquareLattice.py b/Examples/scatter2d/Interference2DCenteredSquareLattice.py
index 23e5f2d262d6b5c08935128bebfae9567a48fae6..d56deadc3ca34f85a50c50437d04b6f0d6af1517 100644
--- a/Examples/scatter2d/Interference2DCenteredSquareLattice.py
+++ b/Examples/scatter2d/Interference2DCenteredSquareLattice.py
@@ -3,7 +3,7 @@
 """
 import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
+from bornagain import deg, nm, kvector_t
 
 
 def get_sample():
diff --git a/Examples/scatter2d/Interference2DLatticeSumOfRotated.py b/Examples/scatter2d/Interference2DLatticeSumOfRotated.py
index ffd4fd6daa91e3d9c781155d1234d2d13f35372e..ac2a432a9da7e192364851d3d6305380eaaab203 100644
--- a/Examples/scatter2d/Interference2DLatticeSumOfRotated.py
+++ b/Examples/scatter2d/Interference2DLatticeSumOfRotated.py
@@ -1,6 +1,6 @@
 import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
+from bornagain import deg, nm
 
 
 def get_sample():
diff --git a/Examples/scatter2d/Interference2DParaCrystal.py b/Examples/scatter2d/Interference2DParaCrystal.py
index e29de7d9a65f9a8411117348b22378c6b3248180..013a65d4b458132bea58253b2009e6eea19ecc72 100644
--- a/Examples/scatter2d/Interference2DParaCrystal.py
+++ b/Examples/scatter2d/Interference2DParaCrystal.py
@@ -3,7 +3,7 @@
 """
 import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
+from bornagain import deg, nm
 
 
 def get_sample():
diff --git a/Examples/scatter2d/Interference2DRotatedSquareLattice.py b/Examples/scatter2d/Interference2DRotatedSquareLattice.py
index ac18356a1f3dc4712292cb12e578b5cd77a35078..f428b4409637b43848f0d11f24738b227de4494d 100644
--- a/Examples/scatter2d/Interference2DRotatedSquareLattice.py
+++ b/Examples/scatter2d/Interference2DRotatedSquareLattice.py
@@ -3,7 +3,7 @@ Cylinders on a rotated 2D lattice
 """
 import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
+from bornagain import deg, nm
 
 
 def get_sample():
diff --git a/Examples/scatter2d/Interference2DSquareFiniteLattice.py b/Examples/scatter2d/Interference2DSquareFiniteLattice.py
index d028bbb99937973108e27caf531b79dfe083ce7d..5fc42c81f4853502b19feb56df8f809956b47f50 100644
--- a/Examples/scatter2d/Interference2DSquareFiniteLattice.py
+++ b/Examples/scatter2d/Interference2DSquareFiniteLattice.py
@@ -3,7 +3,7 @@ Cylinders on a 2D square lattice
 """
 import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
+from bornagain import deg, nm, nm2
 
 
 def get_sample():
diff --git a/Examples/scatter2d/Interference2DSquareLattice.py b/Examples/scatter2d/Interference2DSquareLattice.py
index 1997f7acdc49eb2287d1b358849058484b607c9b..cb64a4e5415e4e78e4b21bd2883147c6f3a03859 100644
--- a/Examples/scatter2d/Interference2DSquareLattice.py
+++ b/Examples/scatter2d/Interference2DSquareLattice.py
@@ -3,7 +3,7 @@ Cylinders on a 2D square lattice
 """
 import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
+from bornagain import deg, nm
 
 
 def get_sample():
diff --git a/Examples/scatter2d/MagneticSpheres.py b/Examples/scatter2d/MagneticSpheres.py
index 6b688d8a99e1b37d573efeb30a080be9370a3ef2..ef10a3bf24d124e387b2659f3c98bea79e3937e5 100644
--- a/Examples/scatter2d/MagneticSpheres.py
+++ b/Examples/scatter2d/MagneticSpheres.py
@@ -3,7 +3,7 @@ Simulation demo: magnetic spheres in substrate
 """
 import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
+from bornagain import deg, nm, kvector_t
 
 
 def get_sample():
diff --git a/Examples/scatter2d/MesoCrystal.py b/Examples/scatter2d/MesoCrystal.py
index f2dc9d1ce244b7706b7c12af949eec738a44378d..c238716ca40d635712e945e67c0a2f217c124cc6 100644
--- a/Examples/scatter2d/MesoCrystal.py
+++ b/Examples/scatter2d/MesoCrystal.py
@@ -3,7 +3,7 @@ Cylindrical mesocrystal on a substrate
 """
 import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
+from bornagain import deg, nm, kvector_t
 
 
 def get_sample():
diff --git a/Examples/scatter2d/ParticlesCrossingInterface.py b/Examples/scatter2d/ParticlesCrossingInterface.py
index 0aeebe7b83f8523ff88d103800b03ee8e18cb5d0..a3be3b18f528851b620ad51d834f92e1a768bf2e 100644
--- a/Examples/scatter2d/ParticlesCrossingInterface.py
+++ b/Examples/scatter2d/ParticlesCrossingInterface.py
@@ -18,7 +18,7 @@ will be thrown when trying to simulate such geometries).
 """
 import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
+from bornagain import deg, nm, kvector_t
 
 
 def get_sample():
diff --git a/Examples/scatter2d/RotatedPyramids.py b/Examples/scatter2d/RotatedPyramids.py
index b56a7506accd9f57c003482e838b6077715afabd..a4e73820ca7c2b52e4464c1112b5aaf77f4ca06c 100644
--- a/Examples/scatter2d/RotatedPyramids.py
+++ b/Examples/scatter2d/RotatedPyramids.py
@@ -3,7 +3,7 @@ Rotated pyramids on top of substrate
 """
 import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
+from bornagain import deg, nm
 
 
 def get_sample():
diff --git a/Examples/scatter2d/SpheresAtHexLattice.py b/Examples/scatter2d/SpheresAtHexLattice.py
index 6397071faa399b9e76de27f2604b49353c1e65ae..961213c595bfb5553f4b08f3e02a0a7378265d44 100644
--- a/Examples/scatter2d/SpheresAtHexLattice.py
+++ b/Examples/scatter2d/SpheresAtHexLattice.py
@@ -3,7 +3,7 @@ Spheres on a hexagonal lattice
 """
 import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
+from bornagain import deg, nm
 
 
 def get_sample():
diff --git a/Examples/scatter2d/TriangularRipple.py b/Examples/scatter2d/TriangularRipple.py
index c1dc7cb3777d41a3e418628e9c050f7d1f59b671..628e36846f17f7cb7155402f7e996c4cf1556311 100644
--- a/Examples/scatter2d/TriangularRipple.py
+++ b/Examples/scatter2d/TriangularRipple.py
@@ -3,7 +3,7 @@
 """
 import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
+from bornagain import deg, nm
 
 
 def get_sample():
diff --git a/Examples/scatter2d/TwoTypesOfCylindersWithSizeDistribution.py b/Examples/scatter2d/TwoTypesOfCylindersWithSizeDistribution.py
index 92cd20b45642f466b478ff4c03d1adc90a5b2c2d..7f4450d0fd4f96d33513ae5479dcc635b1bd1db1 100644
--- a/Examples/scatter2d/TwoTypesOfCylindersWithSizeDistribution.py
+++ b/Examples/scatter2d/TwoTypesOfCylindersWithSizeDistribution.py
@@ -3,7 +3,7 @@ Mixture cylinder particles with different size distribution
 """
 import numpy, sys
 import bornagain as ba
-from bornagain import angstrom, deg, micrometer, nm, nm2, kvector_t
+from bornagain import deg, nm
 
 
 def get_sample():
diff --git a/Tests/Functional/Python/PyPersistence/CMakeLists.txt b/Tests/Examples/CMakeLists.txt
similarity index 60%
rename from Tests/Functional/Python/PyPersistence/CMakeLists.txt
rename to Tests/Examples/CMakeLists.txt
index 2773f9a39f797c29f70831ed5f95563ba6a069b7..13e979ccb6228508398211f96c4575e37267c9fe 100644
--- a/Tests/Functional/Python/PyPersistence/CMakeLists.txt
+++ b/Tests/Examples/CMakeLists.txt
@@ -1,5 +1,54 @@
+###############################################################################
+# Run unmodified examples
+###############################################################################
+
+set(test_script ${TOOL_DIR}/code/batch-plot.py)
+
+set(output_dir ${TEST_OUTPUT_DIR_PY_EXAMPLES})
+file(MAKE_DIRECTORY ${output_dir})
+
+# Run one Python example.
+# Check whether the example passes.
+# Don't check results.
+# If there are matplotlib commands, skip plt.show().
+function(run_example example label)
+    get_filename_component(name ${example} NAME_WE)
+    set(test_name Example.run.${name})
+    add_test(${test_name}
+        env PYTHONPATH=${CMAKE_LIBRARY_OUTPUT_DIRECTORY} NOSHOW=TRUE
+        ${Python3_EXECUTABLE} ${test_script} -s ${example} ${output_dir})
+    if(NOT label STREQUAL "")
+        set_tests_properties(${test_name} PROPERTIES LABELS ${label})
+    endif()
+endfunction()
+
+function(run_examples examples label)
+    foreach(example ${examples})
+        run_example(${example} ${label})
+    endforeach()
+endfunction()
+
+# Examples from scatter2d; all but one prototype get the label XL
+set(prototype ${EXAMPLES_DIR}/scatter2d/CylindersInDWBA.py)
+file(GLOB examples ${EXAMPLES_DIR}/scatter2d/*.py)
+list(REMOVE_ITEM examples ${prototype})
+run_example(${prototype} "")
+run_examples(${examples} "XL")
+
+# Examples from other directories
+file(GLOB examples
+    ${EXAMPLES_DIR}/specular/*.py
+    ${EXAMPLES_DIR}/varia/*.py
+    ${EXAMPLES_DIR}/fit55_Specular/FitSpecularBasics.py)
+if (WIN32)
+    # Convergence problem in Gauss-Kronrod integration
+    list(REMOVE_ITEM examples ${EXAMPLES_DIR}/scatter2d/Interference2DParaCrystal.py)
+endif()
+list(REMOVE_ITEM examples ${EXAMPLES_DIR}/varia/SimulationParameters.py)
+run_examples(${examples} "")
+
 ############################################################################
-# Tests/Functional/PyCore/persistence/CMakeLists.txt
+# Check results of modified examples
 ############################################################################
 
 set(OUTPUT_DIR ${TEST_OUTPUT_DIR_PY_PERSIST})
@@ -11,7 +60,7 @@ function(test_example example tolerance)
     get_filename_component(EXAMPLE_NAME ${script_path} NAME_WE)
     get_filename_component(EXAMPLE_DIR ${script_path} DIRECTORY)
 
-    set(test_name PyPersistence.${EXAMPLE_NAME})
+    set(test_name Example.persist.${EXAMPLE_NAME})
 
     set(PYPERSIST_TOLERANCE ${tolerance})
 
diff --git a/Tests/Functional/Python/PyPersistence/PyPersistence.py.in b/Tests/Examples/PyPersistence.py.in
similarity index 100%
rename from Tests/Functional/Python/PyPersistence/PyPersistence.py.in
rename to Tests/Examples/PyPersistence.py.in
diff --git a/Tests/Functional/Python/CMakeLists.txt b/Tests/Functional/Python/CMakeLists.txt
index 2382fd012d3151a06a2b9a02a27a4b4950c27d79..09415828a2ee0e129f4bed52b0018804c2117f87 100644
--- a/Tests/Functional/Python/CMakeLists.txt
+++ b/Tests/Functional/Python/CMakeLists.txt
@@ -1,8 +1,6 @@
 # Collection of Python related tests
 
 add_subdirectory(Std)
-add_subdirectory(PyExamples)
 add_subdirectory(PyCore)
 add_subdirectory(PyFit)
-add_subdirectory(PyPersistence)
 add_subdirectory(PyEmbedded)
diff --git a/Tests/Functional/Python/PyEmbedded/Tests.cpp b/Tests/Functional/Python/PyEmbedded/Tests.cpp
index 09b944e020ff3d9ee0d08238d4cfce861ab2f472..1ee49ed20661749911b49ff1490200dd805fd302 100644
--- a/Tests/Functional/Python/PyEmbedded/Tests.cpp
+++ b/Tests/Functional/Python/PyEmbedded/Tests.cpp
@@ -362,14 +362,12 @@ TEST_F(PyEmbedded, ExportToPythonAndBack) {
     SampleBuilderFactory factory;
     std::unique_ptr<MultiLayer> sample(factory.createSampleByName("CylindersAndPrismsBuilder"));
 
-    auto code = ExportToPython::sampleCode(*sample);
+    const std::string code = ExportToPython::sampleCode(*sample);
+    const std::string snippet = pyfmt::preambled(code);
 
-    std::stringstream snippet;
-    snippet << pyfmt::scriptPreamble() << code;
-
-    auto multilayer =
-        PyImport::createFromPython(snippet.str(), "get_sample", BABuild::buildLibDir());
-    auto new_code = ExportToPython::sampleCode(*multilayer);
+    const auto multilayer =
+        PyImport::createFromPython(snippet, "get_sample", BABuild::buildLibDir());
+    const std::string new_code = ExportToPython::sampleCode(*multilayer);
 
     EXPECT_TRUE(code == new_code);
 }
diff --git a/Tests/Functional/Python/PyExamples/CMakeLists.txt b/Tests/Functional/Python/PyExamples/CMakeLists.txt
deleted file mode 100644
index 74149d49dde4584bc52b8a504a0b876797857498..0000000000000000000000000000000000000000
--- a/Tests/Functional/Python/PyExamples/CMakeLists.txt
+++ /dev/null
@@ -1,30 +0,0 @@
-###############################################################################
-# Tests/Functional/PyExamples/CMakeLists.txt
-#
-# Runs check_functionality.py over example scripts found in PY_EXAMPLES_DIR
-#
-###############################################################################
-
-set(output_dir ${TEST_OUTPUT_DIR_PY_EXAMPLES})
-file(MAKE_DIRECTORY ${output_dir})
-
-file(GLOB examples
-    ${EXAMPLES_DIR}/scatter2d/*.py
-    ${EXAMPLES_DIR}/specular/*.py
-    ${EXAMPLES_DIR}/varia/*.py
-    ${EXAMPLES_DIR}/fit55_Specular/FitSpecularBasics.py)
-if (WIN32)
-    # Convergence problem in Gauss-Kronrod integration
-    list(REMOVE_ITEM examples ${EXAMPLES_DIR}/scatter2d/Interference2DParaCrystal.py)
-endif()
-list(REMOVE_ITEM examples ${EXAMPLES_DIR}/varia/SimulationParameters.py)
-
-set(test_script ${TOOL_DIR}/code/batch-plot.py)
-
-foreach(example_path ${examples})
-    get_filename_component(example_name ${example_path} NAME_WE)
-    set(test_name PyExamples.${example_name})
-    add_test(${test_name}
-        env PYTHONPATH=${CMAKE_LIBRARY_OUTPUT_DIRECTORY}
-        ${Python3_EXECUTABLE} ${test_script} -s ${example_path} ${output_dir})
-endforeach()
diff --git a/Wrap/Python/plot_utils.py b/Wrap/Python/plot_utils.py
index 2cf3c73e9f0d2f0e1cbf2edf1f056a9e104876cf..0f770b92b2554f19df20d4e62513379ec9f9b0f8 100644
--- a/Wrap/Python/plot_utils.py
+++ b/Wrap/Python/plot_utils.py
@@ -23,6 +23,8 @@ except Exception as e:
 
 label_fontsize = 16
 
+CMAP = "CMAP" in os.environ and os.environ["CMAP"] or 'inferno'
+
 
 def get_axes_limits(result, units):
     """
@@ -100,10 +102,10 @@ def plot_array(array, axes_limits=None, **kwargs):
     xlabel = kwargs.pop('xlabel', None)
     ylabel = kwargs.pop('ylabel', None)
     zlabel = kwargs.pop('zlabel', "Intensity")
-
     title = kwargs.pop('title', None)
+    cmap = kwargs.pop('cmap', CMAP)
 
-    im = plt.imshow(array, norm=norm, extent=axes_limits, **kwargs)
+    im = plt.imshow(array, cmap=cmap, norm=norm, extent=axes_limits, **kwargs)
     cb = plt.colorbar(im, pad=0.025)
 
     if xlabel:
diff --git a/auto/Wrap/doxygenCore.i b/auto/Wrap/doxygenCore.i
index afa3f4f9c1fb51d8e4c1e3959c66f1cc8ad2d1f2..479b47904f7e200e58aa704a47d114b82870fcbf 100644
--- a/auto/Wrap/doxygenCore.i
+++ b/auto/Wrap/doxygenCore.i
@@ -1723,6 +1723,12 @@ Sets qz resolution values via IRangedDistribution and values of relative deviati
 Sets qz resolution values via IRangedDistribution and values of standard deviations.  std_dev can be either single-valued or a numpy array. In the latter case the length of the array should coinside with the length of the qz-axis. 
 ";
 
+%feature("docstring")  QSpecScan::setOffset "void QSpecScan::setOffset(double offset)
+";
+
+%feature("docstring")  QSpecScan::offset "double QSpecScan::offset() const
+";
+
 
 // File: classRelativeDifferenceMetric.xml
 %feature("docstring") RelativeDifferenceMetric "
@@ -2400,7 +2406,7 @@ Returns default metric name.
 
 
 // File: namespacepyfmt.xml
-%feature("docstring")  pyfmt::scriptPreamble "std::string pyfmt::scriptPreamble()
+%feature("docstring")  pyfmt::preambled "std::string pyfmt::preambled(const std::string &code)
 ";
 
 %feature("docstring")  pyfmt::printBool "std::string pyfmt::printBool(double value)
diff --git a/auto/Wrap/libBornAgainCore.py b/auto/Wrap/libBornAgainCore.py
index 41bfd1defca874425258fe94dfcce45a500d4b56..3fa2db3ef9eb044c59b94bbe7f06d094c8c04a66 100644
--- a/auto/Wrap/libBornAgainCore.py
+++ b/auto/Wrap/libBornAgainCore.py
@@ -3449,11 +3449,19 @@ class QSpecScan(object):
         return _libBornAgainCore.QSpecScan_setAbsoluteQResolution(self, *args)
 
     def setOffset(self, offset):
-        r"""setOffset(QSpecScan self, double offset)"""
+        r"""
+        setOffset(QSpecScan self, double offset)
+        void QSpecScan::setOffset(double offset)
+
+        """
         return _libBornAgainCore.QSpecScan_setOffset(self, offset)
 
     def offset(self):
-        r"""offset(QSpecScan self) -> double"""
+        r"""
+        offset(QSpecScan self) -> double
+        double QSpecScan::offset() const
+
+        """
         return _libBornAgainCore.QSpecScan_offset(self)
 
 # Register QSpecScan in _libBornAgainCore:
diff --git a/auto/Wrap/libBornAgainCore_wrap.cpp b/auto/Wrap/libBornAgainCore_wrap.cpp
index 4c09c1ce5fe7f533e6ac4df0be6e3201fd6b31c6..ae87a6c4ce4a752f94043a257d5a53163de6e372 100644
--- a/auto/Wrap/libBornAgainCore_wrap.cpp
+++ b/auto/Wrap/libBornAgainCore_wrap.cpp
@@ -44602,8 +44602,16 @@ static PyMethodDef SwigMethods[] = {
 		"Sets qz resolution values via IRangedDistribution and values of standard deviations.  std_dev can be either single-valued or a numpy array. In the latter case the length of the array should coinside with the length of the qz-axis. \n"
 		"\n"
 		""},
-	 { "QSpecScan_setOffset", _wrap_QSpecScan_setOffset, METH_VARARGS, "QSpecScan_setOffset(QSpecScan self, double offset)"},
-	 { "QSpecScan_offset", _wrap_QSpecScan_offset, METH_O, "QSpecScan_offset(QSpecScan self) -> double"},
+	 { "QSpecScan_setOffset", _wrap_QSpecScan_setOffset, METH_VARARGS, "\n"
+		"QSpecScan_setOffset(QSpecScan self, double offset)\n"
+		"void QSpecScan::setOffset(double offset)\n"
+		"\n"
+		""},
+	 { "QSpecScan_offset", _wrap_QSpecScan_offset, METH_O, "\n"
+		"QSpecScan_offset(QSpecScan self) -> double\n"
+		"double QSpecScan::offset() const\n"
+		"\n"
+		""},
 	 { "QSpecScan_swigregister", QSpecScan_swigregister, METH_O, NULL},
 	 { "QSpecScan_swiginit", QSpecScan_swiginit, METH_VARARGS, NULL},
 	 { "delete_ISimulation", _wrap_delete_ISimulation, METH_O, "\n"
diff --git a/devtools/code/normalize-usercode.py b/devtools/code/normalize-usercode.py
index c218be07171c7d1aa5421373b91836b2a31453a8..36258bd97ca5c29fb4df8949dfedeaa82a754f87 100755
--- a/devtools/code/normalize-usercode.py
+++ b/devtools/code/normalize-usercode.py
@@ -55,7 +55,7 @@ def cycle_text(ti, fname):
     Returns normalized version of script ti as obtained from BornAgain's export-to-Python function.
     """
     s = retrieve_simulation(ti, fname)
-    return ba.generateSimulationCode(s)
+    return ba.simulationPlotCode(s)
 
 
 def normalize_text(ti, fname):
@@ -92,10 +92,9 @@ def normalize_file(fname, inplace):
                 print(f'.. read {len(ti.split())} lines')
 
         m = re.search(r"""if __name__ == '__main__':
-    result = run_simulation\(\)
-    ba.plot_simulation_result\(result.*?\)""", ti)
+    ba.run_and_plot\(get_simulation\(get_sample\(\)\)\)""", ti)
         if not m:
-            print('=> non-standard => SKIPPED')
+            print('=> main is not run_and_plot => SKIPPED')
             return 3
 
         # normalize