diff --git a/Core/CMakeLists.txt b/Core/CMakeLists.txt
index a1d48e8828a730e7c63ea1b6b4144e2408d9c33c..848472e5e2aad555f9c8805026304a2b3480fe1d 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 6b485d549a04c5373f8f24818e99ea01cf4e57b9..475038acfaf1947d76b6818fa1fcf2ac510839b9 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 1bd6e7943fd50d72838845f143663f9d1e7c2719..cde1ee93b72cfc17e3dae00dc1359e35e5a2f534 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 2686f86b8865422898cd2f7067d4d456d7dece01..6a7c182442e4be94410c195450a2d07b661528e8 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 cefa41da8228200471ad68b2b1ecc9516204b776..95c47f396c5b594f1a47dcf47ff21c66ff5fd441 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 5c1335deab389eabe98b4073469c5a956f72a3c1..df65592d9e02437942ab198499bdc002641716f0 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 a069b24920a40affc70b712fc5757e85eeb675b6..06387541b0da5a216b3c9ad187f0aaec9308982b 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 d63255067708829e76b2f0c9623ca626a4cad4bc..ac2849c08d371c04faef803a8bc93c2a617af9c0 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 6c2ca1e52cc6cfcb1e9e9163813124cd552e54b9..f3254acd1192f81981b5c9f2eaec1cd6c5e15c78 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 00cc0cbeba49f3fddd4d556595c8963ca7838d46..a0af5935bfdb8b5c2519b9036f1471c60f0016bd 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 10cd4b2c87bf3542a6a99910ada253e89423823d..a2a8fd5d892f5e5c2da4eb570c88919f52f352d4 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 fd549e6be02f14fa080e508247b5504d4d548f41..9b5d61cc93ee94bb8ce4cce93881a388c6801571 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 8c2779fec36afac01533cce74b9c680132ad091d..3980879060f7c5a0e3578283e8ed1e042c9f680a 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 784826769044d3c69ff07870cfdd72ce656e37c0..5b248db688798fd8253d320a375430594bede62d 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 aa688431c4febf2ca3c818141cfddc58de6048f8..4c78b21b133050a67310e220fe06404dcf881033 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 cd2ed9e273a98427b5abce35e08e54177c8e3aaf..78da8ef6e9f3d01dc12b950a8cfe7358fe6c885d 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 55f72df30ac41e221259148745f60225e50db9e6..b3359474c1a108fd0df8e8651dc108022320084c 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 332e18c67829001629017b11db2b6317469b63cb..a2f8db7a5c2e2d758976496f7b275fcd2f8a7b49 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 6dfcd3050f17a2589a1664bcdd7733137337acfb..bbfc555456df26d72d999c329538a00dd28eeea8 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 4b0507c2cf556029a0f2579433838ba37dd1ae2c..074eb75f171480a6c1be32e7246075de6812352a 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 ed17848a0975de317e422680316a9bcf0172cc37..38cbd8e4fbcc6dd6f375b191c1cef470c426d182 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 6268332b237ee4ef8cdf59855d409ca0d86716b4..7da27bb23a9f25bb1fe8a87435f50674dab9f2a2 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 5910e22f015a02126c1bfc55cc090fb888661369..92c9b115ee1be3f328bdfd36eb860f9b1ea26a81 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 a0486576c2e84d18004013084e6862825d83a25d..d3b6aef158d6fb9296719159af98d6e501d256c2 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 93f99af6e34fbe496c49c61d2086fb7affad0e2f..1fede9b2f0b0eabe8b8ff1a8f6add77f38290e4d 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 6dfdb24970264391faf89506efe2f3cabaffee0b..eb41a7cdb8bb1e24c24c1138a0c8746c45d62504 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 b2bcc020f3bc36e5ad480f6ec5db73c4ba38754b..36d02a2bd88b705e75dec3731664b1efd18d8370 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 002c945f12491f5d15ddcd4c6e2fa74b875b30f6..1215ef132755f4bfd2127b04f064fa03545cb4e8 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 e7fc0ef1cbf9fad3170554d4e00939142f3a995e..c3e0c9c1a3a867525dab1f87cd9895ce0645da82 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 bbe28dbaeab2367d239e157f001b11582447707c..0178f580627a1be478678523eb2f79ceb1418930 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 84ba9f54ce902cdd7cbfe659d3b7853b52f1bba7..635ccd38e360be10123676535c8ef3de5b71fecd 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 44167c48fb40e3c249816955479c6f11aac22c03..a21f5857154f72bdd85b3c2bf5c45a0d7462f1ac 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 a6d2c4073f4838e6a3bfbea7d3b4e9e545abcf0e..5a3745a9bb86b4f8a20b36f8060e8bc944bc3c2b 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 218421abc0645bbcde2d3ae65ac11c73d6ff5a57..1e8fb2786cdadcd5370b5b644d7b6e61a3aa62b2 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 577104f53b4d09c1c34f848bc12b6ec83cd56a7d..76e2446b35af610d54b657ee9cbeb8b6ee862473 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 6f3aabc8e0c62f60737ef4515ba0d62eabcf75c1..7cce11595430c96e4e18e06607a27ac6d36a9b69 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 b9d6ab9036d0d5abca3bc59c5d47fa0985b23f0f..dd7bca4e5f6587792a571ebe8a5d531d26df7e01 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 9b84198af961abab6857509592f585d2404e7642..9acbf973d7c274a5294229ba130089adf51f3408 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 be8d01e9479fe175eb9d32fb37337fc672cab22c..f6ad420ab79e2546d8d1a75f6d2cce8cbd8ae1ae 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 3a74514b4bbd6726391229b879ecb579a01cf46a..cb3d1dbaf05247e5a85046f53e9268c1d8bc44cf 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 799ce9edba99b754b9843e9d393d99270355131b..147840503cc055a0c1a5927f981bc7953d24ac18 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 e6599ea2d73d8bbecf5b41b31d395a634aaa6cf0..7d3a6aee6d6261df30dfddf28b034163a128d3e2 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.