diff --git a/App/src/TestIsGISAXS1.cpp b/App/src/TestIsGISAXS1.cpp
index 4b155fd9079dac02154e8e8c77cd8c6f7c91bf51..0294ffdb3ca217c224f8eeeb581798890e9269b4 100644
--- a/App/src/TestIsGISAXS1.cpp
+++ b/App/src/TestIsGISAXS1.cpp
@@ -63,7 +63,7 @@ void TestIsGISAXS1::initializeSample()
             0.0, 0.5);
     particle_decoration.addNanoParticle(new NanoParticle(n_particle, new FormFactorPrism3(5*Units::nanometer, 5*Units::nanometer)),
             0.0, 0.5);
-    particle_decoration.setInterferenceFunction(new InterferenceFunctionNone());
+    particle_decoration.addInterferenceFunction(new InterferenceFunctionNone());
     LayerDecorator air_layer_decorator(air_layer, particle_decoration);
 
     p_multi_layer->addLayer(air_layer_decorator);
diff --git a/App/src/TestIsGISAXS10.cpp b/App/src/TestIsGISAXS10.cpp
index 53ad2c786ee7f390b1f0f595aecfffa226d6339a..8110e35d520ee163ce728528a0183fbc48cd14a7 100644
--- a/App/src/TestIsGISAXS10.cpp
+++ b/App/src/TestIsGISAXS10.cpp
@@ -59,7 +59,7 @@ void TestIsGISAXS10::initializeSample()
             7*Units::nanometer, 1e7*Units::nanometer);
     NanoParticleDecoration particle_decoration(
                 new NanoParticle(n_particle, new FormFactorCylinder(5*Units::nanometer, 5*Units::nanometer)));
-    particle_decoration.setInterferenceFunction(p_interference_funtion);
+    particle_decoration.addInterferenceFunction(p_interference_funtion);
     LayerDecorator air_layer_decorator(air_layer, particle_decoration);
 
     p_multi_layer->addLayer(air_layer_decorator);
@@ -86,7 +86,7 @@ void TestIsGISAXS10::initializeSample2()
     NanoParticleDecoration particle_decoration(
                 new NanoParticle(n_particle, new FormFactorCylinder(5*Units::nanometer, 5*Units::nanometer)),
                 8*Units::nanometer);
-    particle_decoration.setInterferenceFunction(p_interference_funtion);
+    particle_decoration.addInterferenceFunction(p_interference_funtion);
     LayerDecorator substrate_layer_decorator(substrate_layer, particle_decoration);
 
     p_multi_layer->addLayer(air_layer);
diff --git a/App/src/TestIsGISAXS3.cpp b/App/src/TestIsGISAXS3.cpp
index eabdddda0196c6db9bcb044e10cfaf5ce7f99e8c..5219e3bcccf1bc78dd36d5a8b04e355d3cca80ae 100644
--- a/App/src/TestIsGISAXS3.cpp
+++ b/App/src/TestIsGISAXS3.cpp
@@ -60,7 +60,7 @@ void TestIsGISAXS3::initializeSample()
     substrate_layer.setMaterial(p_substrate_material);
     NanoParticleDecoration particle_decoration(
                 new NanoParticle(n_particle, new FormFactorCylinder(5*Units::nanometer, 5*Units::nanometer)));
-    particle_decoration.setInterferenceFunction(new InterferenceFunctionNone());
+    particle_decoration.addInterferenceFunction(new InterferenceFunctionNone());
     LayerDecorator air_layer_decorator(air_layer, particle_decoration);
 
     p_multi_layer->addLayer(air_layer_decorator);
diff --git a/Core/Algorithms/inc/DWBASimulation.h b/Core/Algorithms/inc/DWBASimulation.h
index f6d442b397b75068964bac1e0eb973f106d013ad..ad10c6db4695fbaaa0a75965d3c6aa3837da7b0c 100644
--- a/Core/Algorithms/inc/DWBASimulation.h
+++ b/Core/Algorithms/inc/DWBASimulation.h
@@ -29,7 +29,6 @@ public:
     OutputData<double> *getDWBAIntensity();
 protected:
     OutputData<double> m_dwba_intensity;
-    OutputData<complex_t> m_dwba_amplitude;
     kvector_t m_ki;
     double m_alpha_i;
 };
diff --git a/Core/Algorithms/inc/DecouplingApproximationStrategy.h b/Core/Algorithms/inc/DecouplingApproximationStrategy.h
new file mode 100644
index 0000000000000000000000000000000000000000..fc4a1dcba73d8a120eb6c5cdfe7f67ffe1212fc9
--- /dev/null
+++ b/Core/Algorithms/inc/DecouplingApproximationStrategy.h
@@ -0,0 +1,35 @@
+#ifndef DECOUPLINGAPPROXIMATIONSTRATEGY_H_
+#define DECOUPLINGAPPROXIMATIONSTRATEGY_H_
+// ********************************************************************
+// * The BornAgain project                                            *
+// * Simulation of neutron and x-ray scattering at grazing incidence  *
+// *                                                                  *
+// * LICENSE AND DISCLAIMER                                           *
+// * Lorem ipsum dolor sit amet, consectetur adipiscing elit.  Mauris *
+// * eget quam orci. Quisque  porta  varius  dui,  quis  posuere nibh *
+// * mollis quis. Mauris commodo rhoncus porttitor.                   *
+// ********************************************************************
+//! @file   DecouplingApproximationStrategy.h
+//! @brief  Definition of DecouplingApproximationStrategy class
+//! @author Scientific Computing Group at FRM II
+//! @date   Jun 29, 2012
+
+#include "IInterferenceFunctionStrategy.h"
+
+class DecouplingApproximationStrategy : public IInterferenceFunctionStrategy
+{
+public:
+    DecouplingApproximationStrategy();
+    virtual ~DecouplingApproximationStrategy() {}
+
+    virtual void init(const std::vector<IFormFactor *> &form_factors,
+            const std::vector<double> &fractions,
+            const std::vector<IInterferenceFunction *> &interference_functions);
+    virtual double evaluateForComplexkz(kvector_t k_i, kvector_t k_f,
+            complex_t k_iz, complex_t k_fz) const;
+private:
+    bool checkVectorSizes();
+};
+
+
+#endif /* DECOUPLINGAPPROXIMATIONSTRATEGY_H_ */
diff --git a/Core/Algorithms/inc/FormFactorSLDDecorator.h b/Core/Algorithms/inc/FormFactorSLDDecorator.h
new file mode 100644
index 0000000000000000000000000000000000000000..215e82483252b351ac38a10acb4e2aa972b7cfcf
--- /dev/null
+++ b/Core/Algorithms/inc/FormFactorSLDDecorator.h
@@ -0,0 +1,72 @@
+#ifndef FORMFACTORSLDDECORATOR_H_
+#define FORMFACTORSLDDECORATOR_H_
+// ********************************************************************
+// * The BornAgain project                                            *
+// * Simulation of neutron and x-ray scattering at grazing incidence  *
+// *                                                                  *
+// * LICENSE AND DISCLAIMER                                           *
+// * Lorem ipsum dolor sit amet, consectetur adipiscing elit.  Mauris *
+// * eget quam orci. Quisque  porta  varius  dui,  quis  posuere nibh *
+// * mollis quis. Mauris commodo rhoncus porttitor.                   *
+// ********************************************************************
+//! @file   FormFactorSLDDecorator.h
+//! @brief  Definition of FormFactorSLDDecorator class
+//! @author Scientific Computing Group at FRM II
+//! @date   Jun 29, 2012
+
+#include "IFormFactor.h"
+
+class FormFactorSLDDecorator : public IFormFactor
+{
+public:
+    FormFactorSLDDecorator(IFormFactor *p_form_factor, complex_t sld);
+    virtual FormFactorSLDDecorator *clone() const;
+    virtual ~FormFactorSLDDecorator();
+
+    virtual complex_t evaluate(kvector_t k_i, kvector_t k_f) const;
+    virtual complex_t evaluateForComplexkz(kvector_t k_i, kvector_t k_f, complex_t k_iz, complex_t k_fz) const;
+    virtual int getNumberOfStochasticParameters() const;
+private:
+    IFormFactor *mp_form_factor;
+    complex_t m_sld;
+};
+
+inline FormFactorSLDDecorator::FormFactorSLDDecorator(
+        IFormFactor* p_form_factor, complex_t sld)
+: mp_form_factor(p_form_factor)
+, m_sld(sld)
+{
+}
+
+inline FormFactorSLDDecorator* FormFactorSLDDecorator::clone() const
+{
+    return new FormFactorSLDDecorator(mp_form_factor->clone(), m_sld);
+}
+
+inline FormFactorSLDDecorator::~FormFactorSLDDecorator()
+{
+    delete mp_form_factor;
+}
+
+inline complex_t FormFactorSLDDecorator::evaluate(kvector_t k_i,
+        kvector_t k_f) const
+{
+    return m_sld*mp_form_factor->evaluate(k_i, k_f);
+}
+
+inline complex_t FormFactorSLDDecorator::evaluateForComplexkz(kvector_t k_i,
+        kvector_t k_f, complex_t k_iz, complex_t k_fz) const
+{
+    return m_sld*mp_form_factor->evaluateForComplexkz(k_i, k_f, k_iz, k_fz);
+}
+
+inline int FormFactorSLDDecorator::getNumberOfStochasticParameters() const
+{
+    return mp_form_factor->getNumberOfStochasticParameters();
+}
+
+
+
+
+
+#endif /* FORMFACTORSLDDECORATOR_H_ */
diff --git a/Core/Algorithms/inc/IInterferenceFunctionStrategy.h b/Core/Algorithms/inc/IInterferenceFunctionStrategy.h
new file mode 100644
index 0000000000000000000000000000000000000000..3d7963c53242be8fcb58f378e6aebb0b97a129f0
--- /dev/null
+++ b/Core/Algorithms/inc/IInterferenceFunctionStrategy.h
@@ -0,0 +1,72 @@
+#ifndef IINTERFERENCEFUNCTIONSTRATEGY_H_
+#define IINTERFERENCEFUNCTIONSTRATEGY_H_
+// ********************************************************************
+// * The BornAgain project                                            *
+// * Simulation of neutron and x-ray scattering at grazing incidence  *
+// *                                                                  *
+// * LICENSE AND DISCLAIMER                                           *
+// * Lorem ipsum dolor sit amet, consectetur adipiscing elit.  Mauris *
+// * eget quam orci. Quisque  porta  varius  dui,  quis  posuere nibh *
+// * mollis quis. Mauris commodo rhoncus porttitor.                   *
+// ********************************************************************
+//! @file   IInterferenceFunctionStrategy.h
+//! @brief  Definition of IInterferenceFunctionStrategy class
+//! @author Scientific Computing Group at FRM II
+//! @date   Jun 29, 2012
+
+#include "IFormFactor.h"
+#include "IInterferenceFunction.h"
+
+#include <vector>
+
+class IInterferenceFunctionStrategy
+{
+public:
+    virtual ~IInterferenceFunctionStrategy();
+    virtual void init(const std::vector<IFormFactor *> &form_factors,
+            const std::vector<double> &fractions,
+            const std::vector<IInterferenceFunction *> &interference_functions);
+    virtual double evaluateForComplexkz(kvector_t k_i, kvector_t k_f,
+            complex_t k_iz, complex_t k_fz) const=0;
+protected:
+    void deleteVectors();
+    std::vector<IFormFactor*> m_form_factors; //!< Includes Scattering Length Density
+    std::vector<double> m_fractions;
+    std::vector<IInterferenceFunction*> m_interference_functions;
+};
+
+inline IInterferenceFunctionStrategy::~IInterferenceFunctionStrategy()
+{
+    deleteVectors();
+}
+
+inline void IInterferenceFunctionStrategy::init(
+        const std::vector<IFormFactor*>& form_factors,
+        const std::vector<double>& fractions,
+        const std::vector<IInterferenceFunction*>& interference_functions)
+{
+    deleteVectors();
+    m_fractions = fractions;
+    for (size_t i=0; i<form_factors.size(); ++i) {
+        m_form_factors.push_back(form_factors[i]->clone());
+    }
+    for (size_t i=0; i<interference_functions.size(); ++i) {
+        m_interference_functions.push_back(interference_functions[i]->clone());
+    }
+}
+
+inline void IInterferenceFunctionStrategy::deleteVectors()
+{
+    for (size_t i=0; i<m_form_factors.size(); ++i) {
+        delete m_form_factors[i];
+    }
+    for (size_t i=0; i<m_interference_functions.size(); ++i) {
+        delete m_interference_functions[i];
+    }
+    m_form_factors.clear();
+    m_interference_functions.clear();
+}
+
+
+
+#endif /* IINTERFERENCEFUNCTIONSTRATEGY_H_ */
diff --git a/Core/Algorithms/inc/LocalMonodisperseApproximationStrategy.h b/Core/Algorithms/inc/LocalMonodisperseApproximationStrategy.h
new file mode 100644
index 0000000000000000000000000000000000000000..c84c786c0240b36dfb93ff13013b67487a333c38
--- /dev/null
+++ b/Core/Algorithms/inc/LocalMonodisperseApproximationStrategy.h
@@ -0,0 +1,36 @@
+#ifndef LOCALMONODISPERSEAPPROXIMATIONSTRATEGY_H_
+#define LOCALMONODISPERSEAPPROXIMATIONSTRATEGY_H_
+// ********************************************************************
+// * The BornAgain project                                            *
+// * Simulation of neutron and x-ray scattering at grazing incidence  *
+// *                                                                  *
+// * LICENSE AND DISCLAIMER                                           *
+// * Lorem ipsum dolor sit amet, consectetur adipiscing elit.  Mauris *
+// * eget quam orci. Quisque  porta  varius  dui,  quis  posuere nibh *
+// * mollis quis. Mauris commodo rhoncus porttitor.                   *
+// ********************************************************************
+//! @file   LocalMonodisperseApproximationStrategy.h
+//! @brief  Definition of LocalMonodisperseApproximationStrategy class
+//! @author Scientific Computing Group at FRM II
+//! @date   Jun 29, 2012
+
+#include "IInterferenceFunctionStrategy.h"
+
+class LocalMonodisperseApproximationStrategy : public IInterferenceFunctionStrategy
+{
+public:
+    LocalMonodisperseApproximationStrategy();
+    virtual ~LocalMonodisperseApproximationStrategy() {}
+
+    virtual void init(const std::vector<IFormFactor *> &form_factors,
+            const std::vector<double> &fractions,
+            const std::vector<IInterferenceFunction *> &interference_functions);
+    virtual double evaluateForComplexkz(kvector_t k_i, kvector_t k_f,
+            complex_t k_iz, complex_t k_fz) const;
+private:
+    bool checkVectorSizes();
+};
+
+
+
+#endif /* LOCALMONODISPERSEAPPROXIMATIONSTRATEGY_H_ */
diff --git a/Core/Algorithms/src/DWBASimulation.cpp b/Core/Algorithms/src/DWBASimulation.cpp
index 8c19c58f6cb2f17e1aab0df49b4ec543e63480bb..7caa0f47995cccc51c8ee97f877cb242fc041420 100644
--- a/Core/Algorithms/src/DWBASimulation.cpp
+++ b/Core/Algorithms/src/DWBASimulation.cpp
@@ -11,12 +11,10 @@ DWBASimulation::~DWBASimulation()
 void DWBASimulation::init(const Experiment& experiment)
 {
     m_dwba_intensity.clear();
-    m_dwba_amplitude.clear();
     Detector detector = experiment.getDetector();
     size_t detector_dimension = detector.getDimension();
     for (size_t dim=0; dim<detector_dimension; ++dim) {
         m_dwba_intensity.addAxis(new NamedVector<double>(detector.getAxis(dim)));
-        m_dwba_amplitude.addAxis(new NamedVector<double>(detector.getAxis(dim)));
     }
     Beam beam = experiment.getBeam();
     m_ki = beam.getCentralK();
diff --git a/Core/Algorithms/src/DecouplingApproximationStrategy.cpp b/Core/Algorithms/src/DecouplingApproximationStrategy.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a512da196f90e1cb558fc11c5560db73e228cad2
--- /dev/null
+++ b/Core/Algorithms/src/DecouplingApproximationStrategy.cpp
@@ -0,0 +1,42 @@
+#include "DecouplingApproximationStrategy.h"
+#include "Exceptions.h"
+
+DecouplingApproximationStrategy::DecouplingApproximationStrategy()
+{
+}
+
+void DecouplingApproximationStrategy::init(
+        const std::vector<IFormFactor*>& form_factors,
+        const std::vector<double>& fractions,
+        const std::vector<IInterferenceFunction*>& interference_functions)
+{
+    IInterferenceFunctionStrategy::init(form_factors, fractions, interference_functions);
+    if (!checkVectorSizes()) {
+        throw ClassInitializationException("Wrong number of formfactors or interference functions for Decoupling Approximation.");
+    }
+}
+
+double DecouplingApproximationStrategy::evaluateForComplexkz(kvector_t k_i,
+        kvector_t k_f, complex_t k_iz, complex_t k_fz) const
+{
+    double intensity = 0.0;
+    complex_t amplitude = complex_t(0.0, 0.0);
+    for (size_t i=0; i<m_form_factors.size(); ++i) {
+        complex_t ff = m_form_factors[i]->evaluateForComplexkz(k_i, k_f, k_iz, k_fz);
+        double fraction = m_fractions[i];
+        amplitude += fraction*ff;
+        intensity += fraction*std::norm(ff);
+    }
+    double amplitude_norm = std::norm(amplitude);
+    double itf_function = m_interference_functions[0]->evaluate(k_i-k_f);
+    return intensity + amplitude_norm*(itf_function-1.0);
+}
+
+bool DecouplingApproximationStrategy::checkVectorSizes()
+{
+    size_t n_ffs = m_form_factors.size();
+    size_t n_frs = m_fractions.size();
+    size_t n_ifs = m_interference_functions.size();
+    return (n_ffs==n_frs && n_ifs==1);
+}
+
diff --git a/Core/Algorithms/src/LayerDecoratorDWBASimulation.cpp b/Core/Algorithms/src/LayerDecoratorDWBASimulation.cpp
index a0337528f75b43e84e049c9e5660a1612a1a7ba0..5c3cdc340a486d4b771284a37ed7dbf58effa7de 100644
--- a/Core/Algorithms/src/LayerDecoratorDWBASimulation.cpp
+++ b/Core/Algorithms/src/LayerDecoratorDWBASimulation.cpp
@@ -1,6 +1,7 @@
 #include "LayerDecoratorDWBASimulation.h"
 #include "LayerDecorator.h"
 #include "DWBAFormFactorConstZ.h"
+#include "FormFactorSLDDecorator.h"
 
 LayerDecoratorDWBASimulation::LayerDecoratorDWBASimulation(
         const LayerDecorator *p_layer_decorator)
@@ -15,51 +16,37 @@ LayerDecoratorDWBASimulation::~LayerDecoratorDWBASimulation()
 
 void LayerDecoratorDWBASimulation::run()
 {
-    // Currently, this is the decoupling approximation
     m_dwba_intensity.resetIndex();
-    m_dwba_intensity.setAllTo(0.0);
-    m_dwba_amplitude.resetIndex();
-    m_dwba_amplitude.setAllTo(complex_t(0.0, 0.0));
     double lambda = 2.0*M_PI/m_ki.mag();
     complex_t k_iz = -mp_kz_function->evaluate(-m_alpha_i);
     const NanoParticleDecoration *p_decoration = mp_layer_decorator->getDecoration();
     complex_t n_layer = mp_layer_decorator->getRefractiveIndex();
     size_t number_of_particles = p_decoration->getNumberOfParticles();
+    std::vector<IFormFactor *> form_factors;
     for (size_t particle_index=0; particle_index<number_of_particles; ++particle_index) {
         NanoParticle *p_particle = p_decoration->getNanoParticle(particle_index);
         double depth = p_decoration->getDepthOfNanoParticle(particle_index);
         double m_abundance_fraction = p_decoration->getAbundanceFractionOfNanoParticle(particle_index);
         complex_t n_decoration = p_particle->getRefractiveIndex();
         complex_t scattering_length_density = (n_layer*n_layer - n_decoration*n_decoration)*M_PI/lambda/lambda;
-        double normalizing_factor = m_abundance_fraction*std::norm(scattering_length_density);
         DWBAFormFactorConstZ dwba_z(p_particle->getFormFactor()->clone(), depth);
         dwba_z.setReflectionFunction(mp_R_function);
         dwba_z.setTransmissionFunction(mp_T_function);
-        while (m_dwba_intensity.hasNext() && m_dwba_amplitude.hasNext())
-        {
-            double phi_f = m_dwba_intensity.getCurrentValueOfAxis<double>("phi_f");
-            double alpha_f = m_dwba_intensity.getCurrentValueOfAxis<double>("alpha_f");
-            kvector_t k_f;
-            complex_t k_fz = mp_kz_function->evaluate(alpha_f);
-            k_f.setLambdaAlphaPhi(lambda, alpha_f, phi_f);
-            complex_t ff = scattering_length_density*dwba_z.evaluateForComplexkz(m_ki, k_f, k_iz, k_fz);
-            m_dwba_intensity.next() += m_abundance_fraction*std::norm(ff);
-            m_dwba_amplitude.next() += m_abundance_fraction*ff;
-        }
-        m_dwba_intensity.resetIndex();
-        m_dwba_amplitude.resetIndex();
+        FormFactorSLDDecorator *p_ff = new FormFactorSLDDecorator(dwba_z.clone(), scattering_length_density);
+        form_factors.push_back(p_ff);
     }
-    if (number_of_particles!=0) {
-        const IInterferenceFunction *p_interference_function = p_decoration->getInterferenceFunction();
-        while (m_dwba_intensity.hasNext())
-        {
-            double phi_f = m_dwba_intensity.getCurrentValueOfAxis<double>("phi_f");
-            double alpha_f = m_dwba_intensity.getCurrentValueOfAxis<double>("alpha_f");
-            kvector_t k_f;
-            k_f.setLambdaAlphaPhi(lambda, alpha_f, phi_f);
-            complex_t mean_amplitude = m_dwba_amplitude.next();
-            double amplitude_term = std::norm(mean_amplitude)*(p_interference_function->evaluate(m_ki-k_f) - 1.0);
-            m_dwba_intensity.next() += amplitude_term;
-        }
+    IInterferenceFunctionStrategy *p_strategy = p_decoration->createStrategy(form_factors);
+    for (size_t i=0; i<form_factors.size(); ++i) {
+        delete form_factors[i];
+    }
+
+    while (m_dwba_intensity.hasNext())
+    {
+        double phi_f = m_dwba_intensity.getCurrentValueOfAxis<double>("phi_f");
+        double alpha_f = m_dwba_intensity.getCurrentValueOfAxis<double>("alpha_f");
+        kvector_t k_f;
+        complex_t k_fz = mp_kz_function->evaluate(alpha_f);
+        k_f.setLambdaAlphaPhi(lambda, alpha_f, phi_f);
+        m_dwba_intensity.next() = p_strategy->evaluateForComplexkz(m_ki, k_f, k_iz, k_fz);
     }
 }
diff --git a/Core/Algorithms/src/LocalMonodisperseApproximationStrategy.cpp b/Core/Algorithms/src/LocalMonodisperseApproximationStrategy.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8e4bf05d24b4d14837bf4ad03a4dd203fa5090ab
--- /dev/null
+++ b/Core/Algorithms/src/LocalMonodisperseApproximationStrategy.cpp
@@ -0,0 +1,38 @@
+#include "LocalMonodisperseApproximationStrategy.h"
+#include "Exceptions.h"
+
+LocalMonodisperseApproximationStrategy::LocalMonodisperseApproximationStrategy()
+{
+}
+
+void LocalMonodisperseApproximationStrategy::init(
+        const std::vector<IFormFactor*>& form_factors,
+        const std::vector<double>& fractions,
+        const std::vector<IInterferenceFunction*>& interference_functions)
+{
+    IInterferenceFunctionStrategy::init(form_factors, fractions, interference_functions);
+    if (!checkVectorSizes()) {
+        throw ClassInitializationException("Wrong number of formfactors or interference functions for Decoupling Approximation.");
+    }
+}
+
+double LocalMonodisperseApproximationStrategy::evaluateForComplexkz(
+        kvector_t k_i, kvector_t k_f, complex_t k_iz, complex_t k_fz) const
+{
+    double intensity = 0.0;
+    for (size_t i=0; i<m_form_factors.size(); ++i) {
+        complex_t ff = m_form_factors[i]->evaluateForComplexkz(k_i, k_f, k_iz, k_fz);
+        double itf_function = m_interference_functions[i]->evaluate(k_i-k_f);
+        double fraction = m_fractions[i];
+        intensity += fraction*itf_function*std::norm(ff);
+    }
+    return intensity;
+}
+
+bool LocalMonodisperseApproximationStrategy::checkVectorSizes()
+{
+    size_t n_ffs = m_form_factors.size();
+    size_t n_frs = m_fractions.size();
+    size_t n_ifs = m_interference_functions.size();
+    return (n_ffs==n_frs && n_ifs==n_ffs);
+}
diff --git a/Core/Core.pro b/Core/Core.pro
index ff7924206000ac40ed07ee8dabb0f7385e8bd242..4fea7c53cebfc075f8310e53a73a570c7612c971 100644
--- a/Core/Core.pro
+++ b/Core/Core.pro
@@ -12,6 +12,7 @@ QMAKE_EXTENSION_SHLIB = so
 
 SOURCES += \
     Algorithms/src/Beam.cpp \
+    Algorithms/src/DecouplingApproximationStrategy.cpp \
     Algorithms/src/Detector.cpp \
     Algorithms/src/DWBADiffuseReflection.cpp \
     Algorithms/src/DWBAFormFactor.cpp \
@@ -21,6 +22,7 @@ SOURCES += \
     Algorithms/src/GISASExperiment.cpp \
     Algorithms/src/LayerDecoratorDWBASimulation.cpp \
     Algorithms/src/LayerDWBASimulation.cpp \
+    Algorithms/src/LocalMonodisperseApproximationStrategy.cpp \
     Algorithms/src/MultiLayerDWBASimulation.cpp \
     Algorithms/src/OpticalFresnel.cpp \
     Algorithms/src/PythonAlgorithmsInterface.cpp \
@@ -63,6 +65,7 @@ SOURCES += \
 
 HEADERS += \
     Algorithms/inc/Beam.h \
+    Algorithms/inc/DecouplingApproximationStrategy.h \
     Algorithms/inc/Detector.h \
     Algorithms/inc/DWBADiffuseReflection.h \
     Algorithms/inc/DWBAFormFactor.h \
@@ -71,8 +74,10 @@ HEADERS += \
     Algorithms/inc/Experiment.h \
     Algorithms/inc/GISASExperiment.h \
     Algorithms/inc/ISimulation.h \
+    Algorithms/inc/IInterferenceFunctionStrategy.h \
     Algorithms/inc/LayerDecoratorDWBASimulation.h \
     Algorithms/inc/LayerDWBASimulation.h \
+    Algorithms/inc/LocalMonodisperseApproximationStrategy.h \
     Algorithms/inc/MultiLayerDWBASimulation.h \
     Algorithms/inc/OpticalFresnel.h \
     Algorithms/inc/PythonAlgorithmsInterface.h \
diff --git a/Core/Samples/inc/NanoParticleDecoration.h b/Core/Samples/inc/NanoParticleDecoration.h
index 3ee3ecc83be49417bb5715cb3b9eb6fa6433018b..b099004e4f195e2b0356732b36eabbd9a0400125 100644
--- a/Core/Samples/inc/NanoParticleDecoration.h
+++ b/Core/Samples/inc/NanoParticleDecoration.h
@@ -17,6 +17,7 @@
 #include "IDecoration.h"
 #include "NanoParticle.h"
 #include "IInterferenceFunction.h"
+#include "IInterferenceFunctionStrategy.h"
 
 //- -------------------------------------------------------------------
 //! @class NanoParticleDecoration
@@ -47,14 +48,14 @@ public:
     /// Get abundance fraction of particle with index
     double getAbundanceFractionOfNanoParticle(size_t index) const;
 
-    /// Set interference function
-    void setInterferenceFunction(IInterferenceFunction* p_interference_function);
+    /// Add interference function
+    void addInterferenceFunction(IInterferenceFunction* p_interference_function);
 
-    /// Get interference function
-    const IInterferenceFunction* getInterferenceFunction() const
-    {
-        return mp_interference_function;
-    }
+    /// Get interference function with index
+    const IInterferenceFunction* getInterferenceFunction(size_t index) const;
+
+    /// Create interference function strategy
+    IInterferenceFunctionStrategy *createStrategy(const std::vector<IFormFactor *> &form_factors) const;
 
 private:
     struct ParticleInfoStruct
@@ -70,7 +71,7 @@ private:
     };
     std::vector<ParticleInfoStruct> m_particles;
     ///< Vector of the types of nano particles
-    IInterferenceFunction* mp_interference_function;
+    std::vector<IInterferenceFunction*> m_interference_functions;
     ///< Currently only a scalar interference function (instead of matrix)
     double m_total_abundance;
     ///< To guarantee that fractions sum up to 1
diff --git a/Core/Samples/src/NanoParticleDecoration.cpp b/Core/Samples/src/NanoParticleDecoration.cpp
index 666b9f928fbbfd1addde11562c8321e8b72565a7..41cf070dcd66f0fa39cf0d8b81f2590b1aab0f86 100644
--- a/Core/Samples/src/NanoParticleDecoration.cpp
+++ b/Core/Samples/src/NanoParticleDecoration.cpp
@@ -1,33 +1,37 @@
 #include "NanoParticleDecoration.h"
 #include "InterferenceFunctionNone.h"
+#include "DecouplingApproximationStrategy.h"
+#include "LocalMonodisperseApproximationStrategy.h"
 
 /* ************************************************************************* */
 NanoParticleDecoration::NanoParticleDecoration()
-: mp_interference_function(0)
-, m_total_abundance(0.0)
+: m_total_abundance(0.0)
 {
 }
 
 NanoParticleDecoration::NanoParticleDecoration(NanoParticle* p_particle, double depth, double abundance)
-: mp_interference_function(new InterferenceFunctionNone)
-, m_total_abundance(0.0)
+: m_total_abundance(0.0)
 {
 	addNanoParticle(p_particle, depth, abundance);
 }
 
 NanoParticleDecoration::~NanoParticleDecoration()
 {
-    delete mp_interference_function;
+    for (size_t i=0; i<m_interference_functions.size(); ++i) {
+        delete m_interference_functions[i];
+    }
 }
 
 NanoParticleDecoration* NanoParticleDecoration::clone() const
 {
     NanoParticleDecoration *p_new = new NanoParticleDecoration();
-    p_new->setInterferenceFunction(mp_interference_function->clone());
     for (size_t index=0; index<m_particles.size(); ++index) {
         p_new->addNanoParticle(m_particles[index].mp_particle->clone(), m_particles[index].m_depth,
                 m_particles[index].m_abundance);
     }
+    for (size_t i=0; i<m_interference_functions.size(); ++i) {
+        p_new->addInterferenceFunction(m_interference_functions[i]->clone());
+    }
     return p_new;
 }
 
@@ -63,13 +67,42 @@ double NanoParticleDecoration::getAbundanceFractionOfNanoParticle(
     throw OutOfBoundsException("Not so many particles in this decoration.");
 }
 
-void NanoParticleDecoration::setInterferenceFunction(
+void NanoParticleDecoration::addInterferenceFunction(
 		IInterferenceFunction* p_interference_function)
 {
-	if (mp_interference_function!=p_interference_function) {
-		delete mp_interference_function;
-		mp_interference_function = p_interference_function;
-	}
+    m_interference_functions.push_back(p_interference_function);
+}
+
+const IInterferenceFunction* NanoParticleDecoration::getInterferenceFunction(
+        size_t index) const
+{
+    if (index<m_interference_functions.size()) {
+        return m_interference_functions[index]->clone();
+    }
+    throw OutOfBoundsException("Not so many interference functions in this decoration.");
+}
+
+IInterferenceFunctionStrategy* NanoParticleDecoration::createStrategy(
+        const std::vector<IFormFactor*>& form_factors) const
+{
+    std::vector<double> fractions;
+    for (size_t i=0; i<m_particles.size(); ++i) {
+        fractions.push_back(getAbundanceFractionOfNanoParticle(i));
+    }
+    IInterferenceFunctionStrategy *p_strategy;
+    size_t n_particles = m_particles.size();
+    size_t n_ifs = m_interference_functions.size();
+    if (n_ifs==1) {
+        p_strategy = new DecouplingApproximationStrategy();
+    }
+    else if (n_ifs==n_particles) {
+        p_strategy = new LocalMonodisperseApproximationStrategy();
+    }
+    else {
+        throw ClassInitializationException("Could not create interference function strategy with given parameters.");
+    }
+    p_strategy->init(form_factors, fractions, m_interference_functions);
+    return p_strategy;
 }
 
 NanoParticleDecoration::ParticleInfoStruct::ParticleInfoStruct(
@@ -108,5 +141,4 @@ const NanoParticleDecoration::ParticleInfoStruct& NanoParticleDecoration::Partic
 
 
 
-
 /* ************************************************************************* */