From ecca5d3ec2910957ca8b40d7167740ed1f03578b Mon Sep 17 00:00:00 2001
From: "Joachim Wuttke (o)" <j.wuttke@fz-juelich.de>
Date: Fri, 5 Aug 2016 16:24:14 +0200
Subject: [PATCH] Converted ex01-ex04; the others temporarily broken

---
 Core/CMakeLists.txt                           |   1 +
 .../ex07_FitAlongSlices/FitAlongSlices.py     |   2 +-
 .../ex10_FitGALAXIData/FitGALAXIData.py       |   2 +-
 .../AllFormFactorsAvailable.py                |   6 +-
 .../ex01_BasicParticles/CylindersAndPrisms.py |  30 +----
 .../ex01_BasicParticles/CylindersInBA.py      |  31 +----
 .../ex01_BasicParticles/CylindersInDWBA.py    |  22 +---
 .../CylindersWithSizeDistribution.py          |  21 +---
 .../ex01_BasicParticles/RotatedPyramids.py    |  21 +---
 ...TwoTypesOfCylindersWithSizeDistribution.py |  21 +---
 .../ex02_LayeredStructures/BuriedParticles.py |  22 +---
 .../CorrelatedRoughness.py                    |  22 +---
 .../ApproximationDA.py                        |  21 +---
 .../ApproximationLMA.py                       |  22 +---
 .../ApproximationSSCA.py                      |  21 +---
 .../CosineRipplesAtRectLattice.py             |  21 +---
 .../Interference1DLattice.py                  |  21 +---
 .../Interference1DRadialParaCrystal.py        |  21 +---
 .../Interference2DCenteredSquareLattice.py    |  21 +---
 .../Interference2DLatticeSumOfRotated.py      |  21 +---
 .../Interference2DParaCrystal.py              |  21 +---
 .../Interference2DRotatedSquareLattice.py     |  21 +---
 .../Interference2DSquareLattice.py            |  21 +---
 .../RectangularGrating.py                     |  21 +---
 .../SpheresAtHexLattice.py                    |  21 +---
 .../TriangularRipple.py                       |  21 +---
 .../CoreShellNanoparticles.py                 |  21 +---
 .../ex04_ComplexShapes/CustomFormFactor.py    |  21 +---
 .../HexagonalLatticesWithBasis.py             |  21 +---
 .../LargeParticlesFormFactor.py               | 110 ++++++++++--------
 .../ex05_BeamAndDetector/BeamDivergence.py    |   6 +-
 .../DetectorResolutionFunction.py             |   6 +-
 .../OffSpecularSimulation.py                  |   6 +-
 .../RectangularDetector.py                    |   4 +-
 .../SpecularSimulation.py                     |   4 +-
 .../AccessingSimulationResults.py             |   6 +-
 .../AxesInDifferentUnits.py                   |  12 +-
 Examples/python/utils/plot_intensity_data.py  |   2 +-
 .../python/utils/plot_intensity_data_diff.py  |   2 +-
 Examples/python/utils/show2d_root.py          |   2 +-
 Wrap/swig/CorePython.py                       |  22 +++-
 auto/Wrap/libBornAgainCore.py                 |  36 +++++-
 42 files changed, 222 insertions(+), 555 deletions(-)

diff --git a/Core/CMakeLists.txt b/Core/CMakeLists.txt
index a1d48e8828a..848472e5e2a 100644
--- a/Core/CMakeLists.txt
+++ b/Core/CMakeLists.txt
@@ -71,6 +71,7 @@ if(BORNAGAIN_PYTHON)
             ${WRAP_DIR}/swig/ignores.i
             ${WRAP_DIR}/swig/shared_pointers.i
             ${WRAP_DIR}/swig/warnings.i
+            ${WRAP_DIR}/swig/CorePython.py
             )
         foreach(FNAM ${swig_dependencies})
             if(NOT EXISTS ${FNAM})
diff --git a/Examples/python/fitting/ex07_FitAlongSlices/FitAlongSlices.py b/Examples/python/fitting/ex07_FitAlongSlices/FitAlongSlices.py
index 6b485d549a0..475038acfaf 100644
--- a/Examples/python/fitting/ex07_FitAlongSlices/FitAlongSlices.py
+++ b/Examples/python/fitting/ex07_FitAlongSlices/FitAlongSlices.py
@@ -9,7 +9,7 @@ import math
 import random
 import bornagain as ba
 from bornagain import deg, angstrom, nm
-import numpy, sys
+import numpy
 
 phi_slice_value = 0.0*deg  # position of vertical slice
 alpha_slice_value = 0.2*deg  # position of horizontal slice
diff --git a/Examples/python/fitting/ex10_FitGALAXIData/FitGALAXIData.py b/Examples/python/fitting/ex10_FitGALAXIData/FitGALAXIData.py
index 1bd6e7943fd..cde1ee93b72 100644
--- a/Examples/python/fitting/ex10_FitGALAXIData/FitGALAXIData.py
+++ b/Examples/python/fitting/ex10_FitGALAXIData/FitGALAXIData.py
@@ -3,7 +3,7 @@ Real life example: experiment at GALAXY
 """
 import matplotlib
 from matplotlib import pyplot as plt
-import numpy, sys
+import numpy
 import bornagain as ba
 from SampleBuilder import MySampleBuilder
 
diff --git a/Examples/python/simulation/ex01_BasicParticles/AllFormFactorsAvailable.py b/Examples/python/simulation/ex01_BasicParticles/AllFormFactorsAvailable.py
index 2686f86b886..6a7c182442e 100644
--- a/Examples/python/simulation/ex01_BasicParticles/AllFormFactorsAvailable.py
+++ b/Examples/python/simulation/ex01_BasicParticles/AllFormFactorsAvailable.py
@@ -1,7 +1,7 @@
 """
 All formfactors available in BornAgain in the Born Approximation
 """
-import numpy, sys
+import numpy
 import bornagain as ba
 from bornagain import deg, angstrom
 
@@ -114,7 +114,7 @@ if __name__ == '__main__':
         for nplot in range(len(formfactors)):
             ff = formfactors[nplot]
             name = ff.__class__.__name__.replace("FormFactor","")
-            print("Generating intensity map.or " + name)
+            print("Generating intensity map for " + name)
             intensities = simulate(ff)
             plot(intensities, nplot+1, name)
         plt.show()
@@ -124,4 +124,4 @@ if __name__ == '__main__':
             intensities = simulate(ff)
             fname = "%s.%s.int" % (sys.argv[1], name)
             ba.IntensityDataIOFactory.writeIntensityData(intensities, fname)
-            print("Stored intensity map.n " + fname)
+            print("Stored intensity map in " + fname)
diff --git a/Examples/python/simulation/ex01_BasicParticles/CylindersAndPrisms.py b/Examples/python/simulation/ex01_BasicParticles/CylindersAndPrisms.py
index cefa41da822..95c47f396c5 100644
--- a/Examples/python/simulation/ex01_BasicParticles/CylindersAndPrisms.py
+++ b/Examples/python/simulation/ex01_BasicParticles/CylindersAndPrisms.py
@@ -1,7 +1,7 @@
 """
 Mixture of cylinders and prisms without interference
 """
-import numpy, sys
+import numpy
 import bornagain as ba
 from bornagain import deg, angstrom, nm
 
@@ -61,32 +61,6 @@ def simulate():
     simulation.runSimulation()
     return simulation.getIntensityData()
 
-def plot(result):
-    """
-    Plots intensity map.
-    """
-    import matplotlib
-    from matplotlib import pyplot as plt
-    im = plt.imshow(
-        result.getArray(),
-        norm=matplotlib.colors.LogNorm(1.0, result.getMaximum()),
-        extent=[result.getXmin()/deg, result.getXmax()/deg,
-                result.getYmin()/deg, result.getYmax()/deg],
-        aspect='auto')
-    cb = plt.colorbar(im)
-    cb.set_label(r'Intensity (arb. u.)', size=16)
-    plt.xlabel(r'$\phi_f (^{\circ})$', fontsize=16)
-    plt.ylabel(r'$\alpha_f (^{\circ})$', fontsize=16)
-    plt.show()
 
 if __name__ == '__main__':
-    if len(sys.argv)<=1:
-        print('Usage:')
-        print('    '+sys.argv[0]+' -p      # to plot result')
-        print('    '+sys.argv[0]+' <fname> # to save result in <fname>.int')
-        sys.exit(1)
-    result = simulate()
-    if sys.argv[1] != '-p':
-        ba.IntensityDataIOFactory.writeIntensityData(result, sys.argv[1]+".int")
-    else:
-        plot(result)
+    ba.simulateThenPlotOrSave(simulate, standardIntensityPlot)
diff --git a/Examples/python/simulation/ex01_BasicParticles/CylindersInBA.py b/Examples/python/simulation/ex01_BasicParticles/CylindersInBA.py
index 5c1335deab3..df65592d9e0 100644
--- a/Examples/python/simulation/ex01_BasicParticles/CylindersInBA.py
+++ b/Examples/python/simulation/ex01_BasicParticles/CylindersInBA.py
@@ -1,7 +1,7 @@
 """
 Cylinder formfactor in Born approximation
 """
-import numpy, sys
+import numpy
 import bornagain as ba
 from bornagain import deg, angstrom, nm
 
@@ -54,32 +54,5 @@ def simulate():
     return simulation.getIntensityData()
 
 
-def plot(result):
-    """
-    Plots intensity map.
-    """
-    import matplotlib
-    from matplotlib import pyplot as plt
-    im = plt.imshow(
-        result.getArray(),
-        norm=matplotlib.colors.LogNorm(1.0, result.getMaximum()),
-        extent=[result.getXmin()/deg, result.getXmax()/deg,
-                result.getYmin()/deg, result.getYmax()/deg],
-        aspect='auto')
-    cb = plt.colorbar(im)
-    cb.set_label(r'Intensity (arb. u.)', size=16)
-    plt.xlabel(r'$\phi_f (^{\circ})$', fontsize=16)
-    plt.ylabel(r'$\alpha_f (^{\circ})$', fontsize=16)
-    plt.show()
-
 if __name__ == '__main__':
-    if len(sys.argv)<=1:
-        print('Usage:')
-        print('    '+sys.argv[0]+' -p      # to plot result')
-        print('    '+sys.argv[0]+' <fname> # to save result in <fname>.int')
-        sys.exit(1)
-    result = simulate()
-    if sys.argv[1] != '-p':
-        ba.IntensityDataIOFactory.writeIntensityData(result, sys.argv[1]+".int")
-    else:
-        plot(result)
+    ba.simulateThenPlotOrSave(simulate, ba.standardIntensityPlot)
diff --git a/Examples/python/simulation/ex01_BasicParticles/CylindersInDWBA.py b/Examples/python/simulation/ex01_BasicParticles/CylindersInDWBA.py
index a069b24920a..06387541b0d 100644
--- a/Examples/python/simulation/ex01_BasicParticles/CylindersInDWBA.py
+++ b/Examples/python/simulation/ex01_BasicParticles/CylindersInDWBA.py
@@ -1,7 +1,7 @@
 """
 Cylinder formfactor in DWBA
 """
-import numpy, sys
+import numpy
 import bornagain as ba
 from bornagain import deg, angstrom, nm
 
@@ -55,24 +55,6 @@ def simulate():
     simulation.runSimulation()
     return simulation.getIntensityData()
 
-def plot(result):
-    """
-    Plots intensity map.
-    """
-    import matplotlib
-    from matplotlib import pyplot as plt
-    im = plt.imshow(
-        result.getArray(),
-        norm=matplotlib.colors.LogNorm(1.0, result.getMaximum()),
-        extent=[result.getXmin()/deg, result.getXmax()/deg,
-                result.getYmin()/deg, result.getYmax()/deg],
-        aspect='auto')
-    cb = plt.colorbar(im)
-    cb.set_label(r'Intensity (arb. u.)', size=16)
-    plt.xlabel(r'$\phi_f (^{\circ})$', fontsize=16)
-    plt.ylabel(r'$\alpha_f (^{\circ})$', fontsize=16)
-    plt.show()
-
 
 if __name__ == '__main__':
-    ba.simulate_then_plot_or_save(simulate, plot)
+    ba.simulateThenPlotOrSave(simulate, standardIntensityPlot)
diff --git a/Examples/python/simulation/ex01_BasicParticles/CylindersWithSizeDistribution.py b/Examples/python/simulation/ex01_BasicParticles/CylindersWithSizeDistribution.py
index d6325506770..ac2849c08d3 100644
--- a/Examples/python/simulation/ex01_BasicParticles/CylindersWithSizeDistribution.py
+++ b/Examples/python/simulation/ex01_BasicParticles/CylindersWithSizeDistribution.py
@@ -1,9 +1,7 @@
 """
 Cylinders with size distribution
 """
-import numpy, sys
-import matplotlib
-from matplotlib import pyplot as plt
+import numpy
 import bornagain as ba
 from bornagain import deg, angstrom, nm
 
@@ -69,21 +67,8 @@ def simulate():
     simulation = get_simulation()
     simulation.setSample(sample)
     simulation.runSimulation()
-    result = simulation.getIntensityData()
-
-    # showing the result
-    im = plt.imshow(
-        result.getArray(),
-        norm=matplotlib.colors.LogNorm(1.0, result.getMaximum()),
-        extent=[result.getXmin()/deg, result.getXmax()/deg,
-                result.getYmin()/deg, result.getYmax()/deg],
-        aspect='auto')
-    cb = plt.colorbar(im)
-    cb.set_label(r'Intensity (arb. u.)', size=16)
-    plt.xlabel(r'$\phi_f (^{\circ})$', fontsize=16)
-    plt.ylabel(r'$\alpha_f (^{\circ})$', fontsize=16)
-    plt.show()
+    return simulation.getIntensityData()
 
 
 if __name__ == '__main__':
-    simulate()
+    ba.simulateThenPlotOrSave(simulate, ba.standardIntensityPlot)
diff --git a/Examples/python/simulation/ex01_BasicParticles/RotatedPyramids.py b/Examples/python/simulation/ex01_BasicParticles/RotatedPyramids.py
index 6c2ca1e52cc..f3254acd119 100644
--- a/Examples/python/simulation/ex01_BasicParticles/RotatedPyramids.py
+++ b/Examples/python/simulation/ex01_BasicParticles/RotatedPyramids.py
@@ -1,9 +1,7 @@
 """
 Rotated pyramids on top of substrate
 """
-import numpy, sys
-import matplotlib
-from matplotlib import pyplot as plt
+import numpy
 import bornagain as ba
 from bornagain import deg, angstrom, nm
 
@@ -57,21 +55,8 @@ def simulate():
     simulation = get_simulation()
     simulation.setSample(sample)
     simulation.runSimulation()
-    result = simulation.getIntensityData()
-
-    # showing the result
-    im = plt.imshow(
-        result.getArray(),
-        norm=matplotlib.colors.LogNorm(1.0, result.getMaximum()),
-        extent=[result.getXmin()/deg, result.getXmax()/deg,
-                result.getYmin()/deg, result.getYmax()/deg],
-        aspect='auto')
-    cb = plt.colorbar(im)
-    cb.set_label(r'Intensity (arb. u.)', size=16)
-    plt.xlabel(r'$\phi_f (^{\circ})$', fontsize=16)
-    plt.ylabel(r'$\alpha_f (^{\circ})$', fontsize=16)
-    plt.show()
+    return simulation.getIntensityData()
 
 
 if __name__ == '__main__':
-    simulate()
+    ba.simulateThenPlotOrSave(simulate, ba.standardIntensityPlot)
diff --git a/Examples/python/simulation/ex01_BasicParticles/TwoTypesOfCylindersWithSizeDistribution.py b/Examples/python/simulation/ex01_BasicParticles/TwoTypesOfCylindersWithSizeDistribution.py
index 00cc0cbeba4..a0af5935bfd 100644
--- a/Examples/python/simulation/ex01_BasicParticles/TwoTypesOfCylindersWithSizeDistribution.py
+++ b/Examples/python/simulation/ex01_BasicParticles/TwoTypesOfCylindersWithSizeDistribution.py
@@ -1,9 +1,7 @@
 """
 Mixture cylinder particles with different size distribution
 """
-import numpy, sys
-import matplotlib
-from matplotlib import pyplot as plt
+import numpy
 import bornagain as ba
 from bornagain import deg, angstrom, nm
 
@@ -83,21 +81,8 @@ def simulate():
     simulation = get_simulation()
     simulation.setSample(sample)
     simulation.runSimulation()
-    result = simulation.getIntensityData()
-
-    # showing the result
-    im = plt.imshow(
-        result.getArray(),
-        norm=matplotlib.colors.LogNorm(1.0, result.getMaximum()),
-        extent=[result.getXmin()/deg, result.getXmax()/deg,
-                result.getYmin()/deg, result.getYmax()/deg],
-        aspect='auto')
-    cb = plt.colorbar(im)
-    cb.set_label(r'Intensity (arb. u.)', size=16)
-    plt.xlabel(r'$\phi_f (^{\circ})$', fontsize=16)
-    plt.ylabel(r'$\alpha_f (^{\circ})$', fontsize=16)
-    plt.show()
+    return simulation.getIntensityData()
 
 
 if __name__ == '__main__':
-    simulate()
+    ba.simulateThenPlotOrSave(simulate, ba.standardIntensityPlot)
diff --git a/Examples/python/simulation/ex02_LayeredStructures/BuriedParticles.py b/Examples/python/simulation/ex02_LayeredStructures/BuriedParticles.py
index 10cd4b2c87b..a2a8fd5d892 100644
--- a/Examples/python/simulation/ex02_LayeredStructures/BuriedParticles.py
+++ b/Examples/python/simulation/ex02_LayeredStructures/BuriedParticles.py
@@ -1,9 +1,7 @@
 """
 Spherical particles embedded in the middle of the layer on top of substrate.
 """
-import numpy, sys
-import matplotlib
-from matplotlib import pyplot as plt
+import numpy
 import bornagain as ba
 from bornagain import deg, angstrom, nm
 
@@ -60,21 +58,7 @@ def simulate():
     simulation = get_simulation()
     simulation.setSample(sample)
     simulation.runSimulation()
-    result = simulation.getIntensityData()
-
-    # showing the result
-    im = plt.imshow(
-        result.getArray(),
-        norm=matplotlib.colors.LogNorm(1.0, result.getMaximum()),
-        extent=[result.getXmin()/deg, result.getXmax()/deg,
-                result.getYmin()/deg, result.getYmax()/deg],
-        aspect='auto')
-    cb = plt.colorbar(im)
-    cb.set_label(r'Intensity (arb. u.)', size=16)
-    plt.xlabel(r'$\phi_f (^{\circ})$', fontsize=16)
-    plt.ylabel(r'$\alpha_f (^{\circ})$', fontsize=16)
-    plt.show()
-
+    return simulation.getIntensityData()
 
 if __name__ == '__main__':
-    simulate()
+    ba.simulateThenPlotOrSave(simulate, ba.standardIntensityPlot)
diff --git a/Examples/python/simulation/ex02_LayeredStructures/CorrelatedRoughness.py b/Examples/python/simulation/ex02_LayeredStructures/CorrelatedRoughness.py
index fd549e6be02..9b5d61cc93e 100644
--- a/Examples/python/simulation/ex02_LayeredStructures/CorrelatedRoughness.py
+++ b/Examples/python/simulation/ex02_LayeredStructures/CorrelatedRoughness.py
@@ -1,9 +1,7 @@
 """
 MultiLayer with correlated roughness
 """
-import numpy, sys
-import matplotlib
-from matplotlib import pyplot as plt
+import numpy
 import bornagain as ba
 from bornagain import deg, angstrom, nm
 
@@ -67,22 +65,8 @@ def simulate():
     simulation = get_simulation()
     simulation.setSample(sample)
     simulation.runSimulation()
-    result = simulation.getIntensityData()
-
-    # showing the result
-    im = plt.imshow(
-        result.getArray(),
-        norm=matplotlib.colors.LogNorm(result.getMaximum()/1000.,
-                                       result.getMaximum()),
-        extent=[result.getXmin()/deg, result.getXmax()/deg,
-                result.getYmin()/deg, result.getYmax()/deg],
-        aspect='auto')
-    cb = plt.colorbar(im)
-    cb.set_label(r'Intensity (arb. u.)', size=16)
-    plt.xlabel(r'$\phi_f (^{\circ})$', fontsize=16)
-    plt.ylabel(r'$\alpha_f (^{\circ})$', fontsize=16)
-    plt.show()
+    return simulation.getIntensityData()
 
 
 if __name__ == '__main__':
-    simulate()
+    ba.simulateThenPlotOrSave(simulate, ba.standardIntensityPlot)
diff --git a/Examples/python/simulation/ex03_InterferenceFunctions/ApproximationDA.py b/Examples/python/simulation/ex03_InterferenceFunctions/ApproximationDA.py
index 8c2779fec36..3980879060f 100644
--- a/Examples/python/simulation/ex03_InterferenceFunctions/ApproximationDA.py
+++ b/Examples/python/simulation/ex03_InterferenceFunctions/ApproximationDA.py
@@ -1,9 +1,7 @@
 """
 Cylinders of two different sizes in Decoupling Approximation
 """
-import numpy, sys
-import matplotlib
-from matplotlib import pyplot as plt
+import numpy
 import bornagain as ba
 from bornagain import deg, angstrom, nm
 
@@ -72,21 +70,8 @@ def simulate():
     simulation = get_simulation()
     simulation.setSample(sample)
     simulation.runSimulation()
-    result = simulation.getIntensityData()
-
-    # showing the result
-    im = plt.imshow(
-        result.getArray(),
-        norm=matplotlib.colors.LogNorm(1.0, result.getMaximum()),
-        extent=[result.getXmin()/deg, result.getXmax()/deg,
-                result.getYmin()/deg, result.getYmax()/deg],
-        aspect='auto')
-    cb = plt.colorbar(im)
-    cb.set_label(r'Intensity (arb. u.)', size=16)
-    plt.xlabel(r'$\phi_f (^{\circ})$', fontsize=16)
-    plt.ylabel(r'$\alpha_f (^{\circ})$', fontsize=16)
-    plt.show()
+    return simulation.getIntensityData()
 
 
 if __name__ == '__main__':
-    simulate()
+    ba.simulateThenPlotOrSave(simulate, ba.standardIntensityPlot)
diff --git a/Examples/python/simulation/ex03_InterferenceFunctions/ApproximationLMA.py b/Examples/python/simulation/ex03_InterferenceFunctions/ApproximationLMA.py
index 78482676904..5b248db6887 100644
--- a/Examples/python/simulation/ex03_InterferenceFunctions/ApproximationLMA.py
+++ b/Examples/python/simulation/ex03_InterferenceFunctions/ApproximationLMA.py
@@ -1,9 +1,7 @@
 """
 Cylinders of two different sizes in Local Monodisperse Approximation
 """
-import numpy, sys
-import matplotlib
-from matplotlib import pyplot as plt
+import numpy
 import bornagain as ba
 from bornagain import deg, angstrom, nm
 
@@ -81,21 +79,7 @@ def simulate():
     simulation = get_simulation()
     simulation.setSample(sample)
     simulation.runSimulation()
-    result = simulation.getIntensityData()
-
-    # showing the result
-    im = plt.imshow(
-        result.getArray(),
-        norm=matplotlib.colors.LogNorm(1.0, result.getMaximum()),
-        extent=[result.getXmin()/deg, result.getXmax()/deg,
-                result.getYmin()/deg, result.getYmax()/deg],
-        aspect='auto')
-    cb = plt.colorbar(im)
-    cb.set_label(r'Intensity (arb. u.)', size=16)
-    plt.xlabel(r'$\phi_f (^{\circ})$', fontsize=16)
-    plt.ylabel(r'$\alpha_f (^{\circ})$', fontsize=16)
-    plt.show()
-
+    return simulation.getIntensityData()
 
 if __name__ == '__main__':
-    simulate()
+    ba.simulateThenPlotOrSave(simulate, ba.standardIntensityPlot)
diff --git a/Examples/python/simulation/ex03_InterferenceFunctions/ApproximationSSCA.py b/Examples/python/simulation/ex03_InterferenceFunctions/ApproximationSSCA.py
index aa688431c4f..4c78b21b133 100644
--- a/Examples/python/simulation/ex03_InterferenceFunctions/ApproximationSSCA.py
+++ b/Examples/python/simulation/ex03_InterferenceFunctions/ApproximationSSCA.py
@@ -1,9 +1,7 @@
 """
 Cylinders of two different sizes in Size-Spacing Coupling Approximation
 """
-import numpy, sys
-import matplotlib
-from matplotlib import pyplot as plt
+import numpy
 import bornagain as ba
 from bornagain import deg, angstrom, nm
 
@@ -74,21 +72,8 @@ def simulate():
     simulation = get_simulation()
     simulation.setSample(sample)
     simulation.runSimulation()
-    result = simulation.getIntensityData()
-
-    # showing the result
-    im = plt.imshow(
-        result.getArray(),
-        norm=matplotlib.colors.LogNorm(1.0, result.getMaximum()),
-        extent=[result.getXmin()/deg, result.getXmax()/deg,
-                result.getYmin()/deg, result.getYmax()/deg],
-        aspect='auto')
-    cb = plt.colorbar(im)
-    cb.set_label(r'Intensity (arb. u.)', size=16)
-    plt.xlabel(r'$\phi_f (^{\circ})$', fontsize=16)
-    plt.ylabel(r'$\alpha_f (^{\circ})$', fontsize=16)
-    plt.show()
+    return simulation.getIntensityData()
 
 
 if __name__ == '__main__':
-    simulate()
+    ba.simulateThenPlotOrSave(simulate, ba.standardIntensityPlot)
diff --git a/Examples/python/simulation/ex03_InterferenceFunctions/CosineRipplesAtRectLattice.py b/Examples/python/simulation/ex03_InterferenceFunctions/CosineRipplesAtRectLattice.py
index cd2ed9e273a..78da8ef6e9f 100644
--- a/Examples/python/simulation/ex03_InterferenceFunctions/CosineRipplesAtRectLattice.py
+++ b/Examples/python/simulation/ex03_InterferenceFunctions/CosineRipplesAtRectLattice.py
@@ -1,9 +1,7 @@
 """
 Cosine ripple on a 2D lattice
 """
-import numpy, sys
-import matplotlib
-from matplotlib import pyplot as plt
+import numpy
 import bornagain as ba
 from bornagain import deg, angstrom, nm
 
@@ -65,21 +63,8 @@ def simulate():
     simulation = get_simulation()
     simulation.setSample(sample)
     simulation.runSimulation()
-    result = simulation.getIntensityData()
-
-    # showing the result
-    im = plt.imshow(
-        result.getArray(),
-        norm=matplotlib.colors.LogNorm(1.0, result.getMaximum()),
-        extent=[result.getXmin()/deg, result.getXmax()/deg,
-                result.getYmin()/deg, result.getYmax()/deg],
-        aspect='auto')
-    cb = plt.colorbar(im)
-    cb.set_label(r'Intensity (arb. u.)', size=16)
-    plt.xlabel(r'$\phi_f (^{\circ})$', fontsize=16)
-    plt.ylabel(r'$\alpha_f (^{\circ})$', fontsize=16)
-    plt.show()
+    return simulation.getIntensityData()
 
 
 if __name__ == '__main__':
-    simulate()
+    ba.simulateThenPlotOrSave(simulate, ba.standardIntensityPlot)
diff --git a/Examples/python/simulation/ex03_InterferenceFunctions/Interference1DLattice.py b/Examples/python/simulation/ex03_InterferenceFunctions/Interference1DLattice.py
index 55f72df30ac..b3359474c1a 100644
--- a/Examples/python/simulation/ex03_InterferenceFunctions/Interference1DLattice.py
+++ b/Examples/python/simulation/ex03_InterferenceFunctions/Interference1DLattice.py
@@ -1,9 +1,7 @@
 """
 Long boxes on a 1D lattice
 """
-import numpy, sys
-import matplotlib
-from matplotlib import pyplot as plt
+import numpy
 import bornagain as ba
 from bornagain import deg, angstrom, nm
 
@@ -66,21 +64,8 @@ def simulate():
     simulation = get_simulation()
     simulation.setSample(sample)
     simulation.runSimulation()
-    result = simulation.getIntensityData()
-
-    # showing the result
-    im = plt.imshow(
-        result.getArray(),
-        norm=matplotlib.colors.LogNorm(1.0, result.getMaximum()),
-        extent=[result.getXmin()/deg, result.getXmax()/deg,
-                result.getYmin()/deg, result.getYmax()/deg],
-        aspect='auto')
-    cb = plt.colorbar(im)
-    cb.set_label(r'Intensity (arb. u.)', size=16)
-    plt.xlabel(r'$\phi_f (^{\circ})$', fontsize=16)
-    plt.ylabel(r'$\alpha_f (^{\circ})$', fontsize=16)
-    plt.show()
+    return simulation.getIntensityData()
 
 
 if __name__ == '__main__':
-    simulate()
+    ba.simulateThenPlotOrSave(simulate, ba.standardIntensityPlot)
diff --git a/Examples/python/simulation/ex03_InterferenceFunctions/Interference1DRadialParaCrystal.py b/Examples/python/simulation/ex03_InterferenceFunctions/Interference1DRadialParaCrystal.py
index 332e18c6782..a2f8db7a5c2 100644
--- a/Examples/python/simulation/ex03_InterferenceFunctions/Interference1DRadialParaCrystal.py
+++ b/Examples/python/simulation/ex03_InterferenceFunctions/Interference1DRadialParaCrystal.py
@@ -1,9 +1,7 @@
 """
 radial paracrystal
 """
-import numpy, sys
-import matplotlib
-from matplotlib import pyplot as plt
+import numpy
 import bornagain as ba
 from bornagain import deg, angstrom, nm
 
@@ -63,21 +61,8 @@ def simulate():
     simulation = get_simulation()
     simulation.setSample(sample)
     simulation.runSimulation()
-    result = simulation.getIntensityData()
-
-    # showing the result
-    im = plt.imshow(
-        result.getArray(),
-        norm=matplotlib.colors.LogNorm(1.0, result.getMaximum()),
-        extent=[result.getXmin()/deg, result.getXmax()/deg,
-                result.getYmin()/deg, result.getYmax()/deg],
-        aspect='auto')
-    cb = plt.colorbar(im)
-    cb.set_label(r'Intensity (arb. u.)', size=16)
-    plt.xlabel(r'$\phi_f (^{\circ})$', fontsize=16)
-    plt.ylabel(r'$\alpha_f (^{\circ})$', fontsize=16)
-    plt.show()
+    return simulation.getIntensityData()
 
 
 if __name__ == '__main__':
-    simulate()
+    ba.simulateThenPlotOrSave(simulate, ba.standardIntensityPlot)
diff --git a/Examples/python/simulation/ex03_InterferenceFunctions/Interference2DCenteredSquareLattice.py b/Examples/python/simulation/ex03_InterferenceFunctions/Interference2DCenteredSquareLattice.py
index 6dfcd3050f1..bbfc555456d 100644
--- a/Examples/python/simulation/ex03_InterferenceFunctions/Interference2DCenteredSquareLattice.py
+++ b/Examples/python/simulation/ex03_InterferenceFunctions/Interference2DCenteredSquareLattice.py
@@ -1,9 +1,7 @@
 """
 2D lattice with disorder, centered square lattice
 """
-import numpy, sys
-import matplotlib
-from matplotlib import pyplot as plt
+import numpy
 import bornagain as ba
 from bornagain import deg, angstrom, nm
 
@@ -68,21 +66,8 @@ def simulate():
     simulation = get_simulation()
     simulation.setSample(sample)
     simulation.runSimulation()
-    result = simulation.getIntensityData()
-
-    # showing the result
-    im = plt.imshow(
-        result.getArray(),
-        norm=matplotlib.colors.LogNorm(1.0, result.getMaximum()),
-        extent=[result.getXmin()/deg, result.getXmax()/deg,
-                result.getYmin()/deg, result.getYmax()/deg],
-        aspect='auto')
-    cb = plt.colorbar(im)
-    cb.set_label(r'Intensity (arb. u.)', size=16)
-    plt.xlabel(r'$\phi_f (^{\circ})$', fontsize=16)
-    plt.ylabel(r'$\alpha_f (^{\circ})$', fontsize=16)
-    plt.show()
+    return simulation.getIntensityData()
 
 
 if __name__ == '__main__':
-    simulate()
+    ba.simulateThenPlotOrSave(simulate, ba.standardIntensityPlot)
diff --git a/Examples/python/simulation/ex03_InterferenceFunctions/Interference2DLatticeSumOfRotated.py b/Examples/python/simulation/ex03_InterferenceFunctions/Interference2DLatticeSumOfRotated.py
index 4b0507c2cf5..074eb75f171 100644
--- a/Examples/python/simulation/ex03_InterferenceFunctions/Interference2DLatticeSumOfRotated.py
+++ b/Examples/python/simulation/ex03_InterferenceFunctions/Interference2DLatticeSumOfRotated.py
@@ -1,7 +1,5 @@
 # 2D lattice with different disorder (IsGISAXS example #6), sum of rotated lattices
-import numpy, sys
-import matplotlib
-from matplotlib import pyplot as plt
+import numpy
 import bornagain as ba
 from bornagain import deg, angstrom, nm
 
@@ -83,21 +81,8 @@ def simulate():
         OutputData_total += single_output
     OutputData_total.scale(1.0/total_weight)
 
-    result = OutputData_total
-
-    # showing the result
-    im = plt.imshow(
-        result.getArray(),
-        norm=matplotlib.colors.LogNorm(1.0, result.getMaximum()),
-        extent=[result.getXmin()/deg, result.getXmax()/deg,
-                result.getYmin()/deg, result.getYmax()/deg],
-        aspect='auto')
-    cb = plt.colorbar(im)
-    cb.set_label(r'Intensity (arb. u.)', size=16)
-    plt.xlabel(r'$\phi_f (^{\circ})$', fontsize=16)
-    plt.ylabel(r'$\alpha_f (^{\circ})$', fontsize=16)
-    plt.show()
+    return OutputData_total
 
 
 if __name__ == '__main__':
-    simulate()
+    ba.simulateThenPlotOrSave(simulate, ba.standardIntensityPlot)
diff --git a/Examples/python/simulation/ex03_InterferenceFunctions/Interference2DParaCrystal.py b/Examples/python/simulation/ex03_InterferenceFunctions/Interference2DParaCrystal.py
index ed17848a097..38cbd8e4fbc 100644
--- a/Examples/python/simulation/ex03_InterferenceFunctions/Interference2DParaCrystal.py
+++ b/Examples/python/simulation/ex03_InterferenceFunctions/Interference2DParaCrystal.py
@@ -1,9 +1,7 @@
 """
 2D paracrystal
 """
-import numpy, sys
-import matplotlib
-from matplotlib import pyplot as plt
+import numpy
 import bornagain as ba
 from bornagain import deg, angstrom, nm
 from bornagain import micrometer
@@ -64,21 +62,8 @@ def simulate():
     simulation = get_simulation()
     simulation.setSample(sample)
     simulation.runSimulation()
-    result = simulation.getIntensityData()
-
-    # showing the result
-    im = plt.imshow(
-        result.getArray(),
-        norm=matplotlib.colors.LogNorm(1.0, result.getMaximum()),
-        extent=[result.getXmin()/deg, result.getXmax()/deg,
-                result.getYmin()/deg, result.getYmax()/deg],
-        aspect='auto')
-    cb = plt.colorbar(im)
-    cb.set_label(r'Intensity (arb. u.)', size=16)
-    plt.xlabel(r'$\phi_f (^{\circ})$', fontsize=16)
-    plt.ylabel(r'$\alpha_f (^{\circ})$', fontsize=16)
-    plt.show()
+    return simulation.getIntensityData()
 
 
 if __name__ == '__main__':
-    simulate()
+    ba.simulateThenPlotOrSave(simulate, ba.standardIntensityPlot)
diff --git a/Examples/python/simulation/ex03_InterferenceFunctions/Interference2DRotatedSquareLattice.py b/Examples/python/simulation/ex03_InterferenceFunctions/Interference2DRotatedSquareLattice.py
index 6268332b237..7da27bb23a9 100644
--- a/Examples/python/simulation/ex03_InterferenceFunctions/Interference2DRotatedSquareLattice.py
+++ b/Examples/python/simulation/ex03_InterferenceFunctions/Interference2DRotatedSquareLattice.py
@@ -1,9 +1,7 @@
 """
 Cylinders on a rotated 2D lattice
 """
-import numpy, sys
-import matplotlib
-from matplotlib import pyplot as plt
+import numpy
 import bornagain as ba
 from bornagain import deg, angstrom, nm
 
@@ -65,21 +63,8 @@ def simulate():
     simulation = get_simulation()
     simulation.setSample(sample)
     simulation.runSimulation()
-    result = simulation.getIntensityData()
-
-    # showing the result
-    im = plt.imshow(
-        result.getArray(),
-        norm=matplotlib.colors.LogNorm(1.0, result.getMaximum()),
-        extent=[result.getXmin()/deg, result.getXmax()/deg,
-                result.getYmin()/deg, result.getYmax()/deg],
-        aspect='auto')
-    cb = plt.colorbar(im)
-    cb.set_label(r'Intensity (arb. u.)', size=16)
-    plt.xlabel(r'$\phi_f (^{\circ})$', fontsize=16)
-    plt.ylabel(r'$\alpha_f (^{\circ})$', fontsize=16)
-    plt.show()
+    return simulation.getIntensityData()
 
 
 if __name__ == '__main__':
-    simulate()
+    ba.simulateThenPlotOrSave(simulate, ba.standardIntensityPlot)
diff --git a/Examples/python/simulation/ex03_InterferenceFunctions/Interference2DSquareLattice.py b/Examples/python/simulation/ex03_InterferenceFunctions/Interference2DSquareLattice.py
index 5910e22f015..92c9b115ee1 100644
--- a/Examples/python/simulation/ex03_InterferenceFunctions/Interference2DSquareLattice.py
+++ b/Examples/python/simulation/ex03_InterferenceFunctions/Interference2DSquareLattice.py
@@ -1,9 +1,7 @@
 """
 Cylinders on a 2D square lattice
 """
-import numpy, sys
-import matplotlib
-from matplotlib import pyplot as plt
+import numpy
 import bornagain as ba
 from bornagain import deg, angstrom, nm
 
@@ -62,21 +60,8 @@ def simulate():
     simulation = get_simulation()
     simulation.setSample(sample)
     simulation.runSimulation()
-    result = simulation.getIntensityData()
-
-    # showing the result
-    im = plt.imshow(
-        result.getArray(),
-        norm=matplotlib.colors.LogNorm(1.0, result.getMaximum()),
-        extent=[result.getXmin()/deg, result.getXmax()/deg,
-                result.getYmin()/deg, result.getYmax()/deg],
-        aspect='auto')
-    cb = plt.colorbar(im)
-    cb.set_label(r'Intensity (arb. u.)', size=16)
-    plt.xlabel(r'$\phi_f (^{\circ})$', fontsize=16)
-    plt.ylabel(r'$\alpha_f (^{\circ})$', fontsize=16)
-    plt.show()
+    return simulation.getIntensityData()
 
 
 if __name__ == '__main__':
-    simulate()
+    ba.simulateThenPlotOrSave(simulate, ba.standardIntensityPlot)
diff --git a/Examples/python/simulation/ex03_InterferenceFunctions/RectangularGrating.py b/Examples/python/simulation/ex03_InterferenceFunctions/RectangularGrating.py
index a0486576c2e..d3b6aef158d 100644
--- a/Examples/python/simulation/ex03_InterferenceFunctions/RectangularGrating.py
+++ b/Examples/python/simulation/ex03_InterferenceFunctions/RectangularGrating.py
@@ -3,9 +3,7 @@ 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 matplotlib
-from matplotlib import pyplot as plt
+import numpy
 import bornagain as ba
 from bornagain import deg, angstrom, nm
 
@@ -72,21 +70,8 @@ def simulate():
     simulation = get_simulation(monte_carlo_integration=True)
     simulation.setSample(sample)
     simulation.runSimulation()
-    result = simulation.getIntensityData()
-
-    # showing the result
-    im = plt.imshow(
-        result.getArray(),
-        norm=matplotlib.colors.LogNorm(1.0, result.getMaximum()),
-        extent=[result.getXmin()/deg, result.getXmax()/deg,
-                result.getYmin()/deg, result.getYmax()/deg],
-        aspect='auto')
-    cb = plt.colorbar(im)
-    cb.set_label(r'Intensity (arb. u.)', size=16)
-    plt.xlabel(r'$\phi_f (^{\circ})$', fontsize=16)
-    plt.ylabel(r'$\alpha_f (^{\circ})$', fontsize=16)
-    plt.show()
+    return simulation.getIntensityData()
 
 
 if __name__ == '__main__':
-    simulate()
+    ba.simulateThenPlotOrSave(simulate, ba.standardIntensityPlot)
diff --git a/Examples/python/simulation/ex03_InterferenceFunctions/SpheresAtHexLattice.py b/Examples/python/simulation/ex03_InterferenceFunctions/SpheresAtHexLattice.py
index 93f99af6e34..1fede9b2f0b 100644
--- a/Examples/python/simulation/ex03_InterferenceFunctions/SpheresAtHexLattice.py
+++ b/Examples/python/simulation/ex03_InterferenceFunctions/SpheresAtHexLattice.py
@@ -1,9 +1,7 @@
 """
 Spheres on a hexagonal lattice
 """
-import numpy, sys
-import matplotlib
-from matplotlib import pyplot as plt
+import numpy
 import bornagain as ba
 from bornagain import deg, angstrom, nm
 
@@ -59,21 +57,8 @@ def simulate():
     simulation = get_simulation()
     simulation.setSample(sample)
     simulation.runSimulation()
-    result = simulation.getIntensityData()
-
-    # showing the result
-    im = plt.imshow(
-        result.getArray(),
-        norm=matplotlib.colors.LogNorm(1.0, result.getMaximum()),
-        extent=[result.getXmin()/deg, result.getXmax()/deg,
-                result.getYmin()/deg, result.getYmax()/deg],
-        aspect='auto')
-    cb = plt.colorbar(im)
-    cb.set_label(r'Intensity (arb. u.)', size=16)
-    plt.xlabel(r'$\phi_f (^{\circ})$', fontsize=16)
-    plt.ylabel(r'$\alpha_f (^{\circ})$', fontsize=16)
-    plt.show()
+    return simulation.getIntensityData()
 
 
 if __name__ == '__main__':
-    simulate()
+    ba.simulateThenPlotOrSave(simulate, ba.standardIntensityPlot)
diff --git a/Examples/python/simulation/ex03_InterferenceFunctions/TriangularRipple.py b/Examples/python/simulation/ex03_InterferenceFunctions/TriangularRipple.py
index 6dfdb249702..eb41a7cdb8b 100644
--- a/Examples/python/simulation/ex03_InterferenceFunctions/TriangularRipple.py
+++ b/Examples/python/simulation/ex03_InterferenceFunctions/TriangularRipple.py
@@ -1,9 +1,7 @@
 """
  Sample from the article D. Babonneau et. al., Phys. Rev. B 85, 235415, 2012 (Fig.3)
 """
-import numpy, sys
-import matplotlib
-from matplotlib import pyplot as plt
+import numpy
 import bornagain as ba
 from bornagain import deg, angstrom, nm
 
@@ -66,21 +64,8 @@ def simulate():
     simulation = get_simulation()
     simulation.setSample(sample)
     simulation.runSimulation()
-    result = simulation.getIntensityData()
-
-    # showing the result
-    im = plt.imshow(
-        result.getArray(),
-        norm=matplotlib.colors.LogNorm(1.0, result.getMaximum()),
-        extent=[result.getXmin()/deg, result.getXmax()/deg,
-                result.getYmin()/deg, result.getYmax()/deg],
-        aspect='auto')
-    cb = plt.colorbar(im)
-    cb.set_label(r'Intensity (arb. u.)', size=16)
-    plt.xlabel(r'$\phi_f (^{\circ})$', fontsize=16)
-    plt.ylabel(r'$\alpha_f (^{\circ})$', fontsize=16)
-    plt.show()
+    return simulation.getIntensityData()
 
 
 if __name__ == '__main__':
-    simulate()
+    ba.simulateThenPlotOrSave(simulate, ba.standardIntensityPlot)
diff --git a/Examples/python/simulation/ex04_ComplexShapes/CoreShellNanoparticles.py b/Examples/python/simulation/ex04_ComplexShapes/CoreShellNanoparticles.py
index b2bcc020f3b..36d02a2bd88 100644
--- a/Examples/python/simulation/ex04_ComplexShapes/CoreShellNanoparticles.py
+++ b/Examples/python/simulation/ex04_ComplexShapes/CoreShellNanoparticles.py
@@ -1,9 +1,7 @@
 """
 Core shell nanoparticles
 """
-import numpy, sys
-import matplotlib
-from matplotlib import pyplot as plt
+import numpy
 import bornagain as ba
 from bornagain import deg, angstrom, nm
 
@@ -61,21 +59,8 @@ def simulate():
     simulation = get_simulation()
     simulation.setSample(sample)
     simulation.runSimulation()
-    result = simulation.getIntensityData()
-
-    # showing the result
-    im = plt.imshow(
-        result.getArray(),
-        norm=matplotlib.colors.LogNorm(1.0, result.getMaximum()),
-        extent=[result.getXmin()/deg, result.getXmax()/deg,
-                result.getYmin()/deg, result.getYmax()/deg],
-        aspect='auto')
-    cb = plt.colorbar(im)
-    cb.set_label(r'Intensity (arb. u.)', size=16)
-    plt.xlabel(r'$\phi_f (^{\circ})$', fontsize=16)
-    plt.ylabel(r'$\alpha_f (^{\circ})$', fontsize=16)
-    plt.show()
+    return simulation.getIntensityData()
 
 
 if __name__ == '__main__':
-    simulate()
+    ba.simulateThenPlotOrSave(simulate, ba.standardIntensityPlot)
diff --git a/Examples/python/simulation/ex04_ComplexShapes/CustomFormFactor.py b/Examples/python/simulation/ex04_ComplexShapes/CustomFormFactor.py
index 002c945f124..1215ef13275 100644
--- a/Examples/python/simulation/ex04_ComplexShapes/CustomFormFactor.py
+++ b/Examples/python/simulation/ex04_ComplexShapes/CustomFormFactor.py
@@ -1,9 +1,7 @@
 """
 Custom form factor in DWBA.
 """
-import numpy, sys
-import matplotlib
-from matplotlib import pyplot as plt
+import numpy
 import bornagain as ba
 from bornagain import deg, angstrom, nm
 import cmath
@@ -98,21 +96,8 @@ def simulate():
     simulation = get_simulation()
     simulation.setSample(sample)
     simulation.runSimulation()
-    result = simulation.getIntensityData()
-
-    # showing the result
-    im = plt.imshow(
-        result.getArray(),
-        norm=matplotlib.colors.LogNorm(1.0, result.getMaximum()),
-        extent=[result.getXmin()/deg, result.getXmax()/deg,
-                result.getYmin()/deg, result.getYmax()/deg],
-        aspect='auto')
-    cb = plt.colorbar(im)
-    cb.set_label(r'Intensity (arb. u.)', size=16)
-    plt.xlabel(r'$\phi_f (^{\circ})$', fontsize=16)
-    plt.ylabel(r'$\alpha_f (^{\circ})$', fontsize=16)
-    plt.show()
+    return simulation.getIntensityData()
 
 
 if __name__ == '__main__':
-    simulate()
+    ba.simulateThenPlotOrSave(simulate, ba.standardIntensityPlot)
diff --git a/Examples/python/simulation/ex04_ComplexShapes/HexagonalLatticesWithBasis.py b/Examples/python/simulation/ex04_ComplexShapes/HexagonalLatticesWithBasis.py
index e7fc0ef1cbf..c3e0c9c1a3a 100644
--- a/Examples/python/simulation/ex04_ComplexShapes/HexagonalLatticesWithBasis.py
+++ b/Examples/python/simulation/ex04_ComplexShapes/HexagonalLatticesWithBasis.py
@@ -1,9 +1,7 @@
 """
 Spheres on two hexagonal close packed layers
 """
-import numpy, sys
-import matplotlib
-from matplotlib import pyplot as plt
+import numpy
 import bornagain as ba
 from bornagain import deg, angstrom, nm
 
@@ -65,21 +63,8 @@ def simulate():
     simulation = get_simulation()
     simulation.setSample(sample)
     simulation.runSimulation()
-    result = simulation.getIntensityData()
-
-    # showing the result
-    im = plt.imshow(
-        result.getArray(),
-        norm=matplotlib.colors.LogNorm(1.0, result.getMaximum()),
-        extent=[result.getXmin()/deg, result.getXmax()/deg,
-                result.getYmin()/deg, result.getYmax()/deg],
-        aspect='auto')
-    cb = plt.colorbar(im)
-    cb.set_label(r'Intensity (arb. u.)', size=16)
-    plt.xlabel(r'$\phi_f (^{\circ})$', fontsize=16)
-    plt.ylabel(r'$\alpha_f (^{\circ})$', fontsize=16)
-    plt.show()
+    return simulation.getIntensityData()
 
 
 if __name__ == '__main__':
-    simulate()
+    ba.simulateThenPlotOrSave(simulate, ba.standardIntensityPlot)
diff --git a/Examples/python/simulation/ex04_ComplexShapes/LargeParticlesFormFactor.py b/Examples/python/simulation/ex04_ComplexShapes/LargeParticlesFormFactor.py
index bbe28dbaeab..0178f580627 100644
--- a/Examples/python/simulation/ex04_ComplexShapes/LargeParticlesFormFactor.py
+++ b/Examples/python/simulation/ex04_ComplexShapes/LargeParticlesFormFactor.py
@@ -7,8 +7,6 @@ oscillates rapidly within one detector bin and analytical calculations
 In this case Monte-Carlo integration over detector bin should be used.
 """
 import numpy, sys
-import matplotlib
-from matplotlib import pyplot as plt
 import bornagain as ba
 from bornagain import deg, angstrom, nm
 
@@ -57,59 +55,79 @@ def get_simulation(integration_flag):
     return simulation
 
 
-def simulate():
+def simulate(condi):
     """
-    Runs simulation and plots 4 results:
-for small and large cylinders, with and without integration
+    Runs simulation and returns result.
     """
+    scale = condi['scale']
+    integration_flag = condi['integration']
+    sample = get_sample(default_cylinder_radius*scale,
+                        default_cylinder_height*scale)
+    simulation = get_simulation(integration_flag)
+    simulation.setSample(sample)
+    simulation.runSimulation()
+    return simulation.getIntensityData()
+
+def plot(result, nframe, title):
+    plt.subplot(2, 2, nframe+1)
+    plt.subplots_adjust(wspace=0.3, hspace=0.3)
+    im = plt.imshow(
+        result.getArray(),
+        norm=matplotlib.colors.LogNorm(1.0, result.getMaximum()),
+        extent=[result.getXmin()/deg, result.getXmax()/deg,
+                result.getYmin()/deg, result.getYmax()/deg],
+        aspect='auto')
+    cb = plt.colorbar(im)
+    cb.set_label(r'Intensity (arb. u.)', size=16)
+    plt.xlabel(r'$\phi_f (^{\circ})$', fontsize=16)
+    plt.ylabel(r'$\alpha_f (^{\circ})$', fontsize=16)
+    plt.text(0.0, 2.1,
+             title,
+             horizontalalignment='center',
+             verticalalignment='center',
+             fontsize=13)
 
-    fig = plt.figure(figsize=(12.80, 10.24))
+if __name__ == '__main__':
+    """
+    Runs one simulation for each condition, and plots results on a single canvas.
+    Conditions are small and large cylinders, with and without integration.
+    """
+    if len(sys.argv)<2:
+        print( "Usage:" )
+        print( "  " + sys.argv[0] + " -p                           # to plot results" )
+        print( "  " + sys.argv[0] + " <filename without extension> # to save results" )
+        sys.exit()
 
     # conditions to define cylinders scale factor and Monte-Carlo integration flag
     conditions = [
-        {'title': "Small cylinders, analytical calculations", 'scale': 1,
+        {'name': "SmallAn",
+         'title': "Small cylinders, analytical calculations", 'scale': 1,
          'integration': False, 'max': 1e+08},
-        {'title': "Small cylinders, Monte-Carlo integration", 'scale': 1,
+        {'name': "SmallMC",
+         'title': "Small cylinders, Monte-Carlo integration", 'scale': 1,
          'integration': True,  'max': 1e+08},
-        {'title': "Large cylinders, analytical calculations", 'scale': 100,
+        {'name': "LargeAn",
+         'title': "Large cylinders, analytical calculations", 'scale': 100,
          'integration': False, 'max': 1e+12},
-        {'title': "Large cylinders, Monte-Carlo integration", 'scale': 100,
+        {'name': "LargeMC",
+         'title': "Large cylinders, Monte-Carlo integration", 'scale': 100,
          'integration': True,  'max': 1e+12}
     ]
 
-    # run simulation 4 times and plot results
-    for i_plot in range(0, len(conditions)):
-        scale = conditions[i_plot]['scale']
-        integration_flag = conditions[i_plot]['integration']
-
-        sample = get_sample(default_cylinder_radius*scale,
-                            default_cylinder_height*scale)
-        simulation = get_simulation(integration_flag)
-        simulation.setSample(sample)
-        simulation.runSimulation()
-        result = simulation.getIntensityData()
-
-        # plotting results
-        plt.subplot(2, 2, i_plot+1)
-        plt.subplots_adjust(wspace=0.3, hspace=0.3)
-        im = plt.imshow(
-            result.getArray(),
-            norm=matplotlib.colors.LogNorm(1.0, result.getMaximum()),
-            extent=[result.getXmin()/deg, result.getXmax()/deg,
-                    result.getYmin()/deg, result.getYmax()/deg],
-            aspect='auto')
-        cb = plt.colorbar(im)
-        cb.set_label(r'Intensity (arb. u.)', size=16)
-        plt.xlabel(r'$\phi_f (^{\circ})$', fontsize=16)
-        plt.ylabel(r'$\alpha_f (^{\circ})$', fontsize=16)
-        plt.text(0.0, 2.1,
-                 conditions[i_plot]['title'],
-                 horizontalalignment='center',
-                 verticalalignment='center',
-                 fontsize=13)
-
-    plt.show()
-
-
-if __name__ == '__main__':
-    simulate()
+    if sys.argv[1] == "-p":
+        import matplotlib
+        from matplotlib import pyplot as plt
+        plt.figure(figsize=(12.80, 10.24))
+        for nplot in range(len(conditions)):
+            condi = conditions[nplot]
+            title = condi['title']
+            print("Generating intensity map for " + title)
+            intensities = simulate(condi)
+            plot(intensities, nplot, title)
+        plt.show()
+    else:
+        for condi in conditions:
+            intensities = simulate(condi)
+            fname = "%s.%s.int" % (sys.argv[1], condi['name'])
+            ba.IntensityDataIOFactory.writeIntensityData(intensities, fname)
+            print("Stored intensity map in " + fname)
diff --git a/Examples/python/simulation/ex05_BeamAndDetector/BeamDivergence.py b/Examples/python/simulation/ex05_BeamAndDetector/BeamDivergence.py
index 84ba9f54ce9..635ccd38e36 100644
--- a/Examples/python/simulation/ex05_BeamAndDetector/BeamDivergence.py
+++ b/Examples/python/simulation/ex05_BeamAndDetector/BeamDivergence.py
@@ -1,7 +1,7 @@
 """
 Cylinder form factor in DWBA with beam divergence
 """
-import numpy, sys
+import numpy
 import matplotlib
 from matplotlib import pyplot as plt
 import bornagain as ba
@@ -64,7 +64,7 @@ def simulate():
     simulation.printParameters()
     simulation.runSimulation()
 
-    result = simulation.getIntensityData()
+    return simulation.getIntensityData()
 
     # showing the result
     im = plt.imshow(
@@ -81,4 +81,4 @@ def simulate():
 
 
 if __name__ == '__main__':
-    simulate()
+    ba.simulateThenPlotOrSave(simulate, ba.standardIntensityPlot)
diff --git a/Examples/python/simulation/ex05_BeamAndDetector/DetectorResolutionFunction.py b/Examples/python/simulation/ex05_BeamAndDetector/DetectorResolutionFunction.py
index 44167c48fb4..a21f5857154 100644
--- a/Examples/python/simulation/ex05_BeamAndDetector/DetectorResolutionFunction.py
+++ b/Examples/python/simulation/ex05_BeamAndDetector/DetectorResolutionFunction.py
@@ -1,7 +1,7 @@
 """
 Cylinder form factor in DWBA with detector resolution function applied
 """
-import numpy, sys
+import numpy
 import matplotlib
 from matplotlib import pyplot as plt
 import bornagain as ba
@@ -59,7 +59,7 @@ def simulate():
     simulation.setSample(sample)
     simulation.printParameters()
     simulation.runSimulation()
-    result = simulation.getIntensityData()
+    return simulation.getIntensityData()
 
     # showing the result
     im = plt.imshow(
@@ -76,4 +76,4 @@ def simulate():
 
 
 if __name__ == '__main__':
-    simulate()
+    ba.simulateThenPlotOrSave(simulate, ba.standardIntensityPlot)
diff --git a/Examples/python/simulation/ex05_BeamAndDetector/OffSpecularSimulation.py b/Examples/python/simulation/ex05_BeamAndDetector/OffSpecularSimulation.py
index a6d2c4073f4..5a3745a9bb8 100644
--- a/Examples/python/simulation/ex05_BeamAndDetector/OffSpecularSimulation.py
+++ b/Examples/python/simulation/ex05_BeamAndDetector/OffSpecularSimulation.py
@@ -1,7 +1,7 @@
 """
 Long boxes at 1D lattice, ba.OffSpecular simulation
 """
-import numpy, sys
+import numpy
 import matplotlib
 from matplotlib import pyplot as plt
 import bornagain as ba
@@ -72,7 +72,7 @@ def simulate():
     simulation = get_simulation()
     simulation.setSample(sample)
     simulation.runSimulation()
-    result = simulation.getIntensityData()
+    return simulation.getIntensityData()
 
     # showing the result
     im = plt.imshow(
@@ -89,4 +89,4 @@ def simulate():
 
 
 if __name__ == '__main__':
-    simulate()
+    ba.simulateThenPlotOrSave(simulate, ba.standardIntensityPlot)
diff --git a/Examples/python/simulation/ex05_BeamAndDetector/RectangularDetector.py b/Examples/python/simulation/ex05_BeamAndDetector/RectangularDetector.py
index 218421abc06..1e8fb2786cd 100644
--- a/Examples/python/simulation/ex05_BeamAndDetector/RectangularDetector.py
+++ b/Examples/python/simulation/ex05_BeamAndDetector/RectangularDetector.py
@@ -2,7 +2,7 @@
 Simulation with rectangular detector. Pilatus3-1M detector is used as an example.
 Results will be compared against simulation with spherical detector.
 """
-import numpy, sys
+import numpy
 import matplotlib
 from matplotlib import pyplot as plt
 import bornagain as ba
@@ -146,4 +146,4 @@ def simulate():
     plot_results(result_sph, result_rect)
 
 if __name__ == '__main__':
-    simulate()
+    ba.simulateThenPlotOrSave(simulate, ba.standardIntensityPlot)
diff --git a/Examples/python/simulation/ex05_BeamAndDetector/SpecularSimulation.py b/Examples/python/simulation/ex05_BeamAndDetector/SpecularSimulation.py
index 577104f53b4..76e2446b35a 100644
--- a/Examples/python/simulation/ex05_BeamAndDetector/SpecularSimulation.py
+++ b/Examples/python/simulation/ex05_BeamAndDetector/SpecularSimulation.py
@@ -1,7 +1,7 @@
 """
 R and T coefficients in multilayer, ba.Specular simulation.
 """
-import numpy, sys
+import numpy
 import matplotlib
 from matplotlib import pyplot as plt
 import bornagain as ba
@@ -93,4 +93,4 @@ def simulate():
 
 
 if __name__ == '__main__':
-    simulate()
+    ba.simulateThenPlotOrSave(simulate, ba.standardIntensityPlot)
diff --git a/Examples/python/simulation/ex06_Miscellaneous/AccessingSimulationResults.py b/Examples/python/simulation/ex06_Miscellaneous/AccessingSimulationResults.py
index 6f3aabc8e0c..7cce1159543 100644
--- a/Examples/python/simulation/ex06_Miscellaneous/AccessingSimulationResults.py
+++ b/Examples/python/simulation/ex06_Miscellaneous/AccessingSimulationResults.py
@@ -3,7 +3,7 @@ Extended example for simulation results treatment (cropping, slicing, exporting)
 The standard "Cylinders in DWBA" sample is used to setup the simulation.
 """
 import math
-import numpy, sys
+import numpy
 import matplotlib
 import random
 from matplotlib import pyplot as plt
@@ -186,10 +186,10 @@ def simulate():
     simulation = get_simulation()
     simulation.setSample(sample)
     simulation.runSimulation()
-    result = simulation.getIntensityData()
+    return simulation.getIntensityData()
 
     plot_results(result)
 
 
 if __name__ == '__main__':
-    simulate()
+    ba.simulateThenPlotOrSave(simulate, ba.standardIntensityPlot)
diff --git a/Examples/python/simulation/ex06_Miscellaneous/AxesInDifferentUnits.py b/Examples/python/simulation/ex06_Miscellaneous/AxesInDifferentUnits.py
index b9d6ab9036d..dd7bca4e5f6 100644
--- a/Examples/python/simulation/ex06_Miscellaneous/AxesInDifferentUnits.py
+++ b/Examples/python/simulation/ex06_Miscellaneous/AxesInDifferentUnits.py
@@ -2,7 +2,7 @@
 In this example we demonstrate how to plot simulation results with
 axes in different units (nbins, mm, degs and QyQz).
 """
-import numpy, sys
+import numpy
 import matplotlib
 from matplotlib import pyplot as plt
 import bornagain as ba
@@ -89,22 +89,22 @@ def simulate():
 
     plt.subplot(2, 2, 1)
     # default units for rectangular detector are millimeters
-    result = simulation.getIntensityData()
+    return simulation.getIntensityData()
     plot_as_colormap.esult, "In default units",
                      r'$X_{mm}$', r'$Y_{mm}$')
 
     plt.subplot(2, 2, 2)
-    result = simulation.getIntensityData(ba.IDetector2D.NBINS)
+    return simulation.getIntensityData(ba.IDetector2D.NBINS)
     plot_as_colormap.esult, "In number of bins",
                      r'$X_{nbins}$', r'$Y_{nbins}$')
 
     plt.subplot(2, 2, 3)
-    result = simulation.getIntensityData(ba.IDetector2D.DEGREES)
+    return simulation.getIntensityData(ba.IDetector2D.DEGREES)
     plot_as_colormap.esult, "In degs",
                      r'$\phi_f ^{\circ}$', r'$\alpha_f ^{\circ}$')
 
     plt.subplot(2, 2, 4)
-    result = simulation.getIntensityData(ba.IDetector2D.QYQZ)
+    return simulation.getIntensityData(ba.IDetector2D.QYQZ)
     plot_as_colormap.esult, "Q-space",
                      r'$Q_{y} [1/nm]$', r'$Q_{z} [1/nm]$')
 
@@ -113,4 +113,4 @@ def simulate():
 
 
 if __name__ == '__main__':
-    simulate()
+    ba.simulateThenPlotOrSave(simulate, ba.standardIntensityPlot)
diff --git a/Examples/python/utils/plot_intensity_data.py b/Examples/python/utils/plot_intensity_data.py
index 9b84198af96..9acbf973d7c 100755
--- a/Examples/python/utils/plot_intensity_data.py
+++ b/Examples/python/utils/plot_intensity_data.py
@@ -2,7 +2,7 @@
 # Plots intensity data stored in BornAgain "*.int" or "*.int.gz" format
 # Usage: python plot_intensity_data.py intensity_file.int.gz
 
-import matplotlib, numpy, sys
+import matplotlib, numpy
 from matplotlib import pyplot as plt
 import bornagain as ba
 from bornagain import deg, angstrom, nm
diff --git a/Examples/python/utils/plot_intensity_data_diff.py b/Examples/python/utils/plot_intensity_data_diff.py
index be8d01e9479..f6ad420ab79 100755
--- a/Examples/python/utils/plot_intensity_data_diff.py
+++ b/Examples/python/utils/plot_intensity_data_diff.py
@@ -2,7 +2,7 @@
 # Plots intensity data difference stored in BornAgain "*.int" or "*.int.gz" format
 # Usage: python plot_intensity_data.py intensity_reference.int.gz intensity_other.int.gz
 
-import matplotlib, numpy, sys
+import matplotlib, numpy
 from matplotlib import pyplot as plt
 import bornagain as ba
 from bornagain import deg, angstrom, nm
diff --git a/Examples/python/utils/show2d_root.py b/Examples/python/utils/show2d_root.py
index 3a74514b4bb..cb3d1dbaf05 100755
--- a/Examples/python/utils/show2d_root.py
+++ b/Examples/python/utils/show2d_root.py
@@ -5,7 +5,7 @@
 from __future__ import print_function
 import argparse
 import os
-import numpy, sys
+import numpy
 import ROOT
 from pylab import *
 
diff --git a/Wrap/swig/CorePython.py b/Wrap/swig/CorePython.py
index 799ce9edba9..147840503cc 100644
--- a/Wrap/swig/CorePython.py
+++ b/Wrap/swig/CorePython.py
@@ -10,13 +10,13 @@
 #
 #   @homepage  http://apps.jcns.fz-juelich.de/BornAgain
 #   @license   GNU General Public License v3 or higher (see COPYING)
-#   @copyright Forschungszentrum Jülich GmbH 2016
+#   @copyright Forschungszentrum Juelich GmbH 2016
 #   @authors   Scientific Computing Group at MLZ Garching
 #   @authors   J. Fisher, M. Ganeva, G. Pospelov, W. Van Herck, J. Wuttke
 #
 #  **************************************************************************  #
 
-def simulate_then_plot_or_save(simulate, plot):
+def simulateThenPlotOrSave(simulate, plot):
     """
     Runs a simulation. Then plots the function or saves the result, depending on given argument.
     """
@@ -32,4 +32,22 @@ def simulate_then_plot_or_save(simulate, plot):
     else:
         plot(result)
 
+def standardIntensityPlot(result):
+    """
+    Plots intensity map.
+    """
+    import matplotlib
+    from matplotlib import pyplot as plt
+    im = plt.imshow(
+        result.getArray(),
+        norm=matplotlib.colors.LogNorm(1.0, result.getMaximum()),
+        extent=[result.getXmin()/deg, result.getXmax()/deg,
+                result.getYmin()/deg, result.getYmax()/deg],
+        aspect='auto')
+    cb = plt.colorbar(im)
+    cb.set_label(r'Intensity (arb. u.)', size=16)
+    plt.xlabel(r'$\phi_f (^{\circ})$', fontsize=16)
+    plt.ylabel(r'$\alpha_f (^{\circ})$', fontsize=16)
+    plt.show()
+
 %}
diff --git a/auto/Wrap/libBornAgainCore.py b/auto/Wrap/libBornAgainCore.py
index e6599ea2d73..7d3a6aee6d6 100644
--- a/auto/Wrap/libBornAgainCore.py
+++ b/auto/Wrap/libBornAgainCore.py
@@ -24529,7 +24529,23 @@ SimulationFactory_swigregister(SimulationFactory)
 
 
 
-def simulate_then_plot_or_save(simulate, plot):
+#  **************************************************************************  #
+#
+#   BornAgain: simulate and fit scattering at grazing incidence
+#
+#   @file      Wrap/swig/CorePython.py
+#   @brief     Python extensions of the SWIG-genrated Python module bornagain.
+#              This file is included by libBornAgainCore.i.
+#
+#   @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
+#
+#  **************************************************************************  #
+
+def simulateThenPlotOrSave(simulate, plot):
     """
     Runs a simulation. Then plots the function or saves the result, depending on given argument.
     """
@@ -24545,6 +24561,24 @@ def simulate_then_plot_or_save(simulate, plot):
     else:
         plot(result)
 
+def standardIntensityPlot(result):
+    """
+    Plots intensity map.
+    """
+    import matplotlib
+    from matplotlib import pyplot as plt
+    im = plt.imshow(
+        result.getArray(),
+        norm=matplotlib.colors.LogNorm(1.0, result.getMaximum()),
+        extent=[result.getXmin()/deg, result.getXmax()/deg,
+                result.getYmin()/deg, result.getYmax()/deg],
+        aspect='auto')
+    cb = plt.colorbar(im)
+    cb.set_label(r'Intensity (arb. u.)', size=16)
+    plt.xlabel(r'$\phi_f (^{\circ})$', fontsize=16)
+    plt.ylabel(r'$\alpha_f (^{\circ})$', fontsize=16)
+    plt.show()
+
 
 # This file is compatible with both classic and new-style classes.
 
-- 
GitLab