diff --git a/Base/Axis/IAxis.cpp b/Base/Axis/IAxis.cpp
index 5b56e268500e76560ec4944f93193ae9de427001..051c9ba8a9b9cc29ff4f1554e8c02e314ffb372e 100644
--- a/Base/Axis/IAxis.cpp
+++ b/Base/Axis/IAxis.cpp
@@ -37,3 +37,7 @@ bool IAxis::contains(double value) const {
 double IAxis::span() const {
     return upperBound() - lowerBound();
 }
+
+double IAxis::center() const {
+    return (upperBound() + lowerBound())/2;
+}
diff --git a/Base/Axis/IAxis.h b/Base/Axis/IAxis.h
index 21154ced4bad5d5a51e53f486a639e6e02817b53..f4ba5bf2186e587853d7130b6ff48dac95409ff3 100644
--- a/Base/Axis/IAxis.h
+++ b/Base/Axis/IAxis.h
@@ -56,6 +56,9 @@ public:
     //! Returns distance from first to last point
     double span() const;
 
+    //! Returns midpoint of axis
+    double center() const;
+
     virtual double binCenter(size_t index) const = 0;
 
     //! find bin index which is best match for given value
diff --git a/CHANGELOG b/CHANGELOG
index 4778aa9ea7bf5bf777f1ee67b0c2af1c56e0d06e..fa35811f78963848abfb40e1e60d0c43a1160502 100644
--- a/CHANGELOG
+++ b/CHANGELOG
@@ -1,10 +1,26 @@
 BornAgain-1.18.99, ongoing development
   > API changes:
-    1) Lattice has been renamed to Lattice3D, for clearer disambuation from Lattice2D.
-    2) Retain only those interference functions constructors that take a Lattice2D
-       or Lattice3D parameter.
+    * Compact constructors like GISASimulation(beam, sample, detector)
+    * Removed "get" from several accessor functions
+    * Class Instrument removed  from Python API (instead use simulation.beam() etc)
+    * Lattice has been renamed to Lattice3D, for clearer disambuation from Lattice2D
+    * Removed some interference functions constructors
+    * Remove old R&T computations from API (now in Code/Legacy and Tests)
+    * Python plot API entirely keyword based
   > Fixes of unreported bugs:
-    * For alpha_i, set scattered intensity to 0 (incorrect, but everything else would be worse)
+    * Several GUI bugs that caused crashes
+    * For alpha_i=0, set scattered intensity to 0
+  > Examples:
+    * Uniform API usage in line with exported Python code
+  > Internal refactoring:
+    * Moved code to break directory include cycles
+    * Renamed a few classes (ISimulation)
+    * Simplified and commented std test machinery
+    * Simplified iteration over INode children
+    * Simplified export of get_sample function
+    * Simplified data read/write wrapper
+    * Replaced boost::filesystem by C++17 std
+    * IComputeFF and children drawn out of IFormFactor inheritance tree
 
 BornAgain-1.18.0, released 2020.10.30
   > API changes:
@@ -33,6 +49,7 @@ BornAgain-1.18.0, released 2020.10.30
     6) Scalar Fresnel computation unified with polarized implementation
     7) Directory Core broken into a hierarchy: Core, Device, Sample, Param, Base
     8) Request C++17
+    9) Error handling: core decoupled from Qt
   > Issues from external (GitHub) tracker:
     * #918:  Specular peak location wrong when using flat detector
     * #946:  CTest points to broken visualization scripts
diff --git a/Core/Export/SimulationToPython.cpp b/Core/Export/SimulationToPython.cpp
index 9257436cac5a5d5f7acad4059f418bb5fe2a28bd..178fb875ed9c8484d86b323e1af3a6f513f06c94 100644
--- a/Core/Export/SimulationToPython.cpp
+++ b/Core/Export/SimulationToPython.cpp
@@ -28,6 +28,7 @@
 #include "Core/Simulation/SpecularSimulation.h"
 #include "Device/Beam/FootprintGauss.h"
 #include "Device/Beam/FootprintSquare.h"
+#include "Device/Detector/DetectorUtils.h"
 #include "Device/Detector/RectangularDetector.h"
 #include "Device/Detector/RegionOfInterest.h"
 #include "Device/Detector/SphericalDetector.h"
@@ -127,13 +128,22 @@ std::string defineDetector(const ISimulation* simulation) {
     result << std::setprecision(12);
 
     if (const auto* const det = dynamic_cast<const SphericalDetector*>(detector)) {
-        result << indent() << "detector = ba.SphericalDetector(";
-        for (size_t index = 0; index < det->dimension(); ++index) {
-            if (index != 0)
-                result << ", ";
-            result << det->axis(index).size() << ", "
-                   << pyfmt::printDegrees(det->axis(index).lowerBound()) << ", "
-                   << pyfmt::printDegrees(det->axis(index).upperBound());
+        ASSERT(det->dimension()==2);
+        if (DetectorUtils::isQuadratic(*det)) {
+            result << indent() << "nbin = " << det->axis(0).size() << "\n";
+            result << indent() << "detector = ba.SphericalDetector(nbin, "
+                   << pyfmt::printDegrees(det->axis(0).span()) << ", "
+                   << pyfmt::printDegrees(det->axis(0).center()) << ", "
+                   << pyfmt::printDegrees(det->axis(1).center());
+        } else {
+            result << indent() << "nx = " << det->axis(0).size() << "\n";
+            result << indent() << "ny = " << det->axis(1).size() << "\n";
+            result << indent() << "detector = ba.SphericalDetector(nx, "
+                   << pyfmt::printDegrees(det->axis(0).lowerBound()) << ", "
+                   << pyfmt::printDegrees(det->axis(0).upperBound()) << ", "
+                   << "ny , "
+                   << pyfmt::printDegrees(det->axis(1).lowerBound()) << ", "
+                   << pyfmt::printDegrees(det->axis(1).upperBound());
         }
         result << ")\n";
     } else if (const auto* const det = dynamic_cast<const RectangularDetector*>(detector)) {
diff --git a/Core/Simulation/ISimulation.cpp b/Core/Simulation/ISimulation.cpp
index 8ae36f605dfa460e8b4cac55ed942eb0d48dac27..402d1fc08121dbaa13a6435c794caa4008b91bc8 100644
--- a/Core/Simulation/ISimulation.cpp
+++ b/Core/Simulation/ISimulation.cpp
@@ -254,7 +254,6 @@ void ISimulation::addParameterDistribution(const std::string& param_name,
 }
 
 void ISimulation::addParameterDistribution(const ParameterDistribution& par_distr) {
-    std::cout << "DEBUG ISimulation::addParameterDistribution" << std::endl;
     validateParametrization(par_distr);
     m_distribution_handler.addParameterDistribution(par_distr);
 }
diff --git a/Device/Detector/DetectorUtils.cpp b/Device/Detector/DetectorUtils.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..fdf1a5af6fbc0bc76ecef5ff515490e42dbc8797
--- /dev/null
+++ b/Device/Detector/DetectorUtils.cpp
@@ -0,0 +1,26 @@
+//  ************************************************************************************************
+//
+//  BornAgain: simulate and fit scattering at grazing incidence
+//
+//! @file      Device/Detector/DetectorUtils.cpp
+//! @brief     Implements namespace DetectorUtils.
+//!
+//! @homepage  http://www.bornagainproject.org
+//! @license   GNU General Public License v3 or higher (see COPYING)
+//! @copyright Forschungszentrum Jülich GmbH 2018
+//! @authors   Scientific Computing Group at MLZ (see CITATION, AUTHORS)
+//
+//  ************************************************************************************************
+
+#include "Device/Detector/DetectorUtils.h"
+#include "Device/Detector/IDetector2D.h"
+
+bool DetectorUtils::isQuadratic(const IDetector2D& det) {
+    ASSERT(det.dimension()==2);
+    if (det.axis(0).size()!=det.axis(1).size())
+        return false;
+    if (std::abs(det.axis(0).span() - det.axis(1).span()) >
+        1e-12*(det.axis(0).span() + det.axis(1).span()))
+        return false;
+    return true;
+}
diff --git a/Device/Detector/DetectorUtils.h b/Device/Detector/DetectorUtils.h
new file mode 100644
index 0000000000000000000000000000000000000000..b23aa3b9ad6bf2ef16592d40460666c0b53b9cd3
--- /dev/null
+++ b/Device/Detector/DetectorUtils.h
@@ -0,0 +1,26 @@
+//  ************************************************************************************************
+//
+//  BornAgain: simulate and fit scattering at grazing incidence
+//
+//! @file      Device/Detector/DetectorUtils.h
+//! @brief     Defines namespace DetectorUtils.
+//!
+//! @homepage  http://www.bornagainproject.org
+//! @license   GNU General Public License v3 or higher (see COPYING)
+//! @copyright Forschungszentrum Jülich GmbH 2018
+//! @authors   Scientific Computing Group at MLZ (see CITATION, AUTHORS)
+//
+//  ************************************************************************************************
+
+#ifndef BORNAGAIN_DEVICE_DETECTOR_DETECTORUTILS_H
+#define BORNAGAIN_DEVICE_DETECTOR_DETECTORUTILS_H
+
+class IDetector2D;
+
+namespace DetectorUtils {
+
+bool isQuadratic(const IDetector2D& det);
+
+} // namespace DetectorUtils
+
+#endif // BORNAGAIN_DEVICE_DETECTOR_DETECTORUTILS_H
diff --git a/Device/Detector/SphericalDetector.cpp b/Device/Detector/SphericalDetector.cpp
index 5266e1c3f51215735ccb01b575e674722062fac0..24ddc90549ed201f5c422239cd92392523e5d913 100644
--- a/Device/Detector/SphericalDetector.cpp
+++ b/Device/Detector/SphericalDetector.cpp
@@ -30,6 +30,10 @@ SphericalDetector::SphericalDetector(size_t n_phi, double phi_min, double phi_ma
     setDetectorParameters(n_phi, phi_min, phi_max, n_alpha, alpha_min, alpha_max);
 }
 
+SphericalDetector::SphericalDetector(size_t n_bin, double width, double phi, double alpha)
+    : SphericalDetector(n_bin, phi-width/2, phi+width/2, n_bin, alpha-width/2, alpha+width/2) {
+}
+
 SphericalDetector::SphericalDetector(const SphericalDetector& other) : IDetector2D(other) {
     setName("SphericalDetector");
 }
diff --git a/Device/Detector/SphericalDetector.h b/Device/Detector/SphericalDetector.h
index 256034574065c3d87fd621092b07975731e226e8..258180ab0c6da8e05c317fb5fa7a2c8b5d374507 100644
--- a/Device/Detector/SphericalDetector.h
+++ b/Device/Detector/SphericalDetector.h
@@ -37,6 +37,13 @@ public:
     SphericalDetector(size_t n_phi, double phi_min, double phi_max, size_t n_alpha,
                       double alpha_min, double alpha_max);
 
+    //! Spherical detector constructor with quadratic angle ranges
+    //! @param n_bin number of bins per direction
+    //! @param width angular width
+    //! @param phi   central phi angle
+    //! @param alpha central alpha angle
+    SphericalDetector(size_t n_bin, double width, double phi, double alpha);
+
     SphericalDetector(const SphericalDetector& other);
 
     SphericalDetector* clone() const override;
diff --git a/Examples/Python/fit54_ExternalMinimizer/lmfit_with_plotting.py b/Examples/Python/fit54_ExternalMinimizer/lmfit_with_plotting.py
index 83a6aac72ee922d771be3204e5cfe3a5c6d93ff8..4871a86cb04063cc6a960c3619afb008658d8eba 100644
--- a/Examples/Python/fit54_ExternalMinimizer/lmfit_with_plotting.py
+++ b/Examples/Python/fit54_ExternalMinimizer/lmfit_with_plotting.py
@@ -27,7 +27,7 @@ def get_sample(params):
     particle_layout.addParticle(sphere)
 
     interference = ba.InterferenceFunction2DLattice(
-        ba.HexagonalLattice2D(lattice_length))
+        ba.HexagonalLattice2D(lattice_length, 0*deg))
     pdf = ba.FTDecayFunction2DCauchy(10*nm, 10*nm, 0)
     interference.setDecayFunction(pdf)
 
@@ -74,13 +74,13 @@ def create_real_data():
     return noisy
 
 
-class Plotter:
+class LMFITPlotter:
     """
     Adapts standard plotter for lmfit minimizer.
     """
     def __init__(self, fit_objective, every_nth=10):
         self.fit_objective = fit_objective
-        self.plotter_gisas = ba.PlotterGISAS()
+        self.plotter_gisas = ba.fit_monitor.PlotterGISAS()
         self.every_nth = every_nth
 
     def __call__(self, params, iter, resid):
@@ -102,7 +102,7 @@ def run_fitting():
     params.add('radius', value=7*nm, min=5*nm, max=8*nm)
     params.add('length', value=10*nm, min=8*nm, max=14*nm)
 
-    plotter = Plotter(fit_objective)
+    plotter = LMFITPlotter(fit_objective)
     result = lmfit.minimize(fit_objective.evaluate_residuals,
                             params,
                             iter_cb=plotter)
diff --git a/Examples/Python/fit55_SpecularIntro/FitSpecularBasics.py b/Examples/Python/fit55_SpecularIntro/FitSpecularBasics.py
index 249ee9daf9f730b0463272d5eee7bbd50626efba..4ce6cb3e1793f1199d7e5f57170c4c2a6d2c20ff 100644
--- a/Examples/Python/fit55_SpecularIntro/FitSpecularBasics.py
+++ b/Examples/Python/fit55_SpecularIntro/FitSpecularBasics.py
@@ -110,7 +110,7 @@ def run_fitting():
     fit_objective = ba.FitObjective()
     fit_objective.addSimulationAndData(get_simulation, real_data, 1.0)
 
-    plot_observer = ba.PlotterSpecular()
+    plot_observer = ba.fit_monitor.PlotterSpecular()
     fit_objective.initPrint(10)
     fit_objective.initPlot(10, plot_observer)
 
diff --git a/Examples/Python/fit55_SpecularIntro/FitWithUncertainties.py b/Examples/Python/fit55_SpecularIntro/FitWithUncertainties.py
index c8f8a4ba6719cb9e20c674f2d1f5af8ad68af098..f8eb0d58963922e2bb96675d55d982c7edaa1915 100644
--- a/Examples/Python/fit55_SpecularIntro/FitWithUncertainties.py
+++ b/Examples/Python/fit55_SpecularIntro/FitWithUncertainties.py
@@ -16,7 +16,6 @@ import bornagain as ba
 from matplotlib import pyplot as plt
 from os import path
 
-
 def get_sample(params):
     """
     Creates a sample and returns it
@@ -118,7 +117,7 @@ def run_fitting():
     fit_objective = ba.FitObjective()
     fit_objective.addSimulationAndData(get_simulation, real_data, uncertainties)
 
-    plot_observer = ba.PlotterSpecular(units=ba.Axes.RQ4)
+    plot_observer = ba.fit_monitor.PlotterSpecular(units=ba.Axes.RQ4)
     fit_objective.initPrint(10)
     fit_objective.initPlot(10, plot_observer)
     fit_objective.setObjectiveMetric("Chi2", "L1")
diff --git a/Examples/Python/sim01_Particles/AllFormFactorsAvailable.py b/Examples/Python/sim01_Particles/AllFormFactorsAvailable.py
index 01330432da73abc2645eeb68dcc03de998fe8140..82ce0f5f93d951b5bd5a01ef4f8660ad91d78f42 100644
--- a/Examples/Python/sim01_Particles/AllFormFactorsAvailable.py
+++ b/Examples/Python/sim01_Particles/AllFormFactorsAvailable.py
@@ -78,7 +78,7 @@ def simulate(ff):
     return simulation.result()
 
 
-def run_simulation():
+def simulate_and_plot():
     """
     Run simulation one by one for every form factor from the list and plot results
     on a single canvas
@@ -113,8 +113,8 @@ def run_simulation():
                  horizontalalignment='center',
                  verticalalignment='center',
                  fontsize=9)
+    plt.show()
 
 
 if __name__ == '__main__':
-    run_simulation()
-    plt.show()
+    simulate_and_plot()
diff --git a/Examples/Python/sim01_Particles/CylindersAndPrisms.py b/Examples/Python/sim01_Particles/CylindersAndPrisms.py
index cbab1b3b6d47c381c6f2adc6dc66d69f608dbf28..795c4f2f1adfa97b9fb120f57f53beb6c2e6fd13 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))
+    nbin = 100
+    detector = ba.SphericalDetector(nbin, 2.0*deg, 0.0*deg, 1.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..4ff161164ef0215c5dca4e1f6a52f23215042dfe 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,17 @@ 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))
+    nx = 200
+    ny = 200
+    detector = ba.SphericalDetector(nx, -2.0*deg, 2.0*deg, ny, 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..73a903eb021ed19559b9200433a9969298bdcdf7 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,17 @@ 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))
+    nx = 200
+    ny = 200
+    detector = ba.SphericalDetector(nx, -2.0*deg, 2.0*deg, ny, 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..81e115809deace253452a6de7e9d30a74f5dcd2d 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))
+    nbin = 200
+    detector = ba.SphericalDetector(nbin, 2.0*deg, 1.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/sim01_Particles/RotatedPyramids.py b/Examples/Python/sim01_Particles/RotatedPyramids.py
index 0077becb7bfe0fbb2de1ba8ed717227810b81960..7c83b39452ae51fb227d738aa9e72b028473d3ca 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,17 @@ 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))
+    nx = 200
+    ny = 200
+    detector = ba.SphericalDetector(nx, -2.0*deg, 2.0*deg, ny, 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..6fa2b1f0786120ee50aaf6fa06923ddf39962b0e 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))
+    nbin = 200
+    detector = ba.SphericalDetector(nbin, 2.0*deg, 1.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/BiMaterialCylinders.py b/Examples/Python/sim02_Complexes/BiMaterialCylinders.py
index c99b9de98ba2a6c4bfefebcc54c150e46627ea53..6b4f501baa909f97000aba19f199337fa625f3bd 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))
+    nbin = 100
+    detector = ba.SphericalDetector(nbin, 2.0*deg, 0.0*deg, 1.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..c6452b607b2fee7671d9b6c15dae326918b202fd 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))
+    nbin = 200
+    detector = ba.SphericalDetector(nbin, 2.0*deg, 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/HexagonalLatticesWithBasis.py b/Examples/Python/sim02_Complexes/HexagonalLatticesWithBasis.py
index 860732d97d96c21cdd8a05304edf2b4bc84b9855..4fc0c6d8f666cb69e5254176103c77e7ae08e283 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,17 @@ 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))
+    nx = 200
+    ny = 200
+    detector = ba.SphericalDetector(nx, -1.0*deg, 1.0*deg, ny, 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/LargeParticlesFormFactor.py b/Examples/Python/sim02_Complexes/LargeParticlesFormFactor.py
index 040309946ec011bfd50c7af3fd36c138d0d6880c..2f2e28ab0a91c5b250a2f2df6bb6075eb8237196 100644
--- a/Examples/Python/sim02_Complexes/LargeParticlesFormFactor.py
+++ b/Examples/Python/sim02_Complexes/LargeParticlesFormFactor.py
@@ -54,7 +54,7 @@ def get_simulation(integration_flag):
     return simulation
 
 
-def run_simulation():
+def simulate_and_plot():
     """
     Run simulation and plot results 4 times: for small and large cylinders,
     with and without integration
@@ -108,8 +108,8 @@ def run_simulation():
         zmin = condition['zmin']
         zmax = condition['zmax']
         ba.plot_colormap(result,
-                         zmin=zmin,
-                         zmax=zmax,
+                         intensity_min=zmin,
+                         intensity_max=zmax,
                          cmap='jet',
                          aspect='auto')
 
@@ -119,8 +119,8 @@ def run_simulation():
                  horizontalalignment='center',
                  verticalalignment='center',
                  fontsize=12)
+    plt.show()
 
 
 if __name__ == '__main__':
-    run_simulation()
-    plt.show()
+    simulate_and_plot()
diff --git a/Examples/Python/sim02_Complexes/MesoCrystal.py b/Examples/Python/sim02_Complexes/MesoCrystal.py
index 3bd203daadd04ab6d752e726b4b8dc616ddf2e31..a3b78a2461fa0fca832355d4ac13c5922eb7a281 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,17 @@ 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))
+    nx = 200
+    ny = 200
+    detector = ba.SphericalDetector(nx, -2.0*deg, 2.0*deg, ny, 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..493b78308c09cdfafad51b77a9dc6a5c0d254405 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))
+    nbin = 100
+    detector = ba.SphericalDetector(nbin, 2.0*deg, 0.0*deg, 1.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..bdac8aa52541d7ea9192a3687d49ee01fdf613f0 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))
+    nbin = 200
+    detector = ba.SphericalDetector(nbin, 2.0*deg, 1.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/ApproximationLMA.py b/Examples/Python/sim03_Structures/ApproximationLMA.py
index 73c8eebb9b3947f23b9da1de47a4608d8b17e872..158032846f186d29d70d70a7d1b1217a312ac908 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))
+    nbin = 200
+    detector = ba.SphericalDetector(nbin, 2.0*deg, 1.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/ApproximationSSCA.py b/Examples/Python/sim03_Structures/ApproximationSSCA.py
index abc4a8cb83ba6dd9afedf41a54b0aa6ac54e53cf..8530cfe49f7126a8a1a4627bd4f1eb11c71085f1 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))
+    nbin = 200
+    detector = ba.SphericalDetector(nbin, 2.0*deg, 1.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/CosineRipplesAtRectLattice.py b/Examples/Python/sim03_Structures/CosineRipplesAtRectLattice.py
index b4ae226890224b83dc62439efe2f945e434d0c20..fc0993f33e756073b3e876780472505ffd71a4b5 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,17 @@ 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))
+    nx = 100
+    ny = 100
+    detector = ba.SphericalDetector(nx, -1.5*deg, 1.5*deg, ny, 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/FindPeaks.py b/Examples/Python/sim03_Structures/FindPeaks.py
index 73d5d3cc323d699c82bce0a2d72721adc9d43647..2ec3cf1487c105693f9d770eceb3fe3c993c7337 100644
--- a/Examples/Python/sim03_Structures/FindPeaks.py
+++ b/Examples/Python/sim03_Structures/FindPeaks.py
@@ -77,7 +77,7 @@ def run_simulation():
     return simulation.result()
 
 
-if __name__ == '__main__':
+def simulate_and_plot():
     result = run_simulation().histogram2d()
     ba.plot_histogram(result, cmap='jet', aspect='auto')
 
@@ -92,3 +92,7 @@ if __name__ == '__main__':
              color='white',
              markersize=10)
     plt.show()
+
+
+if __name__ == '__main__':
+    simulate_and_plot()
diff --git a/Examples/Python/sim03_Structures/Interference1DRadialParaCrystal.py b/Examples/Python/sim03_Structures/Interference1DRadialParaCrystal.py
index 12e550a5a155149517eeca0327417769bcc95125..a411314da8651bc490d2d319d1207db2f37b1854 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,17 @@ 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))
+    nx = 200
+    ny = 200
+    detector = ba.SphericalDetector(nx, -2.0*deg, 2.0*deg, ny, 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..5b5bb61fbc349340c5e105fe4cec728a9c28f5d7 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,17 @@ 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))
+    nx = 200
+    ny = 200
+    detector = ba.SphericalDetector(nx, -2.0*deg, 2.0*deg, ny, 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..fb5b18b520bc400741edd71f58a918e734826739 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))
+    nbin = 100
+    detector = ba.SphericalDetector(nbin, 2.0*deg, 1.0*deg, 1.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..86ad4e6daaa1b431755e70a60b449b0d75d57d0e 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,17 @@ 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))
+    nx = 200
+    ny = 200
+    detector = ba.SphericalDetector(nx, -2.0*deg, 2.0*deg, ny, 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..2b1da4c2d661ff728dfd2935238c9171ad91abe6 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,17 @@ 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))
+    nx = 200
+    ny = 200
+    detector = ba.SphericalDetector(nx, -2.0*deg, 2.0*deg, ny, 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..99a9af46b441b04bd759bedac75df846fb75936d 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,17 @@ 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))
+    nx = 200
+    ny = 200
+    detector = ba.SphericalDetector(nx, -2.0*deg, 2.0*deg, ny, 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..e7968d38fde37fbabf5b2cb8d9e68f073be7f20d 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,17 @@ 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))
+    nx = 200
+    ny = 200
+    detector = ba.SphericalDetector(nx, -2.0*deg, 2.0*deg, ny, 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..b720c9981dbfdd86dfd67314475bf804fccdde90 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,18 @@ 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))
+    nx = 200
+    ny = 200
+    detector = ba.SphericalDetector(nx, -0.5*deg, 0.5*deg, ny, 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..cfff1d0eefef5a98e2e95b6f4629b64b43c88be1 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,17 @@ 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))
+    nx = 200
+    ny = 200
+    detector = ba.SphericalDetector(nx, -1.0*deg, 1.0*deg, ny, 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..4082acffcbb507268762d67d7fa7980939360a81 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,17 @@ 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))
+    nx = 200
+    ny = 200
+    detector = ba.SphericalDetector(nx, -1.5*deg, 1.5*deg, ny, 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..db574d201f13144ac46ebab2715c3575c6283af8 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))
+    nbin = 200
+    detector = ba.SphericalDetector(nbin, 2.0*deg, 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/CorrelatedRoughness.py b/Examples/Python/sim04_Multilayers/CorrelatedRoughness.py
index a32410384e7c997082411d3a6dbf060feee46284..f968871ae5b3ff1fc5792df6b6679ec0641bc6ad 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))
+    nbin = 200
+    detector = ba.SphericalDetector(nbin, 1.0*deg, 0.0*deg, 0.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/CylindersInAverageLayer.py b/Examples/Python/sim04_Multilayers/CylindersInAverageLayer.py
index bfbe4d269e5e858c40f11980d172608c49d41ac1..dcbca06bfa726862b80c00ce790004f56d8deb34 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,18 @@ 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))
+    nx = 100
+    ny = 100
+    detector = ba.SphericalDetector(nx, -2.0*deg, 2.0*deg, ny, 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/GratingMC.py b/Examples/Python/sim04_Multilayers/GratingMC.py
index 91cb58dfa3f6c1264b57abec931f3bc151e8ae57..8e1a89e89f7f2644940cc15419410cab96a34dd0 100644
--- a/Examples/Python/sim04_Multilayers/GratingMC.py
+++ b/Examples/Python/sim04_Multilayers/GratingMC.py
@@ -77,7 +77,7 @@ def run_simulation():
     return simulation.result()
 
 
-if __name__ == '__main__':
+def simulate_and_plot():
     interactive = True
     result = run_simulation().histogram2d()
     ba.plot_histogram(result, cmap='jet', aspect='auto')
@@ -93,3 +93,7 @@ if __name__ == '__main__':
              color='white',
              markersize=10)
     plt.show()
+
+
+if __name__ == '__main__':
+    simulate_and_plot()
diff --git a/Examples/Python/sim04_Multilayers/HalfSpheresInAverageTopLayer.py b/Examples/Python/sim04_Multilayers/HalfSpheresInAverageTopLayer.py
index a691fbab52d14a78ea9c97b1a90129aade4781e1..5c5544d3d3a8ef7dd4d708d0d173f1d318f184db 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,18 @@ 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))
+    nx = 100
+    ny = 100
+    detector = ba.SphericalDetector(nx, -2.0*deg, 2.0*deg, ny, 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..6d1e01bbaca188f25b3c16db0e28a4293dacc4c5 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)
+    nbin = 200
+    detector = ba.SphericalDetector(nbin, 6.0*deg, 0.0*deg, 3.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..088a9ab3e45c52bca19022ca76c3090183c25240 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))
+    nbin = 100
+    detector = ba.SphericalDetector(nbin, 2.0*deg, 1.0*deg, 1.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..0ca4ccd72eb3085323578b0bbfc5179ebc12d43e 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))
+    nbin = 100
+    detector = ba.SphericalDetector(nbin, 2.0*deg, 1.0*deg, 1.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..739ca5512b4303f4798ba39309705904001023b0 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))
+    nbin = 100
+    detector = ba.SphericalDetector(nbin, 2.0*deg, 1.0*deg, 1.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..6845d03748d0ce08a107df721957e1438bfca0ff 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,19 @@ 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))
+    nx = 101
+    ny = 101
+    detector = ba.SphericalDetector(nx, -2.0*deg, 2.0*deg, ny, 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()
 
diff --git a/Examples/Python/sim31_Parameterization/AccessingSimulationResults.py b/Examples/Python/sim31_Parameterization/AccessingSimulationResults.py
index b1e8ad350cd202d2ba4f521599eaecff9ff941ad..4f77425a76831646a8b284da1b91842ec1635b21 100644
--- a/Examples/Python/sim31_Parameterization/AccessingSimulationResults.py
+++ b/Examples/Python/sim31_Parameterization/AccessingSimulationResults.py
@@ -71,15 +71,13 @@ def get_noisy_image(hist):
     return result
 
 
-def plot_histogram(hist, zmin=None, zmax=None):
+def plot_histogram(hist, **kwargs):
     ba.plot_histogram(hist,
                       xlabel=r'$\varphi_f ^{\circ}$',
                       ylabel=r'$\alpha_f ^{\circ}$',
                       zlabel="",
-                      zmin=zmin,
-                      zmax=zmax,
                       cmap='jet',
-                      aspect='auto')
+                      aspect='auto', **kwargs)
 
 
 def get_relative_difference(hist):
@@ -140,7 +138,7 @@ def plot(hist):
 
     plt.subplot(2, 2, 3)
     reldiff = get_relative_difference(hist)
-    plot_histogram(reldiff, zmin=1e-03, zmax=10)
+    plot_histogram(reldiff, intensity_min=1e-03, intensity_max=10)
     plt.title("Relative difference")
 
     plt.subplot(2, 2, 4)
diff --git a/Examples/Python/sim29_DepthProbe/DepthProbe.py b/Examples/Python/special40_DepthProbe/DepthProbe.py
similarity index 100%
rename from Examples/Python/sim29_DepthProbe/DepthProbe.py
rename to Examples/Python/special40_DepthProbe/DepthProbe.py
diff --git a/GUI/coregui/Models/GUIDomainSampleVisitor.cpp b/GUI/coregui/Models/GUIDomainSampleVisitor.cpp
index 996aacc2258e04862450ec0f279676dad9e8d476..f1d6434977bd7e9f498b69cf0e7acf7348e5ca00 100644
--- a/GUI/coregui/Models/GUIDomainSampleVisitor.cpp
+++ b/GUI/coregui/Models/GUIDomainSampleVisitor.cpp
@@ -78,7 +78,7 @@ SessionItem* GUIDomainSampleVisitor::populateSampleModel(SampleModel* sampleMode
     if (m_topSampleName.isEmpty())
         m_topSampleName = sample.getName().c_str();
 
-    for (const auto [child, depth, parent] : NodeUtils::progenyPlus(&sample)) {
+    for (const auto& [child, depth, parent] : NodeUtils::progenyPlus(&sample)) {
         setDepth(depth + 1);
         child->accept(this);
     }
@@ -106,7 +106,7 @@ void GUIDomainSampleVisitor::visit(const Layer* sample) {
     SessionItem* parent = m_levelToParentItem[depth() - 1];
     ASSERT(parent);
 
-    auto multilayer = dynamic_cast<const MultiLayer*>(m_itemToSample[parent]);
+    const auto* multilayer = dynamic_cast<const MultiLayer*>(m_itemToSample[parent]);
     ASSERT(multilayer);
     size_t layer_index = MultiLayerUtils::IndexOfLayer(*multilayer, sample);
     const LayerInterface* top_interface =
diff --git a/Param/Node/NodeUtils.cpp b/Param/Node/NodeUtils.cpp
index 40e8adc60590e2385670e8512ac2e42c930d9f09..efac25fa51986a9d2f0cb0f133a1685e47cbd91a 100644
--- a/Param/Node/NodeUtils.cpp
+++ b/Param/Node/NodeUtils.cpp
@@ -59,12 +59,16 @@ std::string nodeString(const INode& node, int depth) {
 }
 } // namespace
 
+//  ************************************************************************************************
+//  namespace NodeUtils
+//  ************************************************************************************************
+
 std::vector<std::tuple<const INode*, int, const INode*>> NodeUtils::progenyPlus(const INode* node,
                                                                                 int level) {
     std::vector<std::tuple<const INode*, int, const INode*>> result;
     result.push_back({node, level, nullptr});
     for (const auto* child : node->getChildren()) {
-        for (const auto [subchild, sublevel, subparent] : progenyPlus(child, level + 1))
+        for (const auto& [subchild, sublevel, subparent] : progenyPlus(child, level + 1))
             result.push_back({subchild, sublevel, child});
     }
     return result;
@@ -72,7 +76,7 @@ std::vector<std::tuple<const INode*, int, const INode*>> NodeUtils::progenyPlus(
 
 std::string NodeUtils::nodeToString(const INode* node) {
     std::ostringstream result;
-    for (const auto [child, depth, parent] : progenyPlus(node))
+    for (const auto& [child, depth, parent] : progenyPlus(node))
         result << nodeString(*child, depth);
     return result.str();
 }
diff --git a/Tests/Functional/Python/PyCore/CMakeLists.txt b/Tests/Functional/Python/PyCore/CMakeLists.txt
index 9d195c49ec556d2b9f26b5c25a5f7fb7dd0c69a6..ed485d869195f64871f6b7fbb27ed88ec86b39b1 100644
--- a/Tests/Functional/Python/PyCore/CMakeLists.txt
+++ b/Tests/Functional/Python/PyCore/CMakeLists.txt
@@ -5,18 +5,17 @@
 set(OUTPUT_DIR ${TEST_OUTPUT_DIR_PY_CORE})
 file(MAKE_DIRECTORY ${OUTPUT_DIR})
 
-file(GLOB sources RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} "*.py")
+file(GLOB tests RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} "*.py")
+list(REMOVE_ITEM tests utils.py)
 if(NOT BORNAGAIN_TIFF_SUPPORT)
-    list(REMOVE_ITEM sources "intensitydata_io_tiff.py")
+    list(REMOVE_ITEM tests "intensitydata_io_tiff.py")
 endif()
 
-foreach(_src ${sources})
-    configure_file(${CMAKE_CURRENT_SOURCE_DIR}/${_src} ${OUTPUT_DIR}/${_src} @ONLY)
-endforeach()
-
-set(tests ${sources})
-list(REMOVE_ITEM tests utils.py)
+configure_file(${CMAKE_CURRENT_SOURCE_DIR}/utils.py ${OUTPUT_DIR}/utils.py @ONLY)
 
 foreach(_test ${tests})
-    add_test(PyCore.${_test} ${Python3_EXECUTABLE} ${OUTPUT_DIR}/${_test})
+    configure_file(${CMAKE_CURRENT_SOURCE_DIR}/${_test} ${OUTPUT_DIR}/${_test} COPYONLY)
+    add_test(PyCore.${_test}
+        env PYTHONPATH=${CMAKE_LIBRARY_OUTPUT_DIRECTORY}
+        ${Python3_EXECUTABLE} ${OUTPUT_DIR}/${_test})
 endforeach()
diff --git a/Tests/Functional/Python/PyCore/histogram2d.py b/Tests/Functional/Python/PyCore/histogram2d.py
index 84c639fe60acd65fa1b97763f04045b5aac12b9d..c9df823fc106c1d7ae2f36d2d6c95e1131af61fb 100644
--- a/Tests/Functional/Python/PyCore/histogram2d.py
+++ b/Tests/Functional/Python/PyCore/histogram2d.py
@@ -1,6 +1,4 @@
 import numpy, os, sys, unittest
-
-sys.path.append("@CMAKE_LIBRARY_OUTPUT_DIRECTORY@")
 import bornagain as ba
 
 
diff --git a/Tests/Functional/Python/PyCore/intensitydata.py b/Tests/Functional/Python/PyCore/intensitydata.py
index bb6befd4d55f1ca540f9947c51480f6d899e7a58..e32620b5761b4456c2321d1fef2327618f7427ef 100644
--- a/Tests/Functional/Python/PyCore/intensitydata.py
+++ b/Tests/Functional/Python/PyCore/intensitydata.py
@@ -1,8 +1,6 @@
 # Functional test: tests of IntensityData object
 
 import numpy, os, sys, unittest
-
-sys.path.append("@CMAKE_LIBRARY_OUTPUT_DIRECTORY@")
 import bornagain as ba
 
 
diff --git a/Tests/Functional/Python/PyCore/intensitydata_io.py b/Tests/Functional/Python/PyCore/intensitydata_io.py
index f45b73ce23ba2440a47d18f09631a78a36d858c4..d11811068d2c1daa1e9d61dcc9abfda6fb6dba38 100644
--- a/Tests/Functional/Python/PyCore/intensitydata_io.py
+++ b/Tests/Functional/Python/PyCore/intensitydata_io.py
@@ -1,8 +1,6 @@
 # Functional test: tests of IO operations with the IntensityData object
 
 import math, numpy, os, sys, time, unittest
-
-sys.path.append("@CMAKE_LIBRARY_OUTPUT_DIRECTORY@")
 import bornagain as ba
 from bornagain import deg, deg2rad, rad2deg
 
diff --git a/Tests/Functional/Python/PyCore/intensitydata_io_tiff.py b/Tests/Functional/Python/PyCore/intensitydata_io_tiff.py
index 8ae1ad26bb573b3d8aa2cc250fc558bde0d78fa8..4b6e0f29ea8276f31dca5e720182f4d7b66aad09 100644
--- a/Tests/Functional/Python/PyCore/intensitydata_io_tiff.py
+++ b/Tests/Functional/Python/PyCore/intensitydata_io_tiff.py
@@ -1,8 +1,6 @@
 # Functional test: tests of IO operations with the IntensityData object
 
 import math, numpy, os, sys, time, unittest
-
-sys.path.append("@CMAKE_LIBRARY_OUTPUT_DIRECTORY@")
 import bornagain as ba
 
 
diff --git a/Tests/Functional/Python/PyCore/parameterpool.py b/Tests/Functional/Python/PyCore/parameterpool.py
index ca8d9d0252d1e16996bf01be127de90483996458..d61107d8fa06697b54d28a9354fbd02cf84b40db 100644
--- a/Tests/Functional/Python/PyCore/parameterpool.py
+++ b/Tests/Functional/Python/PyCore/parameterpool.py
@@ -1,6 +1,4 @@
 import numpy, os, sys, unittest
-
-sys.path.append("@CMAKE_LIBRARY_OUTPUT_DIRECTORY@")
 import bornagain as ba
 from bornagain import nm
 
diff --git a/Tests/Functional/Python/PyCore/samplebuilder.py b/Tests/Functional/Python/PyCore/samplebuilder.py
index 35dd5ab702679a4a61ab97a7fb1deafd29efbe3d..a8628a3f00297a5e1b1a6eb86271849546ba1c85 100644
--- a/Tests/Functional/Python/PyCore/samplebuilder.py
+++ b/Tests/Functional/Python/PyCore/samplebuilder.py
@@ -1,7 +1,5 @@
 import sys, unittest, ctypes
 import inspect
-
-sys.path.append("@CMAKE_LIBRARY_OUTPUT_DIRECTORY@")
 import bornagain as ba
 
 initial_width = 42
diff --git a/Tests/Functional/Python/PyCore/shape2d.py b/Tests/Functional/Python/PyCore/shape2d.py
index ad92c8211c84231b1c32cead1c399fa953c1e4cf..794763c866023c585212b9d63a17d86795ad9937 100644
--- a/Tests/Functional/Python/PyCore/shape2d.py
+++ b/Tests/Functional/Python/PyCore/shape2d.py
@@ -1,6 +1,4 @@
 import numpy, os, sys, unittest
-
-sys.path.append("@CMAKE_LIBRARY_OUTPUT_DIRECTORY@")
 import bornagain as ba
 
 
diff --git a/Tests/Functional/Python/PyCore/utils.py b/Tests/Functional/Python/PyCore/utils.py
index e17bba8b0543e96a73f88200e09b3111cdee0d36..ab4ff98bd51404693bee29f870da72e4fd25e783 100644
--- a/Tests/Functional/Python/PyCore/utils.py
+++ b/Tests/Functional/Python/PyCore/utils.py
@@ -3,8 +3,6 @@ Collection of utils for testing
 """
 
 import gzip, numpy, sys, os
-
-sys.path.append("@CMAKE_LIBRARY_OUTPUT_DIRECTORY@")
 import bornagain as ba
 from bornagain import deg, angstrom
 
@@ -42,39 +40,7 @@ def get_reference_data(filename):
         os.path.join(REFERENCE_DIR, filename))
 
 
-def get_simulation_MiniGISAS(sample=None):
-    simulation = ba.GISASSimulation()
-    simulation.setDetectorParameters(25, -2.0*deg, 2.0*deg, 25, 0.0*deg,
-                                     2.0*deg)
-    simulation.setBeamParameters(1.0*angstrom, 0.2*deg, 0.0*deg)
-    if sample:
-        simulation.setSample(sample)
-    return simulation
-
-
-def get_simulation_BasicGISAS(sample=None):
-    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)
-    if sample:
-        simulation.setSample(sample)
-    return simulation
-
-
-def plot_int(intensity):
-    import matplotlib, pylab
-    data = intensity.getArray() + 1
-    # data = numpy.abs(intensity.getArray())
-    phi_min = rad2deg(intensity.axis(0).lowerBound())
-    phi_max = rad2deg(intensity.axis(0).upperBound())
-    alpha_min = rad2deg(intensity.axis(1).lowerBound())
-    alpha_max = rad2deg(intensity.axis(1).upperBound())
-    im = pylab.imshow(data,
-                      norm=matplotlib.colors.LogNorm(),
-                      extent=[phi_min, phi_max, alpha_min, alpha_max])
-    cb = pylab.colorbar(im)
-    cb.set_label(r'Intensity (arb. u.)', size=16)
-    pylab.xlabel(r'$\phi_f (^{\circ})$', fontsize=16)
-    pylab.ylabel(r'$\alpha_f (^{\circ})$', fontsize=16)
-    pylab.show()
+def get_simulation_MiniGISAS(sample):
+    detector = ba.SphericalDetector(25, -2*deg, 2*deg, 25, 0*deg, 2*deg)
+    beam = ba.Beam(1., 1*angstrom, ba.Direction(0.2*deg, 0))
+    return ba.GISASSimulation(beam, sample, detector)
diff --git a/Tests/Functional/Python/PyExamples/CMakeLists.txt b/Tests/Functional/Python/PyExamples/CMakeLists.txt
index d7ca6675cf1ce5092063db38b89be494bd52a20d..4b3c7629b4306187be94bbc3534dfdd47f5ed557 100644
--- a/Tests/Functional/Python/PyExamples/CMakeLists.txt
+++ b/Tests/Functional/Python/PyExamples/CMakeLists.txt
@@ -9,17 +9,21 @@ set(output_dir ${TEST_OUTPUT_DIR_PY_EXAMPLES})
 file(MAKE_DIRECTORY ${output_dir})
 
 file(GLOB sim_examples ${PY_EXAMPLES_DIR}/sim*/*.py)
-file(GLOB fit_examples
-    ${PY_EXAMPLES_DIR}/fit55_Specular/FitSpecularBasics.py)
+if (WIN32)
+    # Convergence problem in Gauss-Kronrod integration
+    list(REMOVE_ITEM sim_examples ${PY_EXAMPLES_DIR}/sim03_Structures/Interference2DParaCrystal.py)
+endif()
+
+
+file(GLOB fit_examples ${PY_EXAMPLES_DIR}/fit55_Specular/FitSpecularBasics.py)
 set(examples ${sim_examples} ${fit_examples})
 
-set(test_script ${output_dir}/check_functionality.py)
-configure_file(${CMAKE_CURRENT_SOURCE_DIR}/check_functionality.in.py ${test_script} @ONLY)
+set(test_script ${TOOL_DIR}/code-tools/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} ${Python3_EXECUTABLE} ${test_script} ${example_path})
-    set_tests_properties(${test_name} PROPERTIES LABELS Fullcheck)
+    add_test(${test_name}
+        env PYTHONPATH=${CMAKE_LIBRARY_OUTPUT_DIRECTORY}
+        ${Python3_EXECUTABLE} ${test_script} -s ${example_path} ${output_dir})
 endforeach()
diff --git a/Tests/Functional/Python/PyExamples/check_functionality.in.py b/Tests/Functional/Python/PyExamples/check_functionality.in.py
deleted file mode 100644
index c885bd9275a6d62f4e8a408c6449a5a44854217d..0000000000000000000000000000000000000000
--- a/Tests/Functional/Python/PyExamples/check_functionality.in.py
+++ /dev/null
@@ -1,83 +0,0 @@
-"""
-Usage: python check_functionality.py <path-to-example>/<example>.py
-
-Checks functionality of BornAgain Python example by running it in 'embedded' way.
-
-The check passes successfully if the example runs without exceptions thrown and
-generates non-zero-size intensity image.
-"""
-
-import sys
-import os
-import matplotlib
-matplotlib.use('Agg')
-from matplotlib import pyplot as plt
-
-sys.path.append("@CMAKE_LIBRARY_OUTPUT_DIRECTORY@")
-
-output_dir = "@output_dir@"
-
-
-def get_figure(filename):
-    """
-    Returns pyplot figure of appropriate size
-    """
-    if "AllFormFactorsAvailable" in filename:
-        xsize, ysize = 1024, 768
-    else:
-        xsize, ysize = 640, 480
-
-    dpi = 72.
-    return plt.figure(figsize=(xsize/dpi, ysize/dpi))
-
-
-def exec_full(filepath):
-    """
-    Executes embedded python script.
-    http://stackoverflow.com/questions/436198/what-is-an-alternative-to-execfile-in-python-3
-    """
-    import os
-    global_namespace = {
-        "__file__": filepath,
-        "__name__": "__main__",
-        "__no_terminal__": True
-    }
-    sys.argv = []  # TODO: FIXME after cleanup in plot_utils
-    with open(filepath, 'rb') as file:
-        exec(compile(file.read(), filepath, 'exec'), global_namespace)
-
-
-def run_example(filename):
-    """
-    Tries to run python example and produce a *.png image
-    """
-    if not os.path.exists(filename):
-        raise Exception("Example script '" + filename + "' not found")
-
-    print("Input script: " + filename, flush=True)
-
-    fig = get_figure(filename)
-
-    exec_full(filename)
-    print("Input script completed.", flush=True)
-
-    plot_file_name = os.path.join(
-        output_dir,
-        os.path.splitext(os.path.basename(filename))[0] + ".png")
-    print("Output image: " + plot_file_name)
-    plt.savefig(plot_file_name, bbox_inches='tight')
-    plt.close(fig)
-
-    imgSize = os.path.getsize(plot_file_name)
-    if imgSize == 0:
-        raise Exception("Image file is empty")
-    if imgSize < 20000:
-        raise Exception("Image file is unplausibly small")
-
-
-if __name__ == '__main__':
-
-    if len(sys.argv) != 2:
-        exit("check_functionality called with wrong number of arguments")
-
-    run_example(sys.argv[1])
diff --git a/Tests/Functional/Python/PyFit/CMakeLists.txt b/Tests/Functional/Python/PyFit/CMakeLists.txt
index 81bed49f6c36ba321380085defc919750db1ee14..133fcae57fa9671cad5d9c457779428153694d90 100644
--- a/Tests/Functional/Python/PyFit/CMakeLists.txt
+++ b/Tests/Functional/Python/PyFit/CMakeLists.txt
@@ -5,18 +5,11 @@
 set(output_dir ${TEST_OUTPUT_DIR_PY_FIT})
 file(MAKE_DIRECTORY ${output_dir})
 
-file(GLOB sources RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} "*.py")
-foreach(_src ${sources})
-    configure_file(${CMAKE_CURRENT_SOURCE_DIR}/${_src}
-         ${output_dir}/${_src} @ONLY)
-endforeach()
-
-set(tests
-    "minimizer_api.py"
-    "standalone_fits.py"
-    "fitobjective_api.py"
-)
+file(GLOB tests RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} "*.py")
 
 foreach(_test ${tests})
-    add_test(PyFit.${_test} ${Python3_EXECUTABLE}  ${output_dir}/${_test})
+    configure_file(${CMAKE_CURRENT_SOURCE_DIR}/${_test} ${output_dir}/${_test} COPYONLY)
+    add_test(PyFit.${_test}
+        env PYTHONPATH=${CMAKE_LIBRARY_OUTPUT_DIRECTORY}
+        ${Python3_EXECUTABLE} ${output_dir}/${_test})
 endforeach()
diff --git a/Tests/Functional/Python/PyFit/README.md b/Tests/Functional/Python/PyFit/README.md
index 4221b4ca37838d42966e620b79ff078351be72bc..003cc60e8a4647cc10ded63f04bfe8ca8b62b7fb 100644
--- a/Tests/Functional/Python/PyFit/README.md
+++ b/Tests/Functional/Python/PyFit/README.md
@@ -1,4 +1,4 @@
-Collection of functional tests (Python)
+## Functional tests of libBornAgainFit
 
 Collection of fitting tests for libBornAgainFit library.
 Different geometries, number of fit parameters, variety
diff --git a/Tests/Functional/Python/PyFit/fitobjective_api.py b/Tests/Functional/Python/PyFit/fitobjective_api.py
index bdd91ddca2e7ec69241060340d955282abd73ecf..27ebf571de6cc73fdfd4f2e160efe1f58c15cf39 100644
--- a/Tests/Functional/Python/PyFit/fitobjective_api.py
+++ b/Tests/Functional/Python/PyFit/fitobjective_api.py
@@ -5,8 +5,6 @@ import sys
 import os
 import unittest
 import numpy as np
-
-sys.path.append("@CMAKE_LIBRARY_OUTPUT_DIRECTORY@")
 import bornagain as ba
 
 
diff --git a/Tests/Functional/Python/PyFit/minimizer_api.py b/Tests/Functional/Python/PyFit/minimizer_api.py
index 8eabce41be48ae779115c228a7f641801afda22e..ee5bbe4d0136bfb16bf6743e7823bd24317a4f6b 100644
--- a/Tests/Functional/Python/PyFit/minimizer_api.py
+++ b/Tests/Functional/Python/PyFit/minimizer_api.py
@@ -4,8 +4,6 @@ Testing python specific API for Minimizer related classes.
 import sys
 import os
 import unittest
-
-sys.path.append("@CMAKE_LIBRARY_OUTPUT_DIRECTORY@")
 import bornagain as ba
 
 
diff --git a/Tests/Functional/Python/PyFit/standalone_fits.py b/Tests/Functional/Python/PyFit/standalone_fits.py
index 53f8ff1b39c06af1a76c29189d69cd6ac411e301..e3518b757cac750dc8047eafa1643a1f736fccd8 100644
--- a/Tests/Functional/Python/PyFit/standalone_fits.py
+++ b/Tests/Functional/Python/PyFit/standalone_fits.py
@@ -4,8 +4,6 @@ Fitting scalar and residual based objective functions
 import sys
 import os
 import unittest
-
-sys.path.append("@CMAKE_LIBRARY_OUTPUT_DIRECTORY@")
 import bornagain as ba
 import numpy as np
 
diff --git a/Tests/Functional/Python/PyPersistence/CMakeLists.txt b/Tests/Functional/Python/PyPersistence/CMakeLists.txt
index c1960ad2c3ccd6d2911a54c5d0c6ccb4df7eef34..00a6defb321078447d52915d3704023abf0c6e5e 100644
--- a/Tests/Functional/Python/PyPersistence/CMakeLists.txt
+++ b/Tests/Functional/Python/PyPersistence/CMakeLists.txt
@@ -66,6 +66,6 @@ test_example(sim21_Reflectometry/BeamFullDivergence 2e-10)
 test_example(sim21_Reflectometry/TimeOfFlightReflectometry 2e-10)
 test_example(sim21_Reflectometry/TOFRWithResolution 2e-10)
 
-test_example(sim29_DepthProbe/DepthProbe 2e-10)
+test_example(special40_DepthProbe/DepthProbe 2e-10)
 
 test_example(sim31_Parameterization/SimulationParameters 2e-10)
diff --git a/Wrap/Python/__init__.py.in b/Wrap/Python/__init__.py.in
index 8f8981f68c44978a9ad9aeb6d45f71a0be7fb53b..4e140ee7f8c7afa03e057ed91f5e803c016525fe 100644
--- a/Wrap/Python/__init__.py.in
+++ b/Wrap/Python/__init__.py.in
@@ -26,3 +26,4 @@ from libBornAgainCore import *
 # To prevent inclusion of plotting tools during functional tests:
 if not "NOPLOT" in os.environ:
     from .plot_utils import *
+    import fit_monitor
diff --git a/Wrap/Python/fit_monitor.py b/Wrap/Python/fit_monitor.py
new file mode 100644
index 0000000000000000000000000000000000000000..3496dced8aee3257f470ace6f90663669ccef2bb
--- /dev/null
+++ b/Wrap/Python/fit_monitor.py
@@ -0,0 +1,258 @@
+#  **************************************************************************  #
+"""
+#   BornAgain: simulate and fit scattering at grazing incidence
+#
+#   @file      Wrap/Python/fit_monitor.py
+#   @brief     Plotter classes for monitoring fit progress.
+#
+#   @homepage  http://apps.jcns.fz-juelich.de/BornAgain
+#   @license   GNU General Public License v3 or higher (see COPYING)
+#   @copyright Forschungszentrum Juelich GmbH 2019
+#   @authors   Scientific Computing Group at MLZ (see CITATION, AUTHORS)
+"""
+#  **************************************************************************  #
+
+import bornagain as ba
+try:  # workaround for build servers
+    import numpy as np
+    from matplotlib import pyplot as plt
+    from matplotlib import gridspec, colors
+except Exception as e:
+    print("In plot_utils.py: {:s}".format(str(e)))
+
+label_fontsize = 16
+
+class Plotter:
+    """
+    Draws fit progress. Base class for simulation-specific classes (PlotterGISAS etc).
+    """
+    def __init__(self,
+                 zmin=None,
+                 zmax=None,
+                 xlabel=None,
+                 ylabel=None,
+                 units=ba.Axes.DEFAULT,
+                 aspect=None):
+
+        self._fig = plt.figure(figsize=(10.25, 7.69))
+        self._fig.canvas.draw()
+        self._zmin = zmin
+        self._zmax = zmax
+        self._xlabel = xlabel
+        self._ylabel = ylabel
+        self._units = units
+        self._aspect = aspect
+
+    def reset(self):
+        self._fig.clf()
+
+    def plot(self):
+        self._fig.tight_layout()
+        plt.pause(0.03)
+
+
+class PlotterGISAS(Plotter):
+    """
+    Draws fit progress, for GISAS simulation.
+    """
+    def __init__(self,
+                 zmin=None,
+                 zmax=None,
+                 xlabel=None,
+                 ylabel=None,
+                 units=ba.Axes.DEFAULT,
+                 aspect=None):
+        Plotter.__init__(self, zmin, zmax, xlabel, ylabel, units, aspect)
+
+    @staticmethod
+    def make_subplot(nplot):
+        plt.subplot(2, 2, nplot)
+        plt.subplots_adjust(wspace=0.2, hspace=0.2)
+
+    def plot(self, fit_objective):
+        Plotter.reset(self)
+
+        real_data = fit_objective.experimentalData()
+        sim_data = fit_objective.simulationResult()
+        diff = fit_objective.absoluteDifference()
+
+        self.make_subplot(1)
+
+        # same limits for both plots
+        arr = real_data.array()
+        zmax = np.amax(arr) if self._zmax is None else self._zmax
+        zmin = zmax*1e-6 if self._zmin is None else self._zmin
+
+        ba.plot_colormap(real_data,
+                         title="Experimental data",
+                         zmin=zmin,
+                         zmax=zmax,
+                         units=self._units,
+                         xlabel=self._xlabel,
+                         ylabel=self._ylabel,
+                         zlabel='',
+                         aspect=self._aspect)
+
+        self.make_subplot(2)
+        ba.plot_colormap(sim_data,
+                         title="Simulated data",
+                         zmin=zmin,
+                         zmax=zmax,
+                         units=self._units,
+                         xlabel=self._xlabel,
+                         ylabel=self._ylabel,
+                         zlabel='',
+                         aspect=self._aspect)
+
+        self.make_subplot(3)
+        ba.plot_colormap(diff,
+                         title="Difference",
+                         zmin=zmin,
+                         zmax=zmax,
+                         units=self._units,
+                         xlabel=self._xlabel,
+                         ylabel=self._ylabel,
+                         zlabel='',
+                         aspect=self._aspect)
+
+        self.make_subplot(4)
+        plt.title('Parameters')
+        plt.axis('off')
+
+        iteration_info = fit_objective.iterationInfo()
+
+        plt.text(
+            0.01, 0.85,
+            "Iterations  " + '{:d}'.format(iteration_info.iterationCount()))
+        plt.text(0.01, 0.75,
+                 "Chi2       " + '{:8.4f}'.format(iteration_info.chi2()))
+        index = 0
+        params = iteration_info.parameterMap()
+        for key in params:
+            plt.text(0.01, 0.55 - index*0.1,
+                     '{:30.30s}: {:6.3f}'.format(key, params[key]))
+            index = index + 1
+
+        Plotter.plot(self)
+
+
+class PlotterSpecular(Plotter):
+    """
+    Draws fit progress, for specular simulation.
+    """
+    def __init__(self, units=ba.Axes.DEFAULT):
+        Plotter.__init__(self)
+        self.gs = gridspec.GridSpec(1, 2, width_ratios=[2.5, 1], wspace=0)
+        self.units = units
+
+    def __call__(self, fit_objective):
+        self.plot(fit_objective)
+
+    @staticmethod
+    def as_si(val, ndp):
+        """
+        Fancy print of scientific-formatted values
+        :param val: numeric value
+        :param ndp: number of decimal digits to print
+        :return: a string corresponding to the _val_
+        """
+        s = '{x:0.{ndp:d}e}'.format(x=val, ndp=ndp)
+        m, e = s.split('e')
+        return r'{m:s}\times 10^{{{e:d}}}'.format(m=m, e=int(e))
+
+    @staticmethod
+    def trunc_str(token, length):
+        """
+        Truncates token if it is longer than length.
+
+        Example:
+            trunc_str("123456789", 8) returns "123456.."
+
+            trunc_str("123456789", 9) returns "123456789"
+
+        :param token: input string
+        :param length: max non-truncated length
+        :return:
+        """
+        return (token[:length - 2] + '..') if len(token) > length else token
+
+    def plot_table(self, fit_objective):
+
+        iteration_info = fit_objective.iterationInfo()
+
+        trunc_length = 9  # max string field width in the table
+        n_digits = 1  # number of decimal digits to print
+
+        n_iterations = iteration_info.iterationCount(
+        )  # current number of iterations passed
+        rel_dif = fit_objective.relativeDifference().array().max(
+        )  # maximum relative difference
+        fitted_parameters = iteration_info.parameterMap()
+
+        # creating table content
+        labels = ("Parameter", "Value")
+        table_data = [["Iteration", '${:d}$'.format(n_iterations)],
+                      [
+                          "$d_{r, max}$",
+                          '${:s}$'.format(self.as_si(rel_dif, n_digits))
+                      ]]
+        for key, value in fitted_parameters.iteritems():
+            table_data.append([
+                '{:s}'.format(self.trunc_str(key, trunc_length)),
+                '${:s}$'.format(self.as_si(value, n_digits))
+            ])
+
+        # creating table
+        axs = plt.subplot(self.gs[1])
+        axs.axis('tight')
+        axs.axis('off')
+        table = plt.table(cellText=table_data,
+                          colLabels=labels,
+                          cellLoc='center',
+                          loc='bottom left',
+                          bbox=[0.0, 0.0, 1.0, 1.0])
+
+    def plot_graph(self, fit_objective):
+        # retrieving data from fit suite
+        real_data = fit_objective.experimentalData()
+        sim_data = fit_objective.simulationResult()
+        unc_data = fit_objective.uncertaintyData()
+
+        # data values
+        sim_values = sim_data.array(self.units)
+        real_values = real_data.array(self.units)
+        unc_values = None if unc_data is None else unc_data.array(self.units)
+
+        # default font properties dictionary to use
+        font = {'family': 'serif', 'weight': 'normal', 'size': label_fontsize}
+
+        plt.subplot(self.gs[0])
+        plt.semilogy(sim_data.axis(), sim_values, 'b', real_data.axis(),
+                     real_values, 'k--')
+        if unc_values is not None:
+            plt.semilogy(real_data.axis(),
+                         real_values - unc_values,
+                         'xkcd:grey',
+                         alpha=0.6)
+            plt.semilogy(real_data.axis(),
+                         real_values + unc_values,
+                         'xkcd:grey',
+                         alpha=0.6)
+        plt.ylim((0.5*np.min(real_values), 5*np.max(real_values)))
+
+        xlabel = ba.plot_utils.get_axes_labels(real_data, self.units)[0]
+        legend = ['BornAgain', 'Data']
+        if unc_values is not None:
+            legend = ['BornAgain', 'Data', r'Data $\pm \sigma$']
+        plt.legend(legend, loc='upper right', prop=font)
+        plt.xlabel(xlabel, fontdict=font)
+        plt.ylabel("Intensity", fontdict=font)
+        plt.title("Specular data fitting", fontdict=font)
+
+    def plot(self, fit_objective):
+        Plotter.reset(self)
+
+        self.plot_graph(fit_objective)
+        self.plot_table(fit_objective)
+
+        Plotter.plot(self)
diff --git a/Wrap/Python/plot_utils.py b/Wrap/Python/plot_utils.py
index a19327dd31162e29775f5b03ccab21356e61d642..6b569d757f72ee5de8ad7934d57b9876cfc1a151 100644
--- a/Wrap/Python/plot_utils.py
+++ b/Wrap/Python/plot_utils.py
@@ -2,14 +2,13 @@
 """
 #   BornAgain: simulate and fit scattering at grazing incidence
 #
-#   @file      Wrap/Python.plot_utils
+#   @file      Wrap/Python/plot_utils.py
 #   @brief     Python extensions of the SWIG-generated Python module bornagain.
 #
 #   @homepage  http://apps.jcns.fz-juelich.de/BornAgain
 #   @license   GNU General Public License v3 or higher (see COPYING)
 #   @copyright Forschungszentrum Juelich GmbH 2016
-#   @authors   Scientific Computing Group at MLZ Garching
-#   @authors   J. Fisher, M. Ganeva, G. Pospelov, W. Van Herck, J. Wuttke
+#   @authors   Scientific Computing Group at MLZ (see CITATION, AUTHORS)
 """
 #  **************************************************************************  #
 
@@ -48,19 +47,19 @@ def translate_axis_label(label):
     """
     label_dict = {
         'X [nbins]': r'$X \; $(bins)',
+        'X [mm]': r'$X \; $(mm)',
+        'Y [nbins]': r'$Y \; $(bins)',
+        'Y [mm]': r'$Y \; $(mm)',
         'phi_f [rad]': r'$\varphi_f \; $(rad)',
         'phi_f [deg]': r'$\varphi_f \; $(deg)',
         'alpha_i [rad]': r'$\alpha_i \; $(rad)',
         'alpha_i [deg]': r'$\alpha_i \; $(deg)',
-        'X [mm]': r'$X \; $(mm)',
-        'Qx [1/nm]': r'$Q_x \; $(nm$^{-1}$)',
-        'Qy [1/nm]': r'$Q_y \; $(nm$^{-1}$)',
-        'Q [1/nm]': r'$Q \; $(nm$^{-1}$)',
-        'Y [nbins]': r'$Y \; $(bins)',
         'alpha_f [rad]': r'$\alpha_f \; $(rad)',
         'alpha_f [deg]': r'$\alpha_f \; $(deg)',
-        'Y [mm]': r'$Y \; $(mm)',
+        'Qx [1/nm]': r'$Q_x \; $(nm$^{-1}$)',
+        'Qy [1/nm]': r'$Q_y \; $(nm$^{-1}$)',
         'Qz [1/nm]': r'$Q_z \; $(nm$^{-1}$)',
+        'Q [1/nm]': r'$Q \; $(nm$^{-1}$)',
         'Position [nm]': r'$Position \; $(nm)'
     }
     if label in label_dict.keys():
@@ -85,38 +84,32 @@ def get_axes_labels(result, units):
     return labels
 
 
-def plot_array(array,
-               zmin=None,
-               zmax=None,
-               xlabel=None,
-               ylabel=None,
-               zlabel=None,
-               title=None,
-               axes_limits=None,
-               **kwargs):
+def plot_array(array, axes_limits=None, **kwargs):
     """
     Plots numpy array as a colormap in log scale.
     """
 
-    zlabel = "Intensity" if zlabel is None else zlabel
-
-    zmax = np.amax(array) if zmax is None else zmax
-    zmin = 1e-6*zmax if zmin is None else zmin
+    zmax = kwargs.pop('intensity_max', np.amax(array))
+    zmin = kwargs.pop('intensity_min', 1e-6*zmax)
 
     if zmin == zmax == 0.0:
         norm = colors.Normalize(0, 1)
     else:
         norm = colors.LogNorm(zmin, zmax)
 
+    xlabel = kwargs.pop('xlabel', None)
+    ylabel = kwargs.pop('ylabel', None)
+    zlabel = kwargs.pop('zlabel', "Intensity")
+
+    title = kwargs.pop('title', None)
+
     im = plt.imshow(array, norm=norm, extent=axes_limits, **kwargs)
     cb = plt.colorbar(im, pad=0.025)
 
     if xlabel:
         plt.xlabel(xlabel, fontsize=label_fontsize)
-
     if ylabel:
         plt.ylabel(ylabel, fontsize=label_fontsize)
-
     if zlabel:
         cb.set_label(zlabel, size=label_fontsize)
 
@@ -124,14 +117,7 @@ def plot_array(array,
         plt.title(title)
 
 
-def plot_histogram(hist,
-                   zmin=None,
-                   zmax=None,
-                   xlabel=None,
-                   ylabel=None,
-                   zlabel=None,
-                   title=None,
-                   **kwargs):
+def plot_histogram(hist, **kwargs):
     """
     Plots intensity data as color map
     :param intensity: Histogram2D object obtained from GISASSimulation
@@ -139,11 +125,10 @@ def plot_histogram(hist,
     :param zmax: Max value on amplitude's color bar
     """
 
-    if not xlabel:
-        xlabel = translate_axis_label(hist.xAxis().getName())
-
-    if not ylabel:
-        ylabel = translate_axis_label(hist.yAxis().getName())
+    if not 'xlabel' in kwargs:
+        kwargs['xlabel'] = translate_axis_label(hist.xAxis().getName())
+    if not 'ylabel' in kwargs:
+        kwargs['ylabel'] = translate_axis_label(hist.yAxis().getName())
 
     axes_limits = [
         hist.getXmin(),
@@ -152,26 +137,12 @@ def plot_histogram(hist,
         hist.getYmax()
     ]
 
-    plot_array(hist.array(),
-               zmin=zmin,
-               zmax=zmax,
-               xlabel=xlabel,
-               ylabel=ylabel,
-               zlabel=zlabel,
-               title=title,
-               axes_limits=axes_limits,
-               **kwargs)
-
-
-def plot_colormap(result,
-                  zmin=None,
-                  zmax=None,
-                  units=ba.Axes.DEFAULT,
-                  xlabel=None,
-                  ylabel=None,
-                  zlabel=None,
-                  title=None,
-                  **kwargs):
+    title = kwargs.pop('title', None)
+
+    plot_array(hist.array(), title=title, axes_limits=axes_limits, **kwargs)
+
+
+def plot_colormap(result, **kwargs):
     """
     Plots intensity data as color map
     :param result: SimulationResult from GISAS/OffSpecSimulation
@@ -179,30 +150,19 @@ def plot_colormap(result,
     :param zmax: Max value on amplitude's color bar
     """
 
+    units = kwargs.pop('units', ba.Axes.DEFAULT)
     axes_limits = get_axes_limits(result, units)
     axes_labels = get_axes_labels(result, units)
-    xlabel = axes_labels[0] if xlabel is None else xlabel
-    ylabel = axes_labels[1] if ylabel is None else ylabel
-
-    plot_array(result.array(),
-               zmin=zmin,
-               zmax=zmax,
-               xlabel=xlabel,
-               ylabel=ylabel,
-               zlabel=zlabel,
-               title=title,
-               axes_limits=axes_limits,
-               **kwargs)
-
-
-def plot_specular_simulation_result(result,
-                                    ymin=None,
-                                    ymax=None,
-                                    units=ba.Axes.DEFAULT,
-                                    xlabel=None,
-                                    ylabel=None,
-                                    title=None,
-                                    **kwargs):
+
+    if not 'xlabel' in kwargs:
+        kwargs['xlabel'] = axes_labels[0]
+    if not 'ylabel' in kwargs:
+        kwargs['ylabel'] = axes_labels[1]
+
+    plot_array(result.array(), axes_limits=axes_limits, **kwargs)
+
+
+def plot_specular_simulation_result(result, **kwargs):
     """
     Plots intensity data for specular simulation result
     :param result: SimulationResult from SpecularSimulation
@@ -211,13 +171,16 @@ def plot_specular_simulation_result(result,
     :param units: units on the x-axis
     """
 
+    units = kwargs.pop('units', ba.Axes.DEFAULT)
     intensity = result.array(units)
     x_axis = result.axis(units)
-    ymax = np.amax(intensity)*2.0 if ymax is None else ymax
-    ymin = max(np.amin(intensity)*0.5, 1e-18*ymax) if ymin is None else ymin
 
-    xlabel = get_axes_labels(result, units)[0] if xlabel is None else xlabel
-    ylabel = "Intensity" if ylabel is None else ylabel
+    ymax = kwargs.pop('intensity_max', np.amax(np.amax(intensity)*2))
+    ymin = kwargs.pop('intensity_min', max(np.amin(intensity)*0.5, 1e-18*ymax))
+
+    xlabel = kwargs.pop('xlabel', get_axes_labels(result, units)[0])
+    ylabel = kwargs.pop('ylabel', "Intensity")
+    title = kwargs.pop('title', None)
 
     plt.semilogy(x_axis, intensity, **kwargs)
 
@@ -225,7 +188,6 @@ def plot_specular_simulation_result(result,
 
     if xlabel:
         plt.xlabel(xlabel, fontsize=label_fontsize)
-
     if ylabel:
         plt.ylabel(ylabel, fontsize=label_fontsize)
 
@@ -233,15 +195,7 @@ def plot_specular_simulation_result(result,
         plt.title(title)
 
 
-def plot_simulation_result(result,
-                           intensity_min=None,
-                           intensity_max=None,
-                           units=ba.Axes.DEFAULT,
-                           xlabel=None,
-                           ylabel=None,
-                           postpone_show=False,
-                           title=None,
-                           **kwargs):
+def plot_simulation_result(result, **kwargs):
     """
     Draws simulation result and (optionally) shows the plot.
     :param result_: SimulationResult object obtained from GISAS/OffSpec/SpecularSimulation
@@ -250,249 +204,12 @@ def plot_simulation_result(result,
     :param units: units for plot axes
     :param postpone_show: postpone showing the plot for later tuning (False by default)
     """
+    postpone_show = kwargs.pop('postpone_show', False)
+
     if len(result.array().shape) == 1:  # 1D data, specular simulation assumed
-        plot_specular_simulation_result(result, intensity_min, intensity_max,
-                                        units, **kwargs)
+        plot_specular_simulation_result(result, **kwargs)
     else:
-        plot_colormap(result,
-                      intensity_min,
-                      intensity_max,
-                      units,
-                      xlabel,
-                      ylabel,
-                      title=title,
-                      **kwargs)
+        plot_colormap(result, **kwargs)
     plt.tight_layout()
-    if not postpone_show:
+    if not (postpone_show):
         plt.show()
-
-
-class Plotter:
-    def __init__(self,
-                 zmin=None,
-                 zmax=None,
-                 xlabel=None,
-                 ylabel=None,
-                 units=ba.Axes.DEFAULT,
-                 aspect=None):
-
-        self._fig = plt.figure(figsize=(10.25, 7.69))
-        self._fig.canvas.draw()
-        self._zmin = zmin
-        self._zmax = zmax
-        self._xlabel = xlabel
-        self._ylabel = ylabel
-        self._units = units
-        self._aspect = aspect
-
-    def reset(self):
-        self._fig.clf()
-
-    def plot(self):
-        self._fig.tight_layout()
-        plt.pause(0.03)
-
-
-class PlotterGISAS(Plotter):
-    def __init__(self,
-                 zmin=None,
-                 zmax=None,
-                 xlabel=None,
-                 ylabel=None,
-                 units=ba.Axes.DEFAULT,
-                 aspect=None):
-        Plotter.__init__(self, zmin, zmax, xlabel, ylabel, units, aspect)
-
-    @staticmethod
-    def make_subplot(nplot):
-        plt.subplot(2, 2, nplot)
-        plt.subplots_adjust(wspace=0.2, hspace=0.2)
-
-    def plot(self, fit_objective):
-        Plotter.reset(self)
-
-        real_data = fit_objective.experimentalData()
-        sim_data = fit_objective.simulationResult()
-        diff = fit_objective.absoluteDifference()
-
-        self.make_subplot(1)
-
-        # same limits for both plots
-        arr = real_data.array()
-        zmax = np.amax(arr) if self._zmax is None else self._zmax
-        zmin = zmax*1e-6 if self._zmin is None else self._zmin
-
-        ba.plot_colormap(real_data,
-                         title="Experimental data",
-                         zmin=zmin,
-                         zmax=zmax,
-                         units=self._units,
-                         xlabel=self._xlabel,
-                         ylabel=self._ylabel,
-                         zlabel='',
-                         aspect=self._aspect)
-
-        self.make_subplot(2)
-        ba.plot_colormap(sim_data,
-                         title="Simulated data",
-                         zmin=zmin,
-                         zmax=zmax,
-                         units=self._units,
-                         xlabel=self._xlabel,
-                         ylabel=self._ylabel,
-                         zlabel='',
-                         aspect=self._aspect)
-
-        self.make_subplot(3)
-        ba.plot_colormap(diff,
-                         title="Difference",
-                         zmin=zmin,
-                         zmax=zmax,
-                         units=self._units,
-                         xlabel=self._xlabel,
-                         ylabel=self._ylabel,
-                         zlabel='',
-                         aspect=self._aspect)
-
-        self.make_subplot(4)
-        plt.title('Parameters')
-        plt.axis('off')
-
-        iteration_info = fit_objective.iterationInfo()
-
-        plt.text(
-            0.01, 0.85,
-            "Iterations  " + '{:d}'.format(iteration_info.iterationCount()))
-        plt.text(0.01, 0.75,
-                 "Chi2       " + '{:8.4f}'.format(iteration_info.chi2()))
-        index = 0
-        params = iteration_info.parameterMap()
-        for key in params:
-            plt.text(0.01, 0.55 - index*0.1,
-                     '{:30.30s}: {:6.3f}'.format(key, params[key]))
-            index = index + 1
-
-        Plotter.plot(self)
-
-
-class PlotterSpecular(Plotter):
-    """
-    Draws fit progress. Intended specifically for observing
-    specular data fit.
-    """
-    def __init__(self, units=ba.Axes.DEFAULT):
-        Plotter.__init__(self)
-        self.gs = gridspec.GridSpec(1, 2, width_ratios=[2.5, 1], wspace=0)
-        self.units = units
-
-    def __call__(self, fit_objective):
-        self.plot(fit_objective)
-
-    @staticmethod
-    def as_si(val, ndp):
-        """
-        Fancy print of scientific-formatted values
-        :param val: numeric value
-        :param ndp: number of decimal digits to print
-        :return: a string corresponding to the _val_
-        """
-        s = '{x:0.{ndp:d}e}'.format(x=val, ndp=ndp)
-        m, e = s.split('e')
-        return r'{m:s}\times 10^{{{e:d}}}'.format(m=m, e=int(e))
-
-    @staticmethod
-    def trunc_str(token, length):
-        """
-        Truncates token if it is longer than length.
-
-        Example:
-            trunc_str("123456789", 8) returns "123456.."
-
-            trunc_str("123456789", 9) returns "123456789"
-
-        :param token: input string
-        :param length: max non-truncated length
-        :return:
-        """
-        return (token[:length - 2] + '..') if len(token) > length else token
-
-    def plot_table(self, fit_objective):
-
-        iteration_info = fit_objective.iterationInfo()
-
-        trunc_length = 9  # max string field width in the table
-        n_digits = 1  # number of decimal digits to print
-
-        n_iterations = iteration_info.iterationCount(
-        )  # current number of iterations passed
-        rel_dif = fit_objective.relativeDifference().array().max(
-        )  # maximum relative difference
-        fitted_parameters = iteration_info.parameterMap()
-
-        # creating table content
-        labels = ("Parameter", "Value")
-        table_data = [["Iteration", '${:d}$'.format(n_iterations)],
-                      [
-                          "$d_{r, max}$",
-                          '${:s}$'.format(self.as_si(rel_dif, n_digits))
-                      ]]
-        for key, value in fitted_parameters.iteritems():
-            table_data.append([
-                '{:s}'.format(self.trunc_str(key, trunc_length)),
-                '${:s}$'.format(self.as_si(value, n_digits))
-            ])
-
-        # creating table
-        axs = plt.subplot(self.gs[1])
-        axs.axis('tight')
-        axs.axis('off')
-        table = plt.table(cellText=table_data,
-                          colLabels=labels,
-                          cellLoc='center',
-                          loc='bottom left',
-                          bbox=[0.0, 0.0, 1.0, 1.0])
-
-    def plot_graph(self, fit_objective):
-        # retrieving data from fit suite
-        real_data = fit_objective.experimentalData()
-        sim_data = fit_objective.simulationResult()
-        unc_data = fit_objective.uncertaintyData()
-
-        # data values
-        sim_values = sim_data.array(self.units)
-        real_values = real_data.array(self.units)
-        unc_values = None if unc_data is None else unc_data.array(self.units)
-
-        # default font properties dictionary to use
-        font = {'family': 'serif', 'weight': 'normal', 'size': label_fontsize}
-
-        plt.subplot(self.gs[0])
-        plt.semilogy(sim_data.axis(), sim_values, 'b', real_data.axis(),
-                     real_values, 'k--')
-        if unc_values is not None:
-            plt.semilogy(real_data.axis(),
-                         real_values - unc_values,
-                         'xkcd:grey',
-                         alpha=0.6)
-            plt.semilogy(real_data.axis(),
-                         real_values + unc_values,
-                         'xkcd:grey',
-                         alpha=0.6)
-        plt.ylim((0.5*np.min(real_values), 5*np.max(real_values)))
-
-        xlabel = get_axes_labels(real_data, self.units)[0]
-        legend = ['BornAgain', 'Data']
-        if unc_values is not None:
-            legend = ['BornAgain', 'Data', r'Data $\pm \sigma$']
-        plt.legend(legend, loc='upper right', prop=font)
-        plt.xlabel(xlabel, fontdict=font)
-        plt.ylabel("Intensity", fontdict=font)
-        plt.title("Specular data fitting", fontdict=font)
-
-    def plot(self, fit_objective):
-        Plotter.reset(self)
-
-        self.plot_graph(fit_objective)
-        self.plot_table(fit_objective)
-
-        Plotter.plot(self)
diff --git a/Wrap/Swig/libBornAgainCore.i b/Wrap/Swig/libBornAgainCore.i
index e3c40850623c75f4a0758465c2c0df1a6036f1a1..69114fd852ee5a329f8089a88f44df97c4548639 100644
--- a/Wrap/Swig/libBornAgainCore.i
+++ b/Wrap/Swig/libBornAgainCore.i
@@ -26,6 +26,8 @@
 
 %ignore ISpecularScan;
 
+%import(module="libBornAgainFit") ""
+
 %rename(setSampleBuilderCpp) ISimulation::setSampleBuilder;
 %rename(setSampleBuilderCpp) SpecularSimulation::setSampleBuilder;
 %rename(addSimulationAndData_cpp) FitObjective::addSimulationAndData;
diff --git a/auto/Wrap/doxygenBase.i b/auto/Wrap/doxygenBase.i
index a28a5913a953f0be27a10e9cce37c1cd665433e3..b9bb4d6f0f66217816bb62387d10017d5bda3046 100644
--- a/auto/Wrap/doxygenBase.i
+++ b/auto/Wrap/doxygenBase.i
@@ -550,6 +550,11 @@ Returns value of last point of axis.
 Returns distance from first to last point. 
 ";
 
+%feature("docstring")  IAxis::center "double IAxis::center() const
+
+Returns midpoint of axis. 
+";
+
 %feature("docstring")  IAxis::binCenter "virtual double IAxis::binCenter(size_t index) const =0
 ";
 
diff --git a/auto/Wrap/doxygenDevice.i b/auto/Wrap/doxygenDevice.i
index 4d751daeffaf49b8043aa92fd57989ff7c7316f2..c1a20b8f1ca65a87e951a4fdea27fb9a59aaf111 100644
--- a/auto/Wrap/doxygenDevice.i
+++ b/auto/Wrap/doxygenDevice.i
@@ -2850,6 +2850,26 @@ alpha_max:
 upper edge of last alpha-bin 
 ";
 
+%feature("docstring")  SphericalDetector::SphericalDetector "SphericalDetector::SphericalDetector(size_t n_bin, double width, double phi, double alpha)
+
+Spherical detector constructor with quadratic angle ranges
+
+Parameters:
+-----------
+
+n_bin: 
+number of bins per direction
+
+width: 
+angular width
+
+phi: 
+central phi angle
+
+alpha: 
+central alpha angle 
+";
+
 %feature("docstring")  SphericalDetector::SphericalDetector "SphericalDetector::SphericalDetector(const SphericalDetector &other)
 ";
 
@@ -3003,16 +3023,16 @@ Returns true if area defined by two bins is inside or on border of polygon (more
 // File: classConvolve_1_1Workspace.xml
 
 
-// File: namespace_0d105.xml
+// File: namespace_0d107.xml
 
 
-// File: namespace_0d33.xml
+// File: namespace_0d35.xml
 
 
-// File: namespace_0d56.xml
+// File: namespace_0d58.xml
 
 
-// File: namespace_0d62.xml
+// File: namespace_0d64.xml
 
 
 // File: namespaceArrayUtils.xml
@@ -3126,6 +3146,11 @@ Parse double values from string to vector of double.
 ";
 
 
+// File: namespaceDetectorUtils.xml
+%feature("docstring")  DetectorUtils::isQuadratic "bool DetectorUtils::isQuadratic(const IDetector2D &det)
+";
+
+
 // File: namespaceIntensityDataFunctions.xml
 %feature("docstring")  IntensityDataFunctions::RelativeDifference "double IntensityDataFunctions::RelativeDifference(const SimulationResult &dat, const SimulationResult &ref)
 
@@ -3287,6 +3312,12 @@ make Swappable
 // File: DetectorMask_8h.xml
 
 
+// File: DetectorUtils_8cpp.xml
+
+
+// File: DetectorUtils_8h.xml
+
+
 // File: IDetector_8cpp.xml
 
 
diff --git a/auto/Wrap/libBornAgainBase.py b/auto/Wrap/libBornAgainBase.py
index 4248754550195b1ef42c9afee34922900443535d..42bc36a22441656ddad66f52c1772fc53170097f 100644
--- a/auto/Wrap/libBornAgainBase.py
+++ b/auto/Wrap/libBornAgainBase.py
@@ -2109,6 +2109,16 @@ class IAxis(object):
         """
         return _libBornAgainBase.IAxis_span(self)
 
+    def center(self):
+        r"""
+        center(IAxis self) -> double
+        double IAxis::center() const
+
+        Returns midpoint of axis. 
+
+        """
+        return _libBornAgainBase.IAxis_center(self)
+
     def binCenter(self, index):
         r"""
         binCenter(IAxis self, size_t index) -> double
diff --git a/auto/Wrap/libBornAgainBase_wrap.cpp b/auto/Wrap/libBornAgainBase_wrap.cpp
index 20b8fd3b0156933b34560c4c169551c396535577..c0aa71b944b19ac95bc887ea291f25d177a4c9c5 100644
--- a/auto/Wrap/libBornAgainBase_wrap.cpp
+++ b/auto/Wrap/libBornAgainBase_wrap.cpp
@@ -26004,6 +26004,29 @@ fail:
 }
 
 
+SWIGINTERN PyObject *_wrap_IAxis_center(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+  PyObject *resultobj = 0;
+  IAxis *arg1 = (IAxis *) 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  PyObject *swig_obj[1] ;
+  double result;
+  
+  if (!args) SWIG_fail;
+  swig_obj[0] = args;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_IAxis, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IAxis_center" "', argument " "1"" of type '" "IAxis const *""'"); 
+  }
+  arg1 = reinterpret_cast< IAxis * >(argp1);
+  result = (double)((IAxis const *)arg1)->center();
+  resultobj = SWIG_From_double(static_cast< double >(result));
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
 SWIGINTERN PyObject *_wrap_IAxis_binCenter(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
   IAxis *arg1 = (IAxis *) 0 ;
@@ -33303,6 +33326,13 @@ static PyMethodDef SwigMethods[] = {
 		"Returns distance from first to last point. \n"
 		"\n"
 		""},
+	 { "IAxis_center", _wrap_IAxis_center, METH_O, "\n"
+		"IAxis_center(IAxis self) -> double\n"
+		"double IAxis::center() const\n"
+		"\n"
+		"Returns midpoint of axis. \n"
+		"\n"
+		""},
 	 { "IAxis_binCenter", _wrap_IAxis_binCenter, METH_VARARGS, "\n"
 		"IAxis_binCenter(IAxis self, size_t index) -> double\n"
 		"virtual double IAxis::binCenter(size_t index) const =0\n"
diff --git a/auto/Wrap/libBornAgainCore.py b/auto/Wrap/libBornAgainCore.py
index f2636b40ee8b12804d84bbc1f946ff217738f66d..b29772e589413c6221be05902c337aa03a72b949 100644
--- a/auto/Wrap/libBornAgainCore.py
+++ b/auto/Wrap/libBornAgainCore.py
@@ -1678,6 +1678,7 @@ class vector_pvacuum_double_t(object):
 # Register vector_pvacuum_double_t in _libBornAgainCore:
 _libBornAgainCore.vector_pvacuum_double_t_swigregister(vector_pvacuum_double_t)
 
+import libBornAgainFit
 import libBornAgainBase
 class kvector_t(object):
     r"""Proxy of C++ BasicVector3D< double > class."""
diff --git a/auto/Wrap/libBornAgainDevice.py b/auto/Wrap/libBornAgainDevice.py
index 476eb2340edadc07de5eab51045a5a0a2744400e..15bb81134b2bff3b13dce40d6c6ce168aea6fa74 100644
--- a/auto/Wrap/libBornAgainDevice.py
+++ b/auto/Wrap/libBornAgainDevice.py
@@ -4919,6 +4919,7 @@ class SphericalDetector(IDetector2D):
         r"""
         __init__(SphericalDetector self) -> SphericalDetector
         __init__(SphericalDetector self, size_t n_phi, double phi_min, double phi_max, size_t n_alpha, double alpha_min, double alpha_max) -> SphericalDetector
+        __init__(SphericalDetector self, size_t n_bin, double width, double phi, double alpha) -> SphericalDetector
         __init__(SphericalDetector self, SphericalDetector other) -> SphericalDetector
         SphericalDetector::SphericalDetector(const SphericalDetector &other)
 
diff --git a/auto/Wrap/libBornAgainDevice_wrap.cpp b/auto/Wrap/libBornAgainDevice_wrap.cpp
index 4a781d0bf52a2d09bc608c072e06df875479dd73..bfa18e94307cadf43cddd59109f8a8298dcaa4ac 100644
--- a/auto/Wrap/libBornAgainDevice_wrap.cpp
+++ b/auto/Wrap/libBornAgainDevice_wrap.cpp
@@ -39530,6 +39530,51 @@ fail:
 
 
 SWIGINTERN PyObject *_wrap_new_SphericalDetector__SWIG_2(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
+  PyObject *resultobj = 0;
+  size_t arg1 ;
+  double arg2 ;
+  double arg3 ;
+  double arg4 ;
+  size_t val1 ;
+  int ecode1 = 0 ;
+  double val2 ;
+  int ecode2 = 0 ;
+  double val3 ;
+  int ecode3 = 0 ;
+  double val4 ;
+  int ecode4 = 0 ;
+  SphericalDetector *result = 0 ;
+  
+  if ((nobjs < 4) || (nobjs > 4)) SWIG_fail;
+  ecode1 = SWIG_AsVal_size_t(swig_obj[0], &val1);
+  if (!SWIG_IsOK(ecode1)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode1), "in method '" "new_SphericalDetector" "', argument " "1"" of type '" "size_t""'");
+  } 
+  arg1 = static_cast< size_t >(val1);
+  ecode2 = SWIG_AsVal_double(swig_obj[1], &val2);
+  if (!SWIG_IsOK(ecode2)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "new_SphericalDetector" "', argument " "2"" of type '" "double""'");
+  } 
+  arg2 = static_cast< double >(val2);
+  ecode3 = SWIG_AsVal_double(swig_obj[2], &val3);
+  if (!SWIG_IsOK(ecode3)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode3), "in method '" "new_SphericalDetector" "', argument " "3"" of type '" "double""'");
+  } 
+  arg3 = static_cast< double >(val3);
+  ecode4 = SWIG_AsVal_double(swig_obj[3], &val4);
+  if (!SWIG_IsOK(ecode4)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode4), "in method '" "new_SphericalDetector" "', argument " "4"" of type '" "double""'");
+  } 
+  arg4 = static_cast< double >(val4);
+  result = (SphericalDetector *)new SphericalDetector(arg1,arg2,arg3,arg4);
+  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_SphericalDetector, SWIG_POINTER_NEW |  0 );
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
+SWIGINTERN PyObject *_wrap_new_SphericalDetector__SWIG_3(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
   SphericalDetector *arg1 = 0 ;
   void *argp1 = 0 ;
@@ -39569,7 +39614,35 @@ SWIGINTERN PyObject *_wrap_new_SphericalDetector(PyObject *self, PyObject *args)
     int res = SWIG_ConvertPtr(argv[0], 0, SWIGTYPE_p_SphericalDetector, SWIG_POINTER_NO_NULL | 0);
     _v = SWIG_CheckState(res);
     if (_v) {
-      return _wrap_new_SphericalDetector__SWIG_2(self, argc, argv);
+      return _wrap_new_SphericalDetector__SWIG_3(self, argc, argv);
+    }
+  }
+  if (argc == 4) {
+    int _v;
+    {
+      int res = SWIG_AsVal_size_t(argv[0], NULL);
+      _v = SWIG_CheckState(res);
+    }
+    if (_v) {
+      {
+        int res = SWIG_AsVal_double(argv[1], NULL);
+        _v = SWIG_CheckState(res);
+      }
+      if (_v) {
+        {
+          int res = SWIG_AsVal_double(argv[2], NULL);
+          _v = SWIG_CheckState(res);
+        }
+        if (_v) {
+          {
+            int res = SWIG_AsVal_double(argv[3], NULL);
+            _v = SWIG_CheckState(res);
+          }
+          if (_v) {
+            return _wrap_new_SphericalDetector__SWIG_2(self, argc, argv);
+          }
+        }
+      }
     }
   }
   if (argc == 6) {
@@ -39618,6 +39691,7 @@ fail:
     "  Possible C/C++ prototypes are:\n"
     "    SphericalDetector::SphericalDetector()\n"
     "    SphericalDetector::SphericalDetector(size_t,double,double,size_t,double,double)\n"
+    "    SphericalDetector::SphericalDetector(size_t,double,double,double)\n"
     "    SphericalDetector::SphericalDetector(SphericalDetector const &)\n");
   return 0;
 }
@@ -47082,6 +47156,7 @@ static PyMethodDef SwigMethods[] = {
 	 { "new_SphericalDetector", _wrap_new_SphericalDetector, METH_VARARGS, "\n"
 		"SphericalDetector()\n"
 		"SphericalDetector(size_t n_phi, double phi_min, double phi_max, size_t n_alpha, double alpha_min, double alpha_max)\n"
+		"SphericalDetector(size_t n_bin, double width, double phi, double alpha)\n"
 		"new_SphericalDetector(SphericalDetector other) -> SphericalDetector\n"
 		"SphericalDetector::SphericalDetector(const SphericalDetector &other)\n"
 		"\n"
diff --git a/cmake/BornAgain/Directories.cmake b/cmake/BornAgain/Directories.cmake
index bc51aee1ac7b5e6402f5d24a561f18837b737bec..03eb4f793603304ead0e2e6b71d7e27107af2b26 100644
--- a/cmake/BornAgain/Directories.cmake
+++ b/cmake/BornAgain/Directories.cmake
@@ -2,6 +2,7 @@
 # source directories
 # -----------------------------------------------------------------------------
 
+set(TOOL_DIR ${CMAKE_SOURCE_DIR}/devtools)
 set(WRAP_DIR ${CMAKE_SOURCE_DIR}/Wrap)
 set(SWIG_DIR ${WRAP_DIR}/Swig)
 set(PY_EXAMPLES_DIR ${CMAKE_SOURCE_DIR}/Examples/Python)
diff --git a/cmake/BornAgain/LineLength.cmake b/cmake/BornAgain/LineLength.cmake
index 311637cedbd93d009f65c27dcf393e8ff2228970..c5025e90e94ab60ca713668e509dce0b285c89e2 100644
--- a/cmake/BornAgain/LineLength.cmake
+++ b/cmake/BornAgain/LineLength.cmake
@@ -4,7 +4,7 @@ set(WEB_LEN_LIM 85) # maximum line length of code for display in web docs
 if(NOT MSVC)
 
     set(LINECOUNT
-        ${Python3_EXECUTABLE} ${CMAKE_SOURCE_DIR}/devtools/line-count/check-line-length.py)
+        ${Python3_EXECUTABLE} ${TOOL_DIR}/line-count/check-line-length.py)
 
     foreach(dir ${AllComponents})
         file(GLOB_RECURSE src1 ${dir}/*.cpp)
diff --git a/cmake/BornAgain/PythonAPI.cmake b/cmake/BornAgain/PythonAPI.cmake
index fd9c0bb25900a15f978e87b43345aaa69f04effb..fd1fb7d8a6e9acde0a9cc8323603265e030f6e55 100644
--- a/cmake/BornAgain/PythonAPI.cmake
+++ b/cmake/BornAgain/PythonAPI.cmake
@@ -4,9 +4,6 @@ if(NOT BORNAGAIN_PYTHON)
     message(FATAL_ERROR PythonAPI included though BORNAGAIN_PYTHON=false)
 endif()
 
-configure_file(${WRAP_DIR}/Python/plot_utils.py
-    ${CMAKE_BINARY_DIR}/lib/bornagain/plot_utils.py COPYONLY)
-
 if(WIN32)
     set(BA_MODULES_IMPORT_PATH ../../bin)
 else()
@@ -16,8 +13,12 @@ if(BORNAGAIN_APPLE_BUNDLE)
     set(BA_MODULES_IMPORT_PATH
         lib/BornAgain-${BornAgain_VERSION_MAJOR}.${BornAgain_VERSION_MINOR})
 endif()
+
 configure_file(${WRAP_DIR}/Python/__init__.py.in
     ${CMAKE_BINARY_DIR}/lib/bornagain/__init__.py @ONLY)
+foreach(mod plot_utils.py fit_monitor.py)
+    configure_file(${WRAP_DIR}/Python/${mod} ${CMAKE_BINARY_DIR}/lib/bornagain/${mod} COPYONLY)
+endforeach()
 
 if(CONFIGURE_BINDINGS)
     add_custom_command(
diff --git a/devtools/code-tools/batch-plot.py b/devtools/code-tools/batch-plot.py
new file mode 100755
index 0000000000000000000000000000000000000000..1059d28d19a3f5f14d6d39b4b17f2eccf38f4e7b
--- /dev/null
+++ b/devtools/code-tools/batch-plot.py
@@ -0,0 +1,100 @@
+#!/usr/bin/env python3
+"""
+Run plotting code in batch mode: generate figure, but don't plot to terminal.
+
+This script is used in BornAgain:
+- with option -s in PyExamples tests, to check functionality of Python examples;
+- with option -l to remake full-size images.
+"""
+
+import matplotlib, os, re, sys
+matplotlib.use('Agg')
+from matplotlib import pyplot as plt
+
+
+def exec_full(script, filename):
+    """
+    Executes embedded python script.
+    http://stackoverflow.com/questions/436198/what-is-an-alternative-to-execfile-in-python-3
+    """
+    import os
+    global_namespace = {
+        "__file__": filename,
+        "__name__": "__main__",
+        "__no_terminal__": True
+    }
+    sys.argv = []
+    exec(compile(script, filename, 'exec'), global_namespace)
+
+def reduce_nbin(t):
+    """
+    Overwrites script lines that set nbin, nx, ny
+    """
+    pat = re.compile(r'(^\s+(nbin|nx|ny) = )(\d+)$')
+    ret = []
+    for l in t.split('\n'):
+        m = re.match(pat, l)
+        if m:
+            oldsize = int(m.group(3))
+            newsize = min(10, oldsize)
+            lout = re.sub(pat, m.group(1)+f'{newsize}', l)
+        else:
+            lout = l
+        ret.append(lout)
+    return '\n'.join(ret)
+
+
+def run_example(mode_short, filename, output_dir):
+    """
+    Tries to run python example and produce a *.png image
+    """
+    # Read script from file.
+    if not os.path.exists(filename):
+        raise Exception("Example script '" + filename + "' not found")
+    print("Input script: " + filename, flush=True)
+    with open(filename, 'r') as file:
+        script = file.read()
+
+    # Detect or impose figure size.
+    m = re.search(r'plt\.figure\(.+?figsize=\((.+?),(.+?)\)', script)
+    if m: # set figure size as in script
+        figsize = (float(m.group(1)), float(m.group(2)))
+    else: # script does not specify figure size
+        figsize = (640/72, 480/72)
+    fig = plt.figure(figsize=(figsize[0], figsize[1]))
+
+    if mode_short:
+        script = reduce_nbin(script)
+
+    # Run modified script.
+    exec_full(script, filename)
+    print("Input script completed.", flush=True)
+
+    # Generate output figure.
+    plot_file_name = os.path.join(
+        output_dir,
+        os.path.splitext(os.path.basename(filename))[0] + ".png")
+    print("Output image: " + plot_file_name)
+    plt.savefig(plot_file_name, bbox_inches='tight')
+    plt.close(fig)
+
+    # Check obtained figure.
+    imgSize = os.path.getsize(plot_file_name)
+    if imgSize == 0:
+        raise Exception("Image file is empty")
+    if imgSize < 20000:
+        raise Exception("Image file is unplausibly small")
+
+
+if __name__ == '__main__':
+
+    if len(sys.argv) != 4:
+        exit("Usage: check_functionality -s|-f <script_to_test> <output_dir>")
+    if sys.argv[1]=='-s':
+        mode_short = True
+    elif sys.argv[1]=='-f':
+        mode_short = False
+    else:
+        exit("Use flag -s|-f for short or long runs")
+
+    run_example(mode_short, sys.argv[2], sys.argv[3])
diff --git a/devtools/code-tools/normalize-usercode.py b/devtools/code-tools/normalize-usercode.py
index c427a89801fe1296ef2e20c71a4c86a676f3a95b..10a0435382c65c525acc2026ac6b5b85c71d903d 100755
--- a/devtools/code-tools/normalize-usercode.py
+++ b/devtools/code-tools/normalize-usercode.py
@@ -69,6 +69,14 @@ def normalize_text(ti, fname):
             f'.. cycled through BA, {len(ti.split())} -> {len(tc.split())} lines'
         )
     tf = restitute_sample(ti, tc)
+
+        # restitute main
+    tf = re.sub(r"if __name__ == '__main__':.+",
+                """if __name__ == '__main__':
+    result = run_simulation()
+    ba.plot_simulation_result(result, cmap='jet', aspect='auto')""",
+                tf, flags=re.S)
+
     if verbose:
         print(f'.. normalized, {len(ti.split())} -> {len(tf.split())} lines')
     # YAPF formatting
@@ -88,6 +96,12 @@ def normalize_file(fname, inplace):
             if verbose:
                 print(f'.. read {len(ti.split())} lines')
 
+        m = re.search(r"""if __name__ == '__main__':
+    result = run_simulation\(\)
+    ba.plot_simulation_result\(result, cmap='jet', aspect='auto'\)""", ti)
+        if not m:
+            return 3
+
         # normalize
         tf = normalize_text(ti, fname)
         if verbose:
@@ -156,7 +170,7 @@ if __name__ == '__main__':
     verbose = args.verbose
     files = args.input_files
 
-    count = [0, 0, 0]
+    count = [0, 0, 0, 0]
     for f in files:
         ret = normalize_file(f, args.in_place)
         count[ret] += 1
@@ -171,4 +185,6 @@ if __name__ == '__main__':
         out.append(f'{count[1]} normalized')
     if count[2] > 0:
         out.append(f'{count[2]} failed')
+    if count[2] > 0:
+        out.append(f'{count[3]} skipped')
     print(f'TOTAL of {len(args.input_files)} files: {", ".join(out)}')