diff --git a/Examples/Demos/animation.mpg b/Examples/Demos/animation.mpg
new file mode 100644
index 0000000000000000000000000000000000000000..b3f73b1c9839585c66cd0d6ab2b8443116d9d8ef
Binary files /dev/null and b/Examples/Demos/animation.mpg differ
diff --git a/Examples/Demos/animation2.mpg b/Examples/Demos/animation2.mpg
new file mode 100644
index 0000000000000000000000000000000000000000..8b9c25094c68e80081da9fcf52d001742e9dfd6c
Binary files /dev/null and b/Examples/Demos/animation2.mpg differ
diff --git a/Examples/Demos/simul_demo_cyl_SSCA.py b/Examples/Demos/simul_demo_cyl_SSCA.py
new file mode 100644
index 0000000000000000000000000000000000000000..3b4cba846bc68c8bb1c3e34bfa8e22b683a522e0
--- /dev/null
+++ b/Examples/Demos/simul_demo_cyl_SSCA.py
@@ -0,0 +1,73 @@
+'''
+Simulation demo: Size Space Coupling Approximation
+'''
+
+import numpy, pylab, matplotlib
+
+from libBornAgainCore import *
+
+# ----------------------------------
+# describe sample and run simulation
+# ----------------------------------
+def RunSimulation():
+    # defining materials
+    mAir = MaterialManager.getHomogeneousMaterial("Air", 0.0, 0.0)
+    mSubstrate = MaterialManager.getHomogeneousMaterial("Substrate", 6e-6, 2e-8)
+    mParticle = MaterialManager.getHomogeneousMaterial("Particle", 6e-4, 2e-8)
+    mLayer = MaterialManager.getHomogeneousMaterial("Layer", 2e-5, 2e-8)
+
+    # collection of particles
+    cylinder_ff1 = FormFactorCylinder(5 * nanometer, 2 * nanometer)
+    cylinder_ff2 = FormFactorCylinder(6 * nanometer, 3 * nanometer)
+    cylinder_ff3 = FormFactorCylinder(7 * nanometer, 4 * nanometer)
+    cylinder1 = Particle(mParticle, cylinder_ff1)
+    cylinder2 = Particle(mParticle, cylinder_ff2)
+    cylinder3 = Particle(mParticle, cylinder_ff3)
+    particle_decoration = ParticleDecoration()
+    particle_decoration.addParticle(cylinder1)
+    particle_decoration.addParticle(cylinder2)
+    particle_decoration.addParticle(cylinder3)
+    interference = InterferenceFunction1DParaCrystal(5 * nanometer, 1 * nanometer)
+    # set coupling between size and space
+    interference.setKappa(4)
+    particle_decoration.addInterferenceFunction(interference)
+
+    # air layer with particles and substrate form multi layer
+    air_layer = Layer(mAir)
+    air_layer.setDecoration(particle_decoration)
+    film_layer = Layer(mLayer, 5 * nanometer)
+    substrate_layer = Layer(mSubstrate)
+    subs2_layer = Layer(mSubstrate, 5 * nanometer)
+    multi_layer = MultiLayer()
+    multi_layer.addLayer(air_layer)
+#    multi_layer.addLayer(film_layer)
+    roughness = LayerRoughness(10 * nanometer, 3, 20 * nanometer)
+    multi_layer.addLayerWithTopRoughness(substrate_layer, roughness)
+
+    # simulation parameters
+    sim_params = SimulationParameters()
+#    sim_params.me_if_approx = SimulationParameters.SSCA
+    sim_params.me_if_approx = SimulationParameters.DA
+    # build and run experiment
+    simulation = Simulation()
+    simulation.setSimulationParameters(sim_params)
+    simulation.setDetectorParameters(100, -4.0 * degree, 4.0 * degree, 100, 0.0 * degree, 8.0 * degree)
+    simulation.setBeamParameters(1.0 * angstrom, 0.2 * degree, 0.0 * degree)
+    simulation.setSample(multi_layer)
+    simulation.runSimulation()
+    # intensity data
+    return GetOutputData(simulation)
+
+
+#-------------------------------------------------------------
+# main()
+#-------------------------------------------------------------
+if __name__ == '__main__':
+    result = RunSimulation() + 1 # for log scale
+    im = pylab.imshow(numpy.rot90(result, 1),
+                 norm=matplotlib.colors.LogNorm(),
+                 extent=[-4.0, 4.0, 0, 8.0])
+    pylab.colorbar(im)
+    pylab.xlabel(r'$\phi_f$', fontsize=20)
+    pylab.ylabel(r'$\alpha_f$', fontsize=20)
+    pylab.show()
diff --git a/Examples/Demos/simul_demo_cylsphere.py b/Examples/Demos/simul_demo_cylsphere.py
new file mode 100644
index 0000000000000000000000000000000000000000..5472ae0d3a59d6d2eac312668b4e83c47acce558
--- /dev/null
+++ b/Examples/Demos/simul_demo_cylsphere.py
@@ -0,0 +1,58 @@
+'''
+Simulation demo: Cylinder and/or sphere on substrate
+'''
+
+import numpy, pylab, matplotlib
+
+from libBornAgainCore import *
+
+# ----------------------------------
+# describe sample and run simulation
+# ----------------------------------
+def RunSimulation():
+    # defining materials
+    mAir = MaterialManager.getHomogeneousMaterial("Air", 0.0, 0.0)
+    mSubstrate = MaterialManager.getHomogeneousMaterial("Substrate", 6e-6, 2e-8)
+    mParticle = MaterialManager.getHomogeneousMaterial("Particle", 6e-4, 2e-8)
+
+    # collection of particles
+    cylinder_ff = FormFactorCylinder(5 * nanometer, 2 * nanometer)
+    cylinder = Particle(mParticle, cylinder_ff)
+#    sphere_ff = FormFactorFullSphere(4 * nanometer)
+#    sphere = Particle(mParticle, sphere_ff)
+    particle_decoration = ParticleDecoration()
+    particle_decoration.addParticle(cylinder)
+#    particle_decoration.addParticle(sphere)
+#    interference = InterferenceFunction1DParaCrystal(20 * nanometer, 2 * nanometer)
+#    particle_decoration.addInterferenceFunction(interference)
+
+    # air layer with particles and substrate form multi layer
+    air_layer = Layer(mAir)
+    air_layer.setDecoration(particle_decoration)
+    substrate_layer = Layer(mSubstrate)
+    multi_layer = MultiLayer()
+    multi_layer.addLayer(air_layer)
+    multi_layer.addLayer(substrate_layer)
+
+    # build and run experiment
+    simulation = Simulation()
+    simulation.setDetectorParameters(100, -4.0 * degree, 4.0 * degree, 100, 0.0 * degree, 8.0 * degree)
+    simulation.setBeamParameters(1.0 * angstrom, 0.2 * degree, 0.0 * degree)
+    simulation.setSample(multi_layer)
+    simulation.runSimulation()
+    # intensity data
+    return GetOutputData(simulation)
+
+
+#-------------------------------------------------------------
+# main()
+#-------------------------------------------------------------
+if __name__ == '__main__':
+    result = RunSimulation() + 1 # for log scale
+    im = pylab.imshow(numpy.rot90(result, 1),
+                 norm=matplotlib.colors.LogNorm(),
+                 extent=[-4.0, 4.0, 0, 8.0])
+    pylab.colorbar(im)
+    pylab.xlabel(r'$\phi_f$', fontsize=20)
+    pylab.ylabel(r'$\alpha_f$', fontsize=20)
+    pylab.show()
diff --git a/Examples/Demos/simul_demo_lattice1.py b/Examples/Demos/simul_demo_lattice1.py
new file mode 100644
index 0000000000000000000000000000000000000000..343f4cc6b940a82103ee1d559e8bfa510bfd32d2
--- /dev/null
+++ b/Examples/Demos/simul_demo_lattice1.py
@@ -0,0 +1,92 @@
+'''
+Simulation demo: 2d lattice structure of particles
+'''
+
+import numpy, pylab, matplotlib
+
+from libBornAgainCore import *
+
+M_PI = numpy.pi
+
+# ----------------------------------
+# describe sample and run simulation
+# ----------------------------------
+def RunSimulation():
+    # defining materials
+    mAmbience = MaterialManager.getHomogeneousMaterial("Air", 0.0, 0.0 )
+    mSubstrate = MaterialManager.getHomogeneousMaterial("Substrate", 6e-6, 2e-8 )
+    mParticle = MaterialManager.getHomogeneousMaterial("Particle", 6e-4, 2e-8 )
+    
+    # particle
+    cylinder_ff = FormFactorCylinder(5*nanometer, 5*nanometer)
+    cylinder = Particle(mParticle, cylinder_ff.clone())
+    position = kvector_t(0.0, 0.0, 0.0)
+    particle_info =  PositionParticleInfo(cylinder, position, 1.0)
+    particle_decoration = ParticleDecoration()
+    particle_decoration.addParticleInfo(particle_info)
+
+    # set lattice parameters
+    lattice_params = Lattice2DIFParameters()
+    lattice_params.m_length_1 = 10.0*nanometer
+    lattice_params.m_length_2 = 10.0*nanometer
+    lattice_params.m_angle = 90.0*degree
+    lattice_params.m_xi = 0.0*degree
+    lattice_params.m_domain_size_1 = 20000.0*nanometer
+    lattice_params.m_domain_size_2 = 20000.0*nanometer
+    lattice_params.m_corr_length_1 = 300.0*nanometer/2.0/M_PI
+    lattice_params.m_corr_length_2 = 100.0*nanometer/2.0/M_PI
+
+    # interference function
+    interference = InterferenceFunction2DLattice(lattice_params)
+    pdf = FTDistribution2DCauchy(300.0*nanometer/2.0/M_PI, 100.0*nanometer/2.0/M_PI)
+    interference.setProbabilityDistribution(pdf)
+    particle_decoration.addInterferenceFunction(interference)
+
+    # top air layer
+    air_layer = Layer(mAmbience)
+    air_layer.setDecoration(particle_decoration)
+
+    # substrate layer
+    substrate_layer = Layer(mSubstrate, 0)
+
+    # multilayer
+    multi_layer = MultiLayer()
+    multi_layer.addLayer(air_layer)
+    multi_layer.addLayer(substrate_layer)
+
+    # build and run experiment
+    simulation = Simulation()
+    simulation.setDetectorParameters(100, -2.0*degree, 2.0*degree, 100, 0.0*degree, 4.0*degree)
+    simulation.setBeamParameters(1.0*angstrom, 0.2*degree, 0.0*degree)
+
+    # simulation parameters
+    sim_params= SimulationParameters()
+    sim_params.me_framework = SimulationParameters.DWBA
+    sim_params.me_if_approx = SimulationParameters.LMA
+    sim_params.me_lattice_type = SimulationParameters.LATTICE
+    simulation.setSimulationParameters(sim_params)
+
+    # run simulation
+    simulation.setSample(multi_layer)
+    simulation.runSimulation()
+    return GetOutputData(simulation)
+
+
+#-------------------------------------------------------------
+# main()
+#-------------------------------------------------------------
+if __name__ == '__main__':
+    result = RunSimulation()
+    im = pylab.imshow(numpy.rot90(result+1, 1),
+                 norm=matplotlib.colors.LogNorm(),
+                 extent=[-2.0, 2.0, 0, 4.0])
+    pylab.colorbar(im)
+    pylab.xlabel(r'$\phi_f$', fontsize=20)
+    pylab.ylabel(r'$\alpha_f$', fontsize=20)
+    pylab.show()
+
+
+
+
+ 
+ 
diff --git a/Examples/Demos/simul_demo_lattice2.py b/Examples/Demos/simul_demo_lattice2.py
new file mode 100644
index 0000000000000000000000000000000000000000..bd6ca22c05a02e656d6b295841bd543317cb97d7
--- /dev/null
+++ b/Examples/Demos/simul_demo_lattice2.py
@@ -0,0 +1,94 @@
+'''
+Simulation demo: Cylinder form factor without interference
+'''
+
+import numpy, pylab, matplotlib
+
+from libBornAgainCore import *
+
+M_PI = numpy.pi
+
+# ----------------------------------
+# describe sample and run simulation 
+# ----------------------------------
+def RunSimulation():
+    # defining materials
+    mAmbience = MaterialManager.getHomogeneousMaterial("Air", 0.0, 0.0 )
+    mSubstrate = MaterialManager.getHomogeneousMaterial("Substrate", 6e-6, 2e-8 )
+    mParticle = MaterialManager.getHomogeneousMaterial("Particle", 6e-4, 2e-8 )
+
+    # particle 1
+    cylinder_ff = FormFactorCylinder(5*nanometer, 5*nanometer)
+    cylinder = Particle(mParticle, cylinder_ff)
+    position = kvector_t(0.0, 0.0, 0.0)
+    particle_info = PositionParticleInfo(cylinder, position, 1.0)
+    particle_decoration = ParticleDecoration()
+    particle_decoration.addParticleInfo(particle_info)
+    # particle 2
+    position_2 = kvector_t(5.0*nanometer, 5.0*nanometer, 0.0)
+    particle_info.setPosition(position_2)
+    particle_decoration.addParticleInfo(particle_info)
+
+    # set lattice parameters
+    lattice_params = Lattice2DIFParameters()
+    lattice_params.m_length_1 = 10.0*nanometer
+    lattice_params.m_length_2 = 10.0*nanometer
+    lattice_params.m_angle = 90.0*degree
+    lattice_params.m_xi = 0.0*degree
+    lattice_params.m_domain_size_1 = 20000.0*nanometer
+    lattice_params.m_domain_size_2 = 20000.0*nanometer
+    lattice_params.m_corr_length_1 = 300.0*nanometer/2.0/M_PI
+    lattice_params.m_corr_length_2 = 100.0*nanometer/2.0/M_PI
+
+    # interference function
+    interference = InterferenceFunction2DLattice(lattice_params)
+    pdf = FTDistribution2DCauchy(300.0*nanometer/2.0/M_PI, 100.0*nanometer/2.0/M_PI)
+    interference.setProbabilityDistribution(pdf)
+    particle_decoration.addInterferenceFunction(interference)
+
+    # top air layer
+    air_layer = Layer(mAmbience)
+    air_layer.setDecoration(particle_decoration)
+
+    # substrate layer
+    substrate_layer = Layer(mSubstrate, 0)
+
+    # multilayer
+    multi_layer = MultiLayer()
+    multi_layer.addLayer(air_layer)
+    multi_layer.addLayer(substrate_layer)
+
+    # build and run experiment
+    simulation = Simulation()
+    simulation.setDetectorParameters(100, -2.0*degree, 2.0*degree, 100, 0.0*degree, 4.0*degree, True)
+    simulation.setBeamParameters(1.0*angstrom, 0.2*degree, 0.0*degree)
+
+    # simulation parameters
+    sim_params= SimulationParameters()
+    sim_params.me_framework = SimulationParameters.DWBA
+    sim_params.me_if_approx = SimulationParameters.LMA
+    sim_params.me_lattice_type = SimulationParameters.LATTICE
+    simulation.setSimulationParameters(sim_params)
+
+    # run simulation
+    simulation.setSample(multi_layer)
+    simulation.runSimulation()
+    return GetOutputData(simulation)
+
+
+#-------------------------------------------------------------
+# main()
+#-------------------------------------------------------------
+if __name__ == '__main__':
+    result = RunSimulation()
+    im = pylab.imshow(numpy.rot90(result+1, 1),
+                 norm=matplotlib.colors.LogNorm(),
+                 extent=[-2.0, 2.0, 0, 4.0])
+    pylab.colorbar(im)
+    pylab.xlabel(r'$\phi_f$', fontsize=20)
+    pylab.ylabel(r'$\alpha_f$', fontsize=20)
+    pylab.show()
+
+
+ 
+ 
diff --git a/Examples/Demos/simul_demo_movie.py b/Examples/Demos/simul_demo_movie.py
new file mode 100644
index 0000000000000000000000000000000000000000..6db316ebb3b9f6e975b62c5ca8243f227c9fd4e8
--- /dev/null
+++ b/Examples/Demos/simul_demo_movie.py
@@ -0,0 +1,83 @@
+'''
+Simulation demo: movie of growing particles on substrate
+'''
+
+import os, sys, numpy, pylab, matplotlib
+
+from libBornAgainCore import *
+
+Nframes = 50
+
+radius = 1
+height = 4
+distance = 5
+
+# ----------------------------------
+# describe sample and run simulation
+# ----------------------------------
+def RunSimulation():
+    # defining materials
+    mAir = MaterialManager.getHomogeneousMaterial("Air", 0.0, 0.0)
+    mSubstrate = MaterialManager.getHomogeneousMaterial("Substrate", 6e-6, 2e-8)
+    mParticle = MaterialManager.getHomogeneousMaterial("Particle", 6e-4, 2e-8)
+
+    # collection of particles
+    cylinder_ff = FormFactorCylinder(height, radius)
+    cylinder = Particle(mParticle, cylinder_ff)
+    particle_decoration = ParticleDecoration()
+    particle_decoration.addParticle(cylinder)
+
+    # interference function
+    interference = InterferenceFunction1DParaCrystal(distance, 3 * nanometer)
+    particle_decoration.addInterferenceFunction(interference)
+
+    # air layer with particles and substrate form multi layer
+    air_layer = Layer(mAir)
+    air_layer.setDecoration(particle_decoration)
+    substrate_layer = Layer(mSubstrate)
+    multi_layer = MultiLayer()
+    multi_layer.addLayer(air_layer)
+    multi_layer.addLayer(substrate_layer)
+
+    # build and run experiment
+    simulation = Simulation()
+    simulation.setDetectorParameters(100, -4.0 * degree, 4.0 * degree, 100, 0.0 * degree, 8.0 * degree)
+    simulation.setBeamParameters(1.0 * angstrom, 0.2 * degree, 0.0 * degree)
+    simulation.setSample(multi_layer)
+    simulation.runSimulation()
+    # intensity data
+    return GetOutputData(simulation)
+
+def SetParameters(i):
+    global radius
+    global height
+    global distance
+    radius = (1. + (3.0/Nframes)*i) * nanometer
+    height = (1. + (4.0/Nframes)*i) * nanometer
+    distance = (10. - (1.0/Nframes)*i) * nanometer
+
+#-------------------------------------------------------------
+# main()
+#-------------------------------------------------------------
+if __name__ == '__main__':
+    files = []
+    fig = pylab.figure(figsize=(5,5))
+    ax = fig.add_subplot(111)
+    for i in range(Nframes):
+        SetParameters(i)
+        result = RunSimulation() + 1 # for log scale
+        ax.cla()
+        im = ax.imshow(numpy.rot90(result, 1), vmax=1e3,
+                 norm=matplotlib.colors.LogNorm(),
+                 extent=[-4.0, 4.0, 0, 8.0])
+        pylab.xlabel(r'$\phi_f$', fontsize=20)
+        pylab.ylabel(r'$\alpha_f$', fontsize=20)
+        if i==0:
+            pylab.colorbar(im)
+        fname = '_tmp%03d.png'%i
+        print 'Saving frame', fname
+        fig.savefig(fname)
+        files.append(fname)
+    os.system("mencoder 'mf://_tmp*.png' -mf type=png:fps=10 -ovc lavc -lavcopts vcodec=wmv2 -oac copy -o animation.mpg")
+    print 'Removing temporary files'
+    os.system("rm _tmp*")
diff --git a/Examples/Demos/simul_demo_movie2.py b/Examples/Demos/simul_demo_movie2.py
new file mode 100644
index 0000000000000000000000000000000000000000..f1bff9b5936b71198086128d78a67e0e0069701b
--- /dev/null
+++ b/Examples/Demos/simul_demo_movie2.py
@@ -0,0 +1,86 @@
+'''
+Simulation demo: make movie
+'''
+
+import os, sys, numpy, pylab, matplotlib
+
+from libBornAgainCore import *
+
+radius = 1
+layer_thickness = 1
+Nframes = 100
+Ngrowframes = 30
+
+# ----------------------------------
+# describe sample and run simulation
+# ----------------------------------
+def RunSimulation():
+    # defining materials
+    mAir = MaterialManager.getHomogeneousMaterial("Air", 0.0, 0.0)
+    mSubstrate = MaterialManager.getHomogeneousMaterial("Substrate", 6e-6, 2e-8)
+    mParticle = MaterialManager.getHomogeneousMaterial("Particle", 6e-4, 2e-8)
+
+    # collection of particles
+    semisphere_ff = FormFactorSphere(radius, radius)
+    semisphere = Particle(mParticle, semisphere_ff)
+    particle_decoration = ParticleDecoration()
+    particle_decoration.addParticle(semisphere)
+
+    # interference function
+    interference = InterferenceFunction1DParaCrystal(6 * nanometer, 1 * nanometer)
+    particle_decoration.addInterferenceFunction(interference)
+
+    # air layer with particles and substrate form multi layer
+    air_layer = Layer(mAir)
+    air_layer.setDecoration(particle_decoration)
+    film_layer = Layer(mParticle, layer_thickness)
+    substrate_layer = Layer(mSubstrate)
+    multi_layer = MultiLayer()
+    multi_layer.addLayer(air_layer)
+    multi_layer.addLayer(film_layer)
+    multi_layer.addLayer(substrate_layer)
+
+    # build and run experiment
+    simulation = Simulation()
+    simulation.setDetectorParameters(100, -4.0 * degree, 4.0 * degree, 100, 0.0 * degree, 8.0 * degree)
+    simulation.setBeamParameters(1.0 * angstrom, 0.2 * degree, 0.0 * degree)
+    simulation.setSample(multi_layer)
+    simulation.runSimulation()
+    # intensity data
+    return GetOutputData(simulation)
+
+def SetParameters(i):
+    global radius
+    global layer_thickness
+    if i<Ngrowframes:
+        radius = (1. + (3.0/Ngrowframes)*i) * nanometer
+        layer_thickness = 0.01 * nanometer
+    else:
+        radius = 4. * nanometer
+        layer_thickness = (0.01 + (0.5/(Nframes-Ngrowframes))*(i-Ngrowframes)) * nanometer
+
+#-------------------------------------------------------------
+# main()
+#-------------------------------------------------------------
+if __name__ == '__main__':
+    files = []
+    fig = pylab.figure(figsize=(5,5))
+    ax = fig.add_subplot(111)
+    for i in range(Nframes):
+        SetParameters(i)
+        result = RunSimulation() + 1 # for log scale
+        ax.cla()
+        im = ax.imshow(numpy.rot90(result, 1), vmax=1e3,
+                 norm=matplotlib.colors.LogNorm(),
+                 extent=[-4.0, 4.0, 0, 8.0])
+        pylab.xlabel(r'$\phi_f$', fontsize=20)
+        pylab.ylabel(r'$\alpha_f$', fontsize=20)
+        if i==0:
+            pylab.colorbar(im)
+        fname = '_tmp%03d.png'%i
+        print 'Saving frame', fname
+        fig.savefig(fname)
+        files.append(fname)
+    os.system("mencoder 'mf://_tmp*.png' -mf type=png:fps=10 -ovc lavc -lavcopts vcodec=wmv2 -oac copy -o animation2.mpg")
+    print 'Removing temporary files'
+    os.system("rm _tmp*")
diff --git a/Examples/Demos/simul_demo_polarization.py b/Examples/Demos/simul_demo_polarization.py
new file mode 100644
index 0000000000000000000000000000000000000000..7acff41c1a8b56ff12b8f5a2117dcf9e2defa7e2
--- /dev/null
+++ b/Examples/Demos/simul_demo_polarization.py
@@ -0,0 +1,94 @@
+'''
+Simulation demo: Cylinder and/or sphere on substrate
+'''
+
+import numpy, pylab, matplotlib
+
+from libBornAgainCore import *
+
+# ----------------------------------
+# describe sample and run simulation
+# ----------------------------------
+def RunSimulation():
+    # defining materials
+    magnetic_field_layer = kvector_t(0.0, 6.0, 0.0)
+    magnetic_field_particle = kvector_t(1.7, 1.7, 1.7)
+    mAir = MaterialManager.getHomogeneousMaterial("Air", 0.0, 0.0)
+    mSubstrate = MaterialManager.getHomogeneousMaterial("Substrate", 6e-6, 2e-8)
+    mLayer = MaterialManager.getHomogeneousMagneticMaterial("Layer", 3e-6, 2e-8, magnetic_field_layer)
+    mParticle = MaterialManager.getHomogeneousMagneticMaterial("Particle", 6e-4, 2e-8, magnetic_field_particle)
+
+    # collection of particles
+    cylinder_ff = FormFactorCylinder(5 * nanometer, 2 * nanometer)
+    cylinder = Particle(mParticle, cylinder_ff)
+#    sphere_ff = FormFactorFullSphere(4 * nanometer)
+#    sphere = Particle(mParticle, sphere_ff)
+    particle_decoration = ParticleDecoration()
+    particle_decoration.addParticle(cylinder)
+#    particle_decoration.addParticle(sphere)
+#    interference = InterferenceFunction1DParaCrystal(20 * nanometer, 2 * nanometer)
+#    particle_decoration.addInterferenceFunction(interference)
+
+    # air layer with particles and substrate form multi layer
+    air_layer = Layer(mAir)
+    intermediate_layer = Layer(mLayer)
+    intermediate_layer.setDecoration(particle_decoration)
+    substrate_layer = Layer(mSubstrate)
+    multi_layer = MultiLayer()
+    multi_layer.addLayer(air_layer)
+    multi_layer.addLayer(intermediate_layer)
+    multi_layer.addLayer(substrate_layer)
+
+    # build and run experiment
+    simulation = Simulation()
+    simulation.setDetectorParameters(100, -4.0 * degree, 4.0 * degree, 100, 0.0 * degree, 8.0 * degree)
+    simulation.setBeamParameters(1.0 * angstrom, 0.2 * degree, 0.0 * degree)
+    simulation.setSample(multi_layer)
+    simulation.runSimulation()
+    # intensity data components
+
+    intensity_pp = GetPolarizedOutputDataComponent(simulation, 0, 0)
+    intensity_pm = GetPolarizedOutputDataComponent(simulation, 0, 1)
+    intensity_mp = GetPolarizedOutputDataComponent(simulation, 1, 0)
+    intensity_mm = GetPolarizedOutputDataComponent(simulation, 1, 1)
+    return intensity_pp, intensity_pm, intensity_mp, intensity_mm
+
+
+#-------------------------------------------------------------
+# main()
+#-------------------------------------------------------------
+if __name__ == '__main__':
+    intensity_pp, intensity_pm, intensity_mp, intensity_mm = RunSimulation()
+    pylab.subplot(2, 2, 1)
+    pylab.title('Intensity ++')
+    im = pylab.imshow(numpy.rot90(intensity_pp+1, 1),
+                 norm=matplotlib.colors.LogNorm(),
+                 extent=[-4.0, 4.0, 0, 8.0])
+    pylab.colorbar(im)
+    pylab.xlabel(r'$\phi_f$', fontsize=20)
+    pylab.ylabel(r'$\alpha_f$', fontsize=20)
+    pylab.subplot(2, 2, 2)
+    pylab.title('Intensity +-')
+    im = pylab.imshow(numpy.rot90(intensity_pm+1, 1),
+                 norm=matplotlib.colors.LogNorm(),
+                 extent=[-4.0, 4.0, 0, 8.0])
+    pylab.colorbar(im)
+    pylab.xlabel(r'$\phi_f$', fontsize=20)
+    pylab.ylabel(r'$\alpha_f$', fontsize=20)
+    pylab.subplot(2, 2, 3)
+    pylab.title('Intensity -+')
+    im = pylab.imshow(numpy.rot90(intensity_mp+1, 1),
+                 norm=matplotlib.colors.LogNorm(),
+                 extent=[-4.0, 4.0, 0, 8.0])
+    pylab.colorbar(im)
+    pylab.xlabel(r'$\phi_f$', fontsize=20)
+    pylab.ylabel(r'$\alpha_f$', fontsize=20)
+    pylab.subplot(2, 2, 4)
+    pylab.title('Intensity --')
+    im = pylab.imshow(numpy.rot90(intensity_mm+1, 1),
+                 norm=matplotlib.colors.LogNorm(),
+                 extent=[-4.0, 4.0, 0, 8.0])
+    pylab.colorbar(im)
+    pylab.xlabel(r'$\phi_f$', fontsize=20)
+    pylab.ylabel(r'$\alpha_f$', fontsize=20)
+    pylab.show()