diff --git a/Core/Computation/MainComputation.cpp b/Core/Computation/MainComputation.cpp
index d3f6085c3247b034798b98e9519141ebf67a79a9..c5c1fb9c5ab844180e73b3ebf43cfa8136a4de2e 100644
--- a/Core/Computation/MainComputation.cpp
+++ b/Core/Computation/MainComputation.cpp
@@ -17,7 +17,6 @@
 #include "ParticleLayoutComputation.h"
 #include "Layer.h"
 #include "IFresnelMap.h"
-#include "Instrument.h"
 #include "MatrixFresnelMap.h"
 #include "MultiLayer.h"
 #include "RoughMultiLayerComputation.h"
@@ -27,21 +26,13 @@
 #include "SimulationElement.h"
 #include "MaterialFactoryFuncs.h"
 
-namespace
-{
-Material CalculateAverageMaterial(const Material& layer_mat, double wavelength,
-                                  const std::vector<HomogeneousRegion>& regions);
-}
-
 MainComputation::MainComputation(
     const MultiLayer& multilayer,
-    const Instrument& instrument,
     const SimulationOptions& options,
     ProgressHandler& progress,
     const std::vector<SimulationElement>::iterator& begin_it,
     const std::vector<SimulationElement>::iterator& end_it)
     : mP_multi_layer(multilayer.cloneSliced(options.useAvgMaterials()))
-    , m_instrument(instrument)
     , m_sim_options(options)
     , m_progress(&progress)
     , m_begin_it(begin_it)
@@ -129,8 +120,7 @@ std::unique_ptr<MultiLayer> MainComputation::getAveragedMultilayer() const
         if (!checkRegions(entry.second))
             throw std::runtime_error("MainComputation::getAveragedMultilayer: "
                                      "total volumetric fraction of particles exceeds 1!");
-        double beam_wavelength = m_instrument.getBeam().getWavelength();
-        auto new_mat = CalculateAverageMaterial(layer_mat, beam_wavelength, entry.second);
+        auto new_mat = createAveragedMaterial(layer_mat, entry.second);
         P_result->setLayerMaterial(i_layer, new_mat);
     }
     return P_result;
@@ -158,25 +148,3 @@ bool MainComputation::checkRegions(const std::vector<HomogeneousRegion>& regions
         total_fraction += region.m_volume;
     return (total_fraction >= 0 && total_fraction <= 1);
 }
-
-namespace
-{
-// TODO: make this procedure correct for all types of materials
-Material CalculateAverageMaterial(const Material& layer_mat, double wavelength,
-                                  const std::vector<HomogeneousRegion>& regions)
-{
-    kvector_t magnetization_layer = layer_mat.magnetization();
-    complex_t refr_index2_layer = layer_mat.refractiveIndex2(wavelength);
-    kvector_t magnetization_avg = magnetization_layer;
-    complex_t refr_index2_avg = refr_index2_layer;
-    for (auto& region : regions)
-    {
-        kvector_t magnetization_region = region.m_material.magnetization();
-        complex_t refr_index2_region = region.m_material.refractiveIndex2(wavelength);
-        magnetization_avg += region.m_volume*(magnetization_region - magnetization_layer);
-        refr_index2_avg += region.m_volume*(refr_index2_region - refr_index2_layer);
-    }
-    complex_t refr_index_avg = std::sqrt(refr_index2_avg);
-    return HomogeneousMaterial(layer_mat.getName() + "_avg", refr_index_avg, magnetization_avg);
-}
-}
diff --git a/Core/Computation/MainComputation.h b/Core/Computation/MainComputation.h
index 76b5a551cbaa8558c76f85b9350eb1d087a6ced2..373b3e2dc4c113426e9a55444820b4b55559eed8 100644
--- a/Core/Computation/MainComputation.h
+++ b/Core/Computation/MainComputation.h
@@ -25,7 +25,6 @@
 
 class IFresnelMap;
 class MultiLayer;
-class Instrument;
 struct HomogeneousRegion;
 class IComputationTerm;
 class ProgressHandler;
@@ -43,7 +42,6 @@ class MainComputation final : public INoncopyable
 public:
     MainComputation(
         const MultiLayer& multilayer,
-        const Instrument& instrument,
         const SimulationOptions& options,
         ProgressHandler& progress,
         const std::vector<SimulationElement>::iterator& begin_it,
@@ -68,7 +66,6 @@ private:
     bool checkRegions(const std::vector<HomogeneousRegion>& regions) const;
 
     std::unique_ptr<MultiLayer> mP_multi_layer;
-    const Instrument& m_instrument;
     SimulationOptions m_sim_options;
     ProgressHandler* m_progress;
     //! these iterators define the span of detector bins this simulation will work on
diff --git a/Core/Export/ExportToPython.cpp b/Core/Export/ExportToPython.cpp
index 401c673df1df4eea89fa419b5b6287218319d310..e9bc17c5744025e3493f1d65c013bb4bf29c683d 100644
--- a/Core/Export/ExportToPython.cpp
+++ b/Core/Export/ExportToPython.cpp
@@ -45,6 +45,7 @@
 #include "ParameterUtils.h"
 #include <iomanip>
 #include <set>
+#include <map>
 #include <functional>
 
 class IFormFactor;
@@ -183,6 +184,10 @@ std::string ExportToPython::defineGetSample() const
         + "\n";
 }
 
+const std::map<MATERIAL_TYPES, std::string> factory_names{
+    {MATERIAL_TYPES::RefractiveMaterial, "HomogeneousMaterial"},
+    {MATERIAL_TYPES::MaterialBySLD, "MaterialBySLD"}};
+
 std::string ExportToPython::defineMaterials() const
 {
     const auto themap = m_label->materialMap();
@@ -197,22 +202,24 @@ std::string ExportToPython::defineMaterials() const
             continue;
         visitedMaterials.insert(it->second);
         const Material* p_material = it->first;
-        complex_t material_data = p_material->materialData();
-        double real = std::real(material_data);
-        double imag = std::imag(material_data);
+        const auto factory_name = factory_names.find(p_material->typeID());
+        if (factory_name == factory_names.cend())
+            throw std::runtime_error(
+                "Error in ExportToPython::defineMaterials(): unknown material type");
+        const complex_t& material_data = p_material->materialData();
         if (p_material->isScalarMaterial()) {
-            result << indent() << m_label->labelMaterial(p_material)
-                   << " = ba.HomogeneousMaterial(\"" << p_material->getName()
-                   << "\", " << printDouble(real) << ", "
-                   << printDouble(imag) << ")\n";
+            result << indent() << m_label->labelMaterial(p_material) << " = ba."
+                   << factory_name->second << "(\"" << p_material->getName() << "\", "
+                   << printDouble(material_data.real()) << ", " << printDouble(material_data.imag())
+                   << ")\n";
         } else {
             kvector_t magnetic_field = p_material->magnetization();
             result << indent() << "magnetic_field = kvector_t(" << magnetic_field.x() << ", "
                    << magnetic_field.y() << ", " << magnetic_field.z() << ")\n";
-            result << indent() << m_label->labelMaterial(p_material)
-                   << " = ba.HomogeneousMaterial(\"" << p_material->getName();
-            result << "\", " << printDouble(real) << ", "
-                   << printDouble(imag) << ", "
+            result << indent() << m_label->labelMaterial(p_material) << " = ba."
+                   << factory_name->second << "(\"" << p_material->getName();
+            result << "\", " << printDouble(material_data.real()) << ", "
+                   << printDouble(material_data.imag()) << ", "
                    << "magnetic_field)\n";
         }
     }
diff --git a/Core/Material/BaseMaterialImpl.h b/Core/Material/BaseMaterialImpl.h
index 959165cbb3ef12b00b5d1b1cf87f565ffa8d8097..a583131f46378db313cc4aaf7c1af2061597b4c0 100644
--- a/Core/Material/BaseMaterialImpl.h
+++ b/Core/Material/BaseMaterialImpl.h
@@ -13,8 +13,8 @@
 //
 // ************************************************************************** //
 
-#ifndef IMATERIALIMPL_H_
-#define IMATERIALIMPL_H_
+#ifndef BASEMATERIALIMPL_H_
+#define BASEMATERIALIMPL_H_
 
 #include "INamed.h"
 #include "Vectors3D.h"
@@ -25,9 +25,8 @@ class Transform3D;
 class WavevectorInfo;
 
 enum class MATERIAL_TYPES {
-    VacuumMaterial = 0,
-    RefractiveCoefMaterial,
-    WavelengthIndependentMaterial
+    RefractiveMaterial = 0,
+    MaterialBySLD
 };
 
 //! @ingroup materials
@@ -79,4 +78,4 @@ public:
     virtual void print(std::ostream &ostr) const = 0;
 };
 
-#endif /* IMATERIALIMPL_H_ */
+#endif /* BASEMATERIALIMPL_H_ */
diff --git a/Core/Material/Material.cpp b/Core/Material/Material.cpp
index 848f2a08104e1663713cfa6bbf05fd82d34cc5b7..9d36b91e6352f7fed74365ee2150d27c9a135f51 100644
--- a/Core/Material/Material.cpp
+++ b/Core/Material/Material.cpp
@@ -69,6 +69,11 @@ complex_t Material::materialData() const
     return m_material_impl->materialData();
 }
 
+bool Material::isDefaultMaterial() const
+{
+    return materialData() == complex_t() && isScalarMaterial();
+}
+
 complex_t Material::scalarSubtrSLD(const WavevectorInfo& wavevectors) const
 {
     return m_material_impl->scalarSubtrSLD(wavevectors);
diff --git a/Core/Material/Material.h b/Core/Material/Material.h
index 9c0a021d3ccaef04522654e321769c469dbf183b..349854f0d70cb4d65ea9d68eb1a619ea238fb14f 100644
--- a/Core/Material/Material.h
+++ b/Core/Material/Material.h
@@ -81,6 +81,10 @@ public:
     //! Returns true if material underlying data is nullptr
     bool isEmpty() const {return !m_material_impl;}
 
+    //! Returns true if material has refractive index of (1.0, 0.0)
+    //! and zero magnetization.
+    bool isDefaultMaterial() const;
+
     //! Returns (\f$ \pi/\lambda^2 \f$ - sld), sld (in \f$nm^{-2}\f$) being the scattering length density
     complex_t scalarSubtrSLD(const WavevectorInfo& wavevectors) const;
 
diff --git a/Core/Material/WavelengthIndependentMaterial.cpp b/Core/Material/MaterialBySLDImpl.cpp
similarity index 56%
rename from Core/Material/WavelengthIndependentMaterial.cpp
rename to Core/Material/MaterialBySLDImpl.cpp
index bd14a58b0b00ff94938105d1c4fe6793cae5f8f5..61718346d057dd64fc6f5b40d6f279fae7461ee0 100644
--- a/Core/Material/WavelengthIndependentMaterial.cpp
+++ b/Core/Material/MaterialBySLDImpl.cpp
@@ -1,4 +1,4 @@
-#include "WavelengthIndependentMaterial.h"
+#include "MaterialBySLDImpl.h"
 #include "WavevectorInfo.h"
 
 namespace
@@ -16,44 +16,41 @@ inline double getWlPrefactor(double wavelength)
 }
 }
 
-WavelengthIndependentMaterial::WavelengthIndependentMaterial(const std::string& name, double sld,
+MaterialBySLDImpl::MaterialBySLDImpl(const std::string& name, double sld,
                                                              double abs_term,
                                                              kvector_t magnetization)
     : MagneticMaterialImpl(name, magnetization), m_sld(sld), m_abs_term(abs_term)
 {}
 
-WavelengthIndependentMaterial::~WavelengthIndependentMaterial()
-{}
-
-WavelengthIndependentMaterial* WavelengthIndependentMaterial::clone() const
+MaterialBySLDImpl* MaterialBySLDImpl::clone() const
 {
-    return new WavelengthIndependentMaterial(*this);
+    return new MaterialBySLDImpl(*this);
 }
 
-complex_t WavelengthIndependentMaterial::refractiveIndex(double wavelength) const
+complex_t MaterialBySLDImpl::refractiveIndex(double wavelength) const
 {
     return std::sqrt(refractiveIndex2(wavelength));
 }
 
-complex_t WavelengthIndependentMaterial::refractiveIndex2(double wavelength) const
+complex_t MaterialBySLDImpl::refractiveIndex2(double wavelength) const
 {
     return 1.0 - getWlPrefactor(wavelength) * getSLD(m_sld, m_abs_term);
 }
 
-complex_t WavelengthIndependentMaterial::materialData() const
+complex_t MaterialBySLDImpl::materialData() const
 {
     return complex_t(m_sld, m_abs_term);
 }
 
-complex_t WavelengthIndependentMaterial::scalarSubtrSLD(const WavevectorInfo& wavevectors) const
+complex_t MaterialBySLDImpl::scalarSubtrSLD(const WavevectorInfo& wavevectors) const
 {
     double wavelength = wavevectors.getWavelength();
     return 1.0 / getWlPrefactor(wavelength) - getSLD(m_sld, m_abs_term);
 }
 
-void WavelengthIndependentMaterial::print(std::ostream& ostr) const
+void MaterialBySLDImpl::print(std::ostream& ostr) const
 {
-    ostr << "WavelengthIndependentMaterial:" << getName() << "<" << this << ">{ "
+    ostr << "MaterialBySLD:" << getName() << "<" << this << ">{ "
          << "sld=" << m_sld << ", absorp_term=" << m_abs_term
          << ", B=" << magnetization() << "}";
 }
diff --git a/Core/Material/WavelengthIndependentMaterial.h b/Core/Material/MaterialBySLDImpl.h
similarity index 73%
rename from Core/Material/WavelengthIndependentMaterial.h
rename to Core/Material/MaterialBySLDImpl.h
index 2192718a17bd9050852df97bb6e09ba6af316456..e7d3df780922eecabe33756eda7ad6d0f8387eb0 100644
--- a/Core/Material/WavelengthIndependentMaterial.h
+++ b/Core/Material/MaterialBySLDImpl.h
@@ -2,8 +2,8 @@
 //
 //  BornAgain: simulate and fit scattering at grazing incidence
 //
-//! @file      Core/Material/WavelengthIndependentMaterial.h
-//! @brief     Defines class WavelengthIndependentMaterial.
+//! @file      Core/Material/MaterialBySLDImpl.h
+//! @brief     Defines class MaterialBySLDImpl.
 //!
 //! @homepage  http://www.bornagainproject.org
 //! @license   GNU General Public License v3 or higher (see COPYING)
@@ -13,8 +13,8 @@
 //
 // ************************************************************************** //
 
-#ifndef WAVELENGTHINDEPENDENTMATERIAL_H_
-#define WAVELENGTHINDEPENDENTMATERIAL_H_
+#ifndef MATERIALBYSLDIMPL_H_
+#define MATERIALBYSLDIMPL_H_
 
 #include "MagneticMaterialImpl.h"
 #include "Material.h"
@@ -22,17 +22,15 @@
 //! Material implementation based on wavelength-independent data (valid for a range of wavelengths)
 //! @ingroup materials
 
-class BA_CORE_API_ WavelengthIndependentMaterial : public MagneticMaterialImpl
+class BA_CORE_API_ MaterialBySLDImpl : public MagneticMaterialImpl
 {
 public:
     friend BA_CORE_API_ Material MaterialBySLD(const std::string&, double, double, kvector_t);
 
-    friend BA_CORE_API_ Material MaterialByAbsCX(const std::string&, double, double, kvector_t);
-
-    virtual ~WavelengthIndependentMaterial();
+    virtual ~MaterialBySLDImpl() = default;
 
     //! Returns pointer to a copy of material
-    virtual WavelengthIndependentMaterial* clone() const override;
+    virtual MaterialBySLDImpl* clone() const override;
 
     //! Returns refractive index
     virtual complex_t refractiveIndex(double wavelength) const override;
@@ -46,7 +44,7 @@ public:
     //! Returns type of material implementation
     virtual MATERIAL_TYPES typeID() const override
     {
-        return MATERIAL_TYPES::WavelengthIndependentMaterial;
+        return MATERIAL_TYPES::MaterialBySLD;
     }
 
     //! Returns (\f$ \pi/\lambda^2 \f$ - sld), sld (in \f$nm^{-2}\f$) being the scattering length density
@@ -56,7 +54,7 @@ public:
     virtual void print(std::ostream &ostr) const override;
 
 private:
-    WavelengthIndependentMaterial(const std::string& name, double sld, double abs_term,
+    MaterialBySLDImpl(const std::string& name, double sld, double abs_term,
                                   kvector_t magnetization);
 
     double m_sld; //!< product of number density and coherent scattering length
@@ -65,4 +63,4 @@ private:
                         //!< absorption cross-section normalized to wavelength
 };
 
-#endif /* WAVELENGTHINDEPENDENTMATERIAL_H_ */
+#endif /* MATERIALBYSLDIMPL_H_ */
diff --git a/Core/Material/MaterialFactoryFuncs.cpp b/Core/Material/MaterialFactoryFuncs.cpp
index caa94c0ac7c2409204ab244d1594408b1b220206..74b97c3e3ad3b4bb56b6d4f353626bde43c11ae4 100644
--- a/Core/Material/MaterialFactoryFuncs.cpp
+++ b/Core/Material/MaterialFactoryFuncs.cpp
@@ -1,6 +1,15 @@
 #include "MaterialFactoryFuncs.h"
-#include "RefractiveCoefMaterial.h"
-#include "WavelengthIndependentMaterial.h"
+#include "MaterialBySLDImpl.h"
+#include "RefractiveMaterialImpl.h"
+#include "SlicedParticle.h"
+
+namespace {
+MATERIAL_TYPES checkMaterialTypes(const Material& layer_mat,
+                                  const std::vector<HomogeneousRegion>& regions);
+
+template <class T>
+T averageData(const Material& layer_mat, const std::vector<HomogeneousRegion>& regions, T (*func_ptr) (const Material&));
+}
 
 Material HomogeneousMaterial(const std::string& name, complex_t refractive_index,
                                  kvector_t magnetization)
@@ -13,8 +22,8 @@ Material HomogeneousMaterial(const std::string& name, complex_t refractive_index
 Material HomogeneousMaterial(const std::string& name, double delta, double beta,
                                  kvector_t magnetization)
 {
-    std::unique_ptr<RefractiveCoefMaterial> mat_impl(
-        new RefractiveCoefMaterial(name, delta, beta, magnetization));
+    std::unique_ptr<RefractiveMaterialImpl> mat_impl(
+        new RefractiveMaterialImpl(name, delta, beta, magnetization));
     return Material(std::move(mat_impl));
 }
 
@@ -26,8 +35,8 @@ Material HomogeneousMaterial()
 Material MaterialBySLD(const std::string& name, double sld, double abs_term,
                        kvector_t magnetization)
 {
-    std::unique_ptr<WavelengthIndependentMaterial> mat_impl(
-        new WavelengthIndependentMaterial(name, sld, abs_term, magnetization));
+    std::unique_ptr<MaterialBySLDImpl> mat_impl(
+        new MaterialBySLDImpl(name, sld, abs_term, magnetization));
     return Material(std::move(mat_impl));
 }
 
@@ -35,12 +44,72 @@ constexpr double basic_wavelength = 0.1798197; // nm, wavelength of 2200 m/s neu
 Material MaterialByAbsCX(const std::string& name, double sld, double abs_cx,
                          kvector_t magnetization)
 {
-    std::unique_ptr<WavelengthIndependentMaterial> mat_impl(
-        new WavelengthIndependentMaterial(name, sld, abs_cx / basic_wavelength, magnetization));
-    return Material(std::move(mat_impl));
+    return MaterialBySLD(name, sld, abs_cx / basic_wavelength, magnetization);
 }
 
 Material MaterialBySLD()
 {
     return MaterialBySLD("vacuum", 0.0, 0.0, kvector_t{});
 }
+
+Material createAveragedMaterial(const Material& layer_mat,
+                                const std::vector<HomogeneousRegion>& regions)
+{
+    // determine the type of returned material
+    const MATERIAL_TYPES avr_material_type = checkMaterialTypes(layer_mat, regions);
+
+    // create the name of returned material
+    const std::string avr_mat_name = layer_mat.getName() + "_avg";
+
+    // calculate averaged magnetization
+    kvector_t (*avrMag)(const Material&)
+        = [](const Material& mat) { return mat.magnetization(); };
+    const kvector_t mag_avr = averageData(layer_mat, regions, avrMag);
+
+    if (avr_material_type == MATERIAL_TYPES::RefractiveMaterial) {
+        // avrData returns (1 - mdc)^2 - 1, where mdc is material data conjugate
+        complex_t (*avrData)(const Material&) = [](const Material& mat) -> complex_t {
+            const complex_t mdc = std::conj(mat.materialData());
+            return mdc * mdc - 2.0 * mdc;
+        };
+        const complex_t avr_mat_data
+            = std::conj(1.0 - std::sqrt(1.0 + averageData(layer_mat, regions, avrData)));
+        return HomogeneousMaterial(avr_mat_name, avr_mat_data.real(), avr_mat_data.imag(), mag_avr);
+    } else if (avr_material_type == MATERIAL_TYPES::MaterialBySLD) {
+        complex_t (*avrData)(const Material&)
+            = [](const Material& mat) { return mat.materialData(); };
+        const complex_t avr_mat_data = averageData(layer_mat, regions, avrData);
+        return MaterialBySLD(avr_mat_name, avr_mat_data.real(), avr_mat_data.imag(), mag_avr);
+    } else
+        throw std::runtime_error("Error in CalculateAverageMaterial: unknown material type.");
+}
+
+namespace
+{
+MATERIAL_TYPES checkMaterialTypes(const Material& layer_mat, const std::vector<HomogeneousRegion>& regions)
+{
+    MATERIAL_TYPES result = layer_mat.typeID();
+    bool isDefault = layer_mat.isDefaultMaterial();
+    for (auto& region : regions) {
+        if (isDefault) {
+            result = region.m_material.typeID();
+            isDefault = region.m_material.isDefaultMaterial();
+            continue;
+        }
+        if (region.m_material.typeID() != result && !region.m_material.isDefaultMaterial())
+            throw std::runtime_error("Error in checkMaterialTypes(): non-default materials of "
+                                     "different types used simultaneously.");
+    }
+    return result;
+}
+
+template <class T>
+T averageData(const Material& layer_mat, const std::vector<HomogeneousRegion>& regions, T (*func_ptr) (const Material&))
+{
+    const T layer_data = func_ptr(layer_mat);
+    T average = layer_data;
+    for (auto& region : regions)
+        average += region.m_volume * (func_ptr(region.m_material) - layer_data);
+    return average;
+}
+}
diff --git a/Core/Material/MaterialFactoryFuncs.h b/Core/Material/MaterialFactoryFuncs.h
index 70507d2db0cfffddd0b576017361b24ecc7a93ec..7bbd47200756462a2d5d93aabb9a21b07a7c635e 100644
--- a/Core/Material/MaterialFactoryFuncs.h
+++ b/Core/Material/MaterialFactoryFuncs.h
@@ -18,6 +18,8 @@
 
 #include "Material.h"
 
+struct HomogeneousRegion;
+
 //! @ingroup materials
 
 //! Constructs a material with _name_ and _refractive_index_.
@@ -35,10 +37,7 @@ BA_CORE_API_ Material HomogeneousMaterial(const std::string& name, double delta,
 
 //! @ingroup materials
 
-//! Constructs vacuum material based on refractive coefficients.
-//! Though in practice there is no difference between vacuum materials
-//! produced with MaterialBySLD() and HomogeneousMaterial(), they are not equal because of
-//! the difference in the type of underlying data
+//! Constructs material with zero refractive coefficients and zero magnetization.
 BA_CORE_API_ Material HomogeneousMaterial();
 
 //! @ingroup materials
@@ -70,10 +69,13 @@ BA_CORE_API_ Material MaterialByAbsCX(const std::string& name, double sld, doubl
 
 //! @ingroup materials
 
-//! Constructs wavelength-independent vacuum material.
-//! Though in practice there is no difference between vacuum materials
-//! produced with MaterialBySLD() and HomogeneousMaterial(), they are not equal because of
-//! the difference in the type of underlying data
+//! Constructs wavelength-independent material with zero sld and zero magnetization.
 BA_CORE_API_ Material MaterialBySLD();
 
+//! @ingroup materials
+
+//! Creates averaged material. Square refractive index of returned material is arithmetic mean over
+//! _regions_ and _layer_mat_. Magnetization (if present) is averaged linearly.
+BA_CORE_API_ Material createAveragedMaterial(const Material& layer_mat, const std::vector<HomogeneousRegion>& regions);
+
 #endif /* MATERIALFACTORYFUNCS_H_ */
diff --git a/Core/Material/RefractiveCoefMaterial.cpp b/Core/Material/RefractiveMaterialImpl.cpp
similarity index 52%
rename from Core/Material/RefractiveCoefMaterial.cpp
rename to Core/Material/RefractiveMaterialImpl.cpp
index 8ef5f37b9460e77c5af283b94021737a7ecbdd1a..0779488d9c6e340a9a22e2373a12d89134625b05 100644
--- a/Core/Material/RefractiveCoefMaterial.cpp
+++ b/Core/Material/RefractiveMaterialImpl.cpp
@@ -1,47 +1,44 @@
-#include "RefractiveCoefMaterial.h"
+#include "RefractiveMaterialImpl.h"
 #include "WavevectorInfo.h"
 
-RefractiveCoefMaterial::RefractiveCoefMaterial(const std::string& name, double delta, double beta,
+RefractiveMaterialImpl::RefractiveMaterialImpl(const std::string& name, double delta, double beta,
                                                  kvector_t magnetization)
     : MagneticMaterialImpl(name, magnetization)
     , m_delta(delta)
     , m_beta(beta)
 {}
 
-RefractiveCoefMaterial::~RefractiveCoefMaterial()
-{}
-
-RefractiveCoefMaterial* RefractiveCoefMaterial::clone() const
+RefractiveMaterialImpl* RefractiveMaterialImpl::clone() const
 {
-    return new RefractiveCoefMaterial(*this);
+    return new RefractiveMaterialImpl(*this);
 }
 
-complex_t RefractiveCoefMaterial::refractiveIndex(double) const
+complex_t RefractiveMaterialImpl::refractiveIndex(double) const
 {
     return complex_t(1.0 - m_delta, m_beta);
 }
 
-complex_t RefractiveCoefMaterial::refractiveIndex2(double) const
+complex_t RefractiveMaterialImpl::refractiveIndex2(double) const
 {
     complex_t result(1.0 - m_delta, m_beta);
     return result * result;
 }
 
-complex_t RefractiveCoefMaterial::materialData() const
+complex_t RefractiveMaterialImpl::materialData() const
 {
     return complex_t(m_delta, m_beta);
 }
 
-complex_t RefractiveCoefMaterial::scalarSubtrSLD(const WavevectorInfo& wavevectors) const
+complex_t RefractiveMaterialImpl::scalarSubtrSLD(const WavevectorInfo& wavevectors) const
 {
     double wavelength = wavevectors.getWavelength();
     double prefactor = M_PI/wavelength/wavelength;
     return prefactor * refractiveIndex2(wavelength);
 }
 
-void RefractiveCoefMaterial::print(std::ostream& ostr) const
+void RefractiveMaterialImpl::print(std::ostream& ostr) const
 {
-    ostr << "RefractiveCoefMaterial:" << getName() << "<" << this << ">{ "
+    ostr << "RefractiveMaterial:" << getName() << "<" << this << ">{ "
          << "delta=" << m_delta << ", beta=" << m_beta
          << ", B=" << magnetization() << "}";
 }
diff --git a/Core/Material/RefractiveCoefMaterial.h b/Core/Material/RefractiveMaterialImpl.h
similarity index 82%
rename from Core/Material/RefractiveCoefMaterial.h
rename to Core/Material/RefractiveMaterialImpl.h
index 3611d948b1b166669213ba80f9c13d84a9abd9d8..56dba40b3968721db47e659b98e14ab28a866113 100644
--- a/Core/Material/RefractiveCoefMaterial.h
+++ b/Core/Material/RefractiveMaterialImpl.h
@@ -2,8 +2,8 @@
 //
 //  BornAgain: simulate and fit scattering at grazing incidence
 //
-//! @file      Core/Material/RefractiveCoefMaterial.h
-//! @brief     Defines class RefractiveCoefMaterial.
+//! @file      Core/Material/RefractiveMaterialImpl.h
+//! @brief     Defines class RefractiveMaterialImpl.
 //!
 //! @homepage  http://www.bornagainproject.org
 //! @license   GNU General Public License v3 or higher (see COPYING)
@@ -13,8 +13,8 @@
 //
 // ************************************************************************** //
 
-#ifndef REFRACTIVECOEFMATERIAL_H_
-#define REFRACTIVECOEFMATERIAL_H_
+#ifndef REFRACTIVEMATERIALIMPL_H_
+#define REFRACTIVEMATERIALIMPL_H_
 
 #include "MagneticMaterialImpl.h"
 #include "Material.h"
@@ -22,16 +22,16 @@
 //! Material implementation based on refractive coefficiencts (valid for one wavelength value only)
 //! @ingroup materials
 
-class BA_CORE_API_ RefractiveCoefMaterial : public MagneticMaterialImpl
+class BA_CORE_API_ RefractiveMaterialImpl : public MagneticMaterialImpl
 {
 public:
     friend BA_CORE_API_ Material HomogeneousMaterial(const std::string&, double, double,
                                                          kvector_t);
 
-    virtual ~RefractiveCoefMaterial();
+    virtual ~RefractiveMaterialImpl() = default;
 
     //! Returns pointer to a copy of material
-    virtual RefractiveCoefMaterial* clone() const override;
+    virtual RefractiveMaterialImpl* clone() const override;
 
     //! Returns refractive index
     //! For this particular implementation returned value does not depend
@@ -49,7 +49,7 @@ public:
     //! Returns type of material implementation
     virtual MATERIAL_TYPES typeID() const override
     {
-        return MATERIAL_TYPES::RefractiveCoefMaterial;
+        return MATERIAL_TYPES::RefractiveMaterial;
     }
 
     //! Returns (\f$ \pi/\lambda^2 \f$ - sld), sld (in \f$nm^{-2}\f$) being the scattering length density.
@@ -62,11 +62,11 @@ public:
     virtual void print(std::ostream &ostr) const override;
 
 private:
-    RefractiveCoefMaterial(const std::string& name, double delta, double beta,
+    RefractiveMaterialImpl(const std::string& name, double delta, double beta,
                            kvector_t magnetization);
 
     double m_delta; //!< \f$\delta\f$ coefficient for refractive index \f$n = 1 - \delta + i \beta\f$
     double m_beta; //!< \f$\beta\f$ coefficient for refractive index \f$n = 1 - \delta + i \beta\f$
 };
 
-#endif /* REFRACTIVECOEFMATERIAL_H_ */
+#endif /* REFRACTIVEMATERIALIMPL_H_ */
diff --git a/Core/Simulation/Simulation.cpp b/Core/Simulation/Simulation.cpp
index 189e980b7703411259c1980d1f85bf915a143327..67f83bc2041ead0c58747cde1689b1525438dbaf 100644
--- a/Core/Simulation/Simulation.cpp
+++ b/Core/Simulation/Simulation.cpp
@@ -215,7 +215,7 @@ void Simulation::runSingleSimulation()
     if (m_options.getNumberOfThreads() == 1) {
         // Single thread.
         std::unique_ptr<MainComputation> P_computation(
-            new MainComputation(*sample(), m_instrument, m_options, m_progress, batch_start, batch_end));
+            new MainComputation(*sample(), m_options, m_progress, batch_start, batch_end));
         P_computation->run(); // the work is done here
         if (!P_computation->isCompleted()) {
             std::string message = P_computation->errorMessage();
@@ -246,7 +246,7 @@ void Simulation::runSingleSimulation()
             else
                 end_it = batch_start + end_thread_index;
             computations.emplace_back(
-                new MainComputation(*sample(), m_instrument, m_options, m_progress, begin_it, end_it));
+                new MainComputation(*sample(), m_options, m_progress, begin_it, end_it));
         }
 
         // Run simulations in n threads.
diff --git a/Tests/UnitTests/Core/Other/MaterialTest.h b/Tests/UnitTests/Core/Other/MaterialTest.h
index d7fc1fd16f3eca0ee8fade5f62b1225f2ac6bf57..7855206827a6e7bee54c03c2c46835059a771ae2 100644
--- a/Tests/UnitTests/Core/Other/MaterialTest.h
+++ b/Tests/UnitTests/Core/Other/MaterialTest.h
@@ -1,6 +1,7 @@
 #include "MaterialFactoryFuncs.h"
-#include "WavelengthIndependentMaterial.h"
-#include "RefractiveCoefMaterial.h"
+#include "MaterialBySLDImpl.h"
+#include "RefractiveMaterialImpl.h"
+#include "SlicedParticle.h"
 #include "WavevectorInfo.h"
 #include "Rotations.h"
 #include "Units.h"
@@ -59,24 +60,14 @@ TEST_F(MaterialTest, MaterialConstruction)
     EXPECT_EQ(material_data, material6.materialData());
     EXPECT_EQ(default_magnetism, material6.magnetization());
 
-    Material material7 = HomogeneousMaterial();
-    EXPECT_EQ(material7.getName(), HomogeneousMaterial().getName());
-    EXPECT_EQ(material7.getName(), MaterialBySLD().getName());
-    EXPECT_EQ(material7.materialData(), HomogeneousMaterial().materialData());
-    EXPECT_EQ(material7.materialData(), MaterialBySLD().materialData());
-    EXPECT_EQ(material7.magnetization(), HomogeneousMaterial().magnetization());
-    EXPECT_EQ(material7.magnetization(), MaterialBySLD().magnetization());
-    EXPECT_TRUE(material7.typeID() == HomogeneousMaterial().typeID());
-    EXPECT_FALSE(material7.typeID() == MaterialBySLD().typeID());
-
     constexpr double basic_wavelength = 0.1798197; // nm
-    Material material8 = MaterialByAbsCX("Material", material_data.real(),
+    Material material7 = MaterialByAbsCX("Material", material_data.real(),
                                          material_data.imag() * basic_wavelength);
-    EXPECT_TRUE(material8.getName() == material6.getName());
-    EXPECT_TRUE(material8.magnetization() == material6.magnetization());
-    EXPECT_DOUBLE_EQ(material8.materialData().real(), material6.materialData().real());
-    EXPECT_DOUBLE_EQ(material8.materialData().imag(), material6.materialData().imag());
-    EXPECT_TRUE(material8.typeID() == material6.typeID());
+    EXPECT_TRUE(material7.getName() == material6.getName());
+    EXPECT_TRUE(material7.magnetization() == material6.magnetization());
+    EXPECT_DOUBLE_EQ(material7.materialData().real(), material6.materialData().real());
+    EXPECT_DOUBLE_EQ(material7.materialData().imag(), material6.materialData().imag());
+    EXPECT_TRUE(material7.typeID() == material6.typeID());
 }
 
 TEST_F(MaterialTest, MaterialTransform)
@@ -102,6 +93,28 @@ TEST_F(MaterialTest, MaterialTransform)
     EXPECT_EQ(transformed_mag, material4.magnetization());
 }
 
+TEST_F(MaterialTest, DefaultMaterials)
+{
+    Material material = HomogeneousMaterial();
+    const double dummy_wavelength = 1.0;
+
+    EXPECT_EQ(material.getName(), std::string("vacuum"));
+    EXPECT_EQ(material.getName(), MaterialBySLD().getName());
+
+    EXPECT_EQ(material.materialData(), complex_t());
+    EXPECT_EQ(material.materialData(), MaterialBySLD().materialData());
+
+    EXPECT_EQ(material.magnetization(), kvector_t{});
+    EXPECT_EQ(material.magnetization(), MaterialBySLD().magnetization());
+
+    EXPECT_EQ(material.refractiveIndex(dummy_wavelength), complex_t(1.0, 0.0));
+    EXPECT_EQ(material.refractiveIndex(dummy_wavelength),
+              MaterialBySLD().refractiveIndex(dummy_wavelength));
+
+    EXPECT_TRUE(material.typeID() == HomogeneousMaterial().typeID());
+    EXPECT_FALSE(material.typeID() == MaterialBySLD().typeID());
+}
+
 TEST_F(MaterialTest, ComputationTest)
 {
     // Reference data for Fe taken from
@@ -140,12 +153,47 @@ TEST_F(MaterialTest, ComputationTest)
     EXPECT_FLOAT_EQ(subtrSLD.imag(), subtrSLDWlIndep.imag());
 }
 
+TEST_F(MaterialTest, AveragedMaterialTest)
+{
+    kvector_t magnetization = kvector_t {1.0, 0.0, 0.0};
+    const Material material = HomogeneousMaterial("Material", 0.5, 0.5, magnetization);
+    const std::vector<HomogeneousRegion> regions = {HomogeneousRegion{0.25, material},
+                                                    HomogeneousRegion{0.25, material}};
+
+    const Material material_avr = createAveragedMaterial(material, regions);
+    EXPECT_EQ(material_avr.getName(), material.getName() + "_avg");
+    EXPECT_EQ(material_avr.magnetization(), magnetization);
+    EXPECT_DOUBLE_EQ(material_avr.materialData().real(), 0.5);
+    EXPECT_DOUBLE_EQ(material_avr.materialData().imag(), 0.5);
+    EXPECT_TRUE(material_avr.typeID() == MATERIAL_TYPES::RefractiveMaterial);
+
+    const Material material2 = MaterialBySLD();
+    const Material material_avr2 = createAveragedMaterial(material2, regions);
+    const complex_t expected_res = std::conj(1.0 - std::sqrt(complex_t(0.5, 0.25)));
+    EXPECT_DOUBLE_EQ(material_avr2.materialData().real(), expected_res.real());
+    EXPECT_DOUBLE_EQ(material_avr2.materialData().imag(), expected_res.imag());
+    EXPECT_EQ(material_avr2.magnetization(), kvector_t(0.5, 0.0, 0.0));
+    EXPECT_TRUE(material_avr2.typeID() == MATERIAL_TYPES::RefractiveMaterial);
+
+    const Material material3 = MaterialBySLD("Material3", 0.5, 0.5, magnetization);
+    EXPECT_THROW(createAveragedMaterial(material3, regions), std::runtime_error);
+
+    const Material material4 = HomogeneousMaterial();
+    const std::vector<HomogeneousRegion> regions2
+        = {HomogeneousRegion{0.25, material3}, HomogeneousRegion{0.25, material3}};
+    const Material material_avr3 = createAveragedMaterial(material4, regions2);
+    EXPECT_DOUBLE_EQ(material_avr3.materialData().real(), 0.25);
+    EXPECT_DOUBLE_EQ(material_avr3.materialData().imag(), 0.25);
+    EXPECT_EQ(material_avr3.magnetization(), kvector_t(0.5, 0.0, 0.0));
+    EXPECT_TRUE(material_avr3.typeID() == MATERIAL_TYPES::MaterialBySLD);
+}
+
 TEST_F(MaterialTest, TypeIdsTest)
 {
     Material material = MaterialBySLD("Material", 1.0, 1.0);
     Material material2 = HomogeneousMaterial("Material", 1.0, 1.0);
-    EXPECT_TRUE(material.typeID() == MATERIAL_TYPES::WavelengthIndependentMaterial);
-    EXPECT_TRUE(material2.typeID() == MATERIAL_TYPES::RefractiveCoefMaterial);
+    EXPECT_TRUE(material.typeID() == MATERIAL_TYPES::MaterialBySLD);
+    EXPECT_TRUE(material2.typeID() == MATERIAL_TYPES::RefractiveMaterial);
     EXPECT_TRUE(material.typeID() != material2.typeID());
     Material material3 = MaterialBySLD("Material", 1.0, 1.0);
     EXPECT_TRUE(material.typeID() == material3.typeID());
diff --git a/auto/Wrap/doxygen_core.i b/auto/Wrap/doxygen_core.i
index fe1e35dfc8e589b62e0cdc8c5e1b217beb31769a..b3e530f50e2120cce655dc20d577d30a57c991ce 100644
--- a/auto/Wrap/doxygen_core.i
+++ b/auto/Wrap/doxygen_core.i
@@ -43,19 +43,14 @@ C++ includes: AdjustMinimizerStrategy.h
 
 
 // File: classBaseMaterialImpl.xml
-%feature("docstring") BaseMaterialImpl "
+%feature("docstring") BaseMaterialImpl "";
 
-A basic implementation for homogeneous and wavelength-independent materials. Incorporates data and methods required to handle material magnetization.
+%feature("docstring")  BaseMaterialImpl::BaseMaterialImpl "BaseMaterialImpl::BaseMaterialImpl(const std::string &name)
 
-C++ includes: BaseMaterialImpl.h
+Constructs basic material with name. 
 ";
 
-%feature("docstring")  BaseMaterialImpl::BaseMaterialImpl "BaseMaterialImpl::BaseMaterialImpl(const std::string &name, kvector_t magnetization)
-
-Constructs basic material with name and magnetization. 
-";
-
-%feature("docstring")  BaseMaterialImpl::~BaseMaterialImpl "BaseMaterialImpl::~BaseMaterialImpl()
+%feature("docstring")  BaseMaterialImpl::~BaseMaterialImpl "virtual BaseMaterialImpl::~BaseMaterialImpl()=default
 ";
 
 %feature("docstring")  BaseMaterialImpl::clone "virtual BaseMaterialImpl* BaseMaterialImpl::clone() const =0
@@ -63,26 +58,30 @@ Constructs basic material with name and magnetization.
 Returns pointer to a copy of material. 
 ";
 
-%feature("docstring")  BaseMaterialImpl::inverted "BaseMaterialImpl * BaseMaterialImpl::inverted() const
+%feature("docstring")  BaseMaterialImpl::inverted "virtual BaseMaterialImpl* BaseMaterialImpl::inverted() const =0
 
 Constructs a material with inverted magnetization. 
 ";
 
 %feature("docstring")  BaseMaterialImpl::refractiveIndex "virtual complex_t BaseMaterialImpl::refractiveIndex(double wavelength) const =0
+
+Returns refractive index. 
 ";
 
 %feature("docstring")  BaseMaterialImpl::refractiveIndex2 "virtual complex_t BaseMaterialImpl::refractiveIndex2(double wavelength) const =0
+
+Returns squared refractive index. 
 ";
 
-%feature("docstring")  BaseMaterialImpl::isScalarMaterial "bool BaseMaterialImpl::isScalarMaterial() const
+%feature("docstring")  BaseMaterialImpl::isScalarMaterial "virtual bool BaseMaterialImpl::isScalarMaterial() const =0
 
 Indicates whether the interaction with the material is scalar. This means that different polarization states will be diffracted equally 
 ";
 
-%feature("docstring")  BaseMaterialImpl::isMagneticMaterial "bool BaseMaterialImpl::isMagneticMaterial() const
+%feature("docstring")  BaseMaterialImpl::isMagneticMaterial "virtual bool BaseMaterialImpl::isMagneticMaterial() const =0
 ";
 
-%feature("docstring")  BaseMaterialImpl::magnetization "kvector_t BaseMaterialImpl::magnetization() const
+%feature("docstring")  BaseMaterialImpl::magnetization "virtual kvector_t BaseMaterialImpl::magnetization() const =0
 
 Returns the magnetization (in A/m) 
 ";
@@ -92,17 +91,22 @@ Returns the magnetization (in A/m)
 Returns underlying material data. 
 ";
 
+%feature("docstring")  BaseMaterialImpl::typeID "virtual MATERIAL_TYPES BaseMaterialImpl::typeID() const =0
+
+Returns type of material implementation. 
+";
+
 %feature("docstring")  BaseMaterialImpl::scalarSubtrSLD "virtual complex_t BaseMaterialImpl::scalarSubtrSLD(const WavevectorInfo &wavevectors) const =0
 
 Returns (  $ \\\\pi/\\\\lambda^2 $ - sld), sld being the scattering length density. 
 ";
 
-%feature("docstring")  BaseMaterialImpl::polarizedSubtrSLD "Eigen::Matrix2cd BaseMaterialImpl::polarizedSubtrSLD(const WavevectorInfo &wavevectors) const
+%feature("docstring")  BaseMaterialImpl::polarizedSubtrSLD "virtual Eigen::Matrix2cd BaseMaterialImpl::polarizedSubtrSLD(const WavevectorInfo &wavevectors) const =0
 
 Returns (  $ \\\\pi/\\\\lambda^2 $ - sld) matrix with magnetization corrections. 
 ";
 
-%feature("docstring")  BaseMaterialImpl::transformedMaterial "BaseMaterialImpl * BaseMaterialImpl::transformedMaterial(const Transform3D &transform) const
+%feature("docstring")  BaseMaterialImpl::transformedMaterial "virtual BaseMaterialImpl* BaseMaterialImpl::transformedMaterial(const Transform3D &transform) const =0
 ";
 
 %feature("docstring")  BaseMaterialImpl::print "virtual void BaseMaterialImpl::print(std::ostream &ostr) const =0
@@ -597,21 +601,21 @@ C++ includes: BoxCompositionBuilder.h
 ";
 
 
-// File: structIntegratorReal_1_1CallBackHolder.xml
-%feature("docstring") IntegratorReal::CallBackHolder "
+// File: structIntegratorMCMiser_1_1CallBackHolder.xml
+%feature("docstring") IntegratorMCMiser::CallBackHolder "
 
 structure holding the object and possible extra parameters
 
-C++ includes: IntegratorReal.h
+C++ includes: IntegratorMCMiser.h
 ";
 
 
-// File: structIntegratorMCMiser_1_1CallBackHolder.xml
-%feature("docstring") IntegratorMCMiser::CallBackHolder "
+// File: structIntegratorReal_1_1CallBackHolder.xml
+%feature("docstring") IntegratorReal::CallBackHolder "
 
 structure holding the object and possible extra parameters
 
-C++ includes: IntegratorMCMiser.h
+C++ includes: IntegratorReal.h
 ";
 
 
@@ -9882,6 +9886,54 @@ C++ includes: MagneticParticlesBuilder.h
 ";
 
 
+// File: classMagneticMaterialImpl.xml
+%feature("docstring") MagneticMaterialImpl "
+
+Basic implementation for magnetized material. Incorporates data and methods required to handle material magnetization.
+
+C++ includes: MagneticMaterialImpl.h
+";
+
+%feature("docstring")  MagneticMaterialImpl::MagneticMaterialImpl "MagneticMaterialImpl::MagneticMaterialImpl(const std::string &name, kvector_t magnetization)
+
+Constructs basic material with name and magnetization. 
+";
+
+%feature("docstring")  MagneticMaterialImpl::~MagneticMaterialImpl "virtual MagneticMaterialImpl::~MagneticMaterialImpl()=default
+";
+
+%feature("docstring")  MagneticMaterialImpl::clone "virtual MagneticMaterialImpl* MagneticMaterialImpl::clone() const override=0
+
+Returns pointer to a copy of material. 
+";
+
+%feature("docstring")  MagneticMaterialImpl::inverted "MagneticMaterialImpl * MagneticMaterialImpl::inverted() const override final
+
+Constructs a material with inverted magnetization. 
+";
+
+%feature("docstring")  MagneticMaterialImpl::isScalarMaterial "bool MagneticMaterialImpl::isScalarMaterial() const override final
+
+Indicates whether the interaction with the material is scalar. This means that different polarization states will be diffracted equally 
+";
+
+%feature("docstring")  MagneticMaterialImpl::isMagneticMaterial "bool MagneticMaterialImpl::isMagneticMaterial() const override final
+";
+
+%feature("docstring")  MagneticMaterialImpl::magnetization "kvector_t MagneticMaterialImpl::magnetization() const override final
+
+Returns the magnetization (in A/m) 
+";
+
+%feature("docstring")  MagneticMaterialImpl::polarizedSubtrSLD "Eigen::Matrix2cd MagneticMaterialImpl::polarizedSubtrSLD(const WavevectorInfo &wavevectors) const override final
+
+Returns (  $ \\\\pi/\\\\lambda^2 $ - sld) matrix with magnetization corrections. 
+";
+
+%feature("docstring")  MagneticMaterialImpl::transformedMaterial "MagneticMaterialImpl * MagneticMaterialImpl::transformedMaterial(const Transform3D &transform) const override final
+";
+
+
 // File: classMagneticParticleZeroFieldBuilder.xml
 %feature("docstring") MagneticParticleZeroFieldBuilder "
 
@@ -9952,7 +10004,7 @@ Controlled by the multi-threading machinery in  Simulation::runSingleSimulation(
 C++ includes: MainComputation.h
 ";
 
-%feature("docstring")  MainComputation::MainComputation "MainComputation::MainComputation(const MultiLayer &multilayer, const Instrument &instrument, const SimulationOptions &options, ProgressHandler &progress, const std::vector< SimulationElement >::iterator &begin_it, const std::vector< SimulationElement >::iterator &end_it)
+%feature("docstring")  MainComputation::MainComputation "MainComputation::MainComputation(const MultiLayer &multilayer, const SimulationOptions &options, ProgressHandler &progress, const std::vector< SimulationElement >::iterator &begin_it, const std::vector< SimulationElement >::iterator &end_it)
 ";
 
 %feature("docstring")  MainComputation::~MainComputation "MainComputation::~MainComputation()
@@ -10019,9 +10071,9 @@ Indicates whether the interaction with the material is scalar. This means that d
 Returns the name of material. 
 ";
 
-%feature("docstring")  Material::dataType "size_t Material::dataType() const
+%feature("docstring")  Material::typeID "MATERIAL_TYPES Material::typeID() const
 
-Returns hash code of underlying material implementation. 
+Returns the type of underlying material implementation. 
 ";
 
 %feature("docstring")  Material::magnetization "kvector_t Material::magnetization() const
@@ -10039,6 +10091,11 @@ Returns underlying material data.
 Returns true if material underlying data is nullptr. 
 ";
 
+%feature("docstring")  Material::isDefaultMaterial "bool Material::isDefaultMaterial() const
+
+Returns true if material has refractive index of (1.0, 0.0) and zero magnetization. 
+";
+
 %feature("docstring")  Material::scalarSubtrSLD "complex_t Material::scalarSubtrSLD(const WavevectorInfo &wavevectors) const
 
 Returns (  $ \\\\pi/\\\\lambda^2 $ - sld), sld (in  $nm^{-2}$) being the scattering length density. 
@@ -10053,6 +10110,53 @@ Returns (  $ \\\\pi/\\\\lambda^2 $ - sld) matrix with magnetization corrections.
 ";
 
 
+// File: classMaterialBySLDImpl.xml
+%feature("docstring") MaterialBySLDImpl "
+
+Material implementation based on wavelength-independent data (valid for a range of wavelengths)
+
+C++ includes: MaterialBySLDImpl.h
+";
+
+%feature("docstring")  MaterialBySLDImpl::~MaterialBySLDImpl "virtual MaterialBySLDImpl::~MaterialBySLDImpl()=default
+";
+
+%feature("docstring")  MaterialBySLDImpl::clone "MaterialBySLDImpl * MaterialBySLDImpl::clone() const override
+
+Returns pointer to a copy of material. 
+";
+
+%feature("docstring")  MaterialBySLDImpl::refractiveIndex "complex_t MaterialBySLDImpl::refractiveIndex(double wavelength) const override
+
+Returns refractive index. 
+";
+
+%feature("docstring")  MaterialBySLDImpl::refractiveIndex2 "complex_t MaterialBySLDImpl::refractiveIndex2(double wavelength) const override
+
+Returns squared refractive index. 
+";
+
+%feature("docstring")  MaterialBySLDImpl::materialData "complex_t MaterialBySLDImpl::materialData() const override
+
+Returns underlying material data. 
+";
+
+%feature("docstring")  MaterialBySLDImpl::typeID "virtual MATERIAL_TYPES MaterialBySLDImpl::typeID() const override
+
+Returns type of material implementation. 
+";
+
+%feature("docstring")  MaterialBySLDImpl::scalarSubtrSLD "complex_t MaterialBySLDImpl::scalarSubtrSLD(const WavevectorInfo &wavevectors) const override
+
+Returns (  $ \\\\pi/\\\\lambda^2 $ - sld), sld (in  $nm^{-2}$) being the scattering length density. 
+";
+
+%feature("docstring")  MaterialBySLDImpl::print "void MaterialBySLDImpl::print(std::ostream &ostr) const override
+
+Prints object data. 
+";
+
+
 // File: classMatrixFresnelMap.xml
 %feature("docstring") MatrixFresnelMap "
 
@@ -12127,43 +12231,48 @@ C++ includes: ParaCrystalBuilder.h
 ";
 
 
-// File: classRefractiveCoefMaterial.xml
-%feature("docstring") RefractiveCoefMaterial "
+// File: classRefractiveMaterialImpl.xml
+%feature("docstring") RefractiveMaterialImpl "
 
 Material implementation based on refractive coefficiencts (valid for one wavelength value only)
 
-C++ includes: RefractiveCoefMaterial.h
+C++ includes: RefractiveMaterialImpl.h
 ";
 
-%feature("docstring")  RefractiveCoefMaterial::~RefractiveCoefMaterial "RefractiveCoefMaterial::~RefractiveCoefMaterial()
+%feature("docstring")  RefractiveMaterialImpl::~RefractiveMaterialImpl "virtual RefractiveMaterialImpl::~RefractiveMaterialImpl()=default
 ";
 
-%feature("docstring")  RefractiveCoefMaterial::clone "RefractiveCoefMaterial * RefractiveCoefMaterial::clone() const override
+%feature("docstring")  RefractiveMaterialImpl::clone "RefractiveMaterialImpl * RefractiveMaterialImpl::clone() const override
 
 Returns pointer to a copy of material. 
 ";
 
-%feature("docstring")  RefractiveCoefMaterial::refractiveIndex "complex_t RefractiveCoefMaterial::refractiveIndex(double wavelength) const override
+%feature("docstring")  RefractiveMaterialImpl::refractiveIndex "complex_t RefractiveMaterialImpl::refractiveIndex(double wavelength) const override
 
 Returns refractive index For this particular implementation returned value does not depend on passed wavelength 
 ";
 
-%feature("docstring")  RefractiveCoefMaterial::refractiveIndex2 "complex_t RefractiveCoefMaterial::refractiveIndex2(double wavelength) const override
+%feature("docstring")  RefractiveMaterialImpl::refractiveIndex2 "complex_t RefractiveMaterialImpl::refractiveIndex2(double wavelength) const override
 
 Returns squared refractive index. For this particular implementation returned value does not depend on passed wavelength. 
 ";
 
-%feature("docstring")  RefractiveCoefMaterial::materialData "complex_t RefractiveCoefMaterial::materialData() const override
+%feature("docstring")  RefractiveMaterialImpl::materialData "complex_t RefractiveMaterialImpl::materialData() const override
 
 Returns underlying material data. 
 ";
 
-%feature("docstring")  RefractiveCoefMaterial::scalarSubtrSLD "complex_t RefractiveCoefMaterial::scalarSubtrSLD(const WavevectorInfo &wavevectors) const override
+%feature("docstring")  RefractiveMaterialImpl::typeID "virtual MATERIAL_TYPES RefractiveMaterialImpl::typeID() const override
+
+Returns type of material implementation. 
+";
+
+%feature("docstring")  RefractiveMaterialImpl::scalarSubtrSLD "complex_t RefractiveMaterialImpl::scalarSubtrSLD(const WavevectorInfo &wavevectors) const override
 
 Returns (  $ \\\\pi/\\\\lambda^2 $ - sld), sld (in  $nm^{-2}$) being the scattering length density. If the wavelength associated with passed wavevector is different from the one associated with refractive coefficients used during the object construction, provided result is inconsistent. 
 ";
 
-%feature("docstring")  RefractiveCoefMaterial::print "void RefractiveCoefMaterial::print(std::ostream &ostr) const override
+%feature("docstring")  RefractiveMaterialImpl::print "void RefractiveMaterialImpl::print(std::ostream &ostr) const override
 
 Prints object data. 
 ";
@@ -14218,48 +14327,6 @@ Returns true if area defined by two bins is inside or on border of polygon (more
 ";
 
 
-// File: classWavelengthIndependentMaterial.xml
-%feature("docstring") WavelengthIndependentMaterial "
-
-Material implementation based on wavelength-independent data (valid for a range of wavelengths)
-
-C++ includes: WavelengthIndependentMaterial.h
-";
-
-%feature("docstring")  WavelengthIndependentMaterial::~WavelengthIndependentMaterial "WavelengthIndependentMaterial::~WavelengthIndependentMaterial()
-";
-
-%feature("docstring")  WavelengthIndependentMaterial::clone "WavelengthIndependentMaterial * WavelengthIndependentMaterial::clone() const override
-
-Returns pointer to a copy of material. 
-";
-
-%feature("docstring")  WavelengthIndependentMaterial::refractiveIndex "complex_t WavelengthIndependentMaterial::refractiveIndex(double wavelength) const override
-
-Returns refractive index. 
-";
-
-%feature("docstring")  WavelengthIndependentMaterial::refractiveIndex2 "complex_t WavelengthIndependentMaterial::refractiveIndex2(double wavelength) const override
-
-Returns squared refractive index. 
-";
-
-%feature("docstring")  WavelengthIndependentMaterial::materialData "complex_t WavelengthIndependentMaterial::materialData() const override
-
-Returns underlying material data. 
-";
-
-%feature("docstring")  WavelengthIndependentMaterial::scalarSubtrSLD "complex_t WavelengthIndependentMaterial::scalarSubtrSLD(const WavevectorInfo &wavevectors) const override
-
-Returns (  $ \\\\pi/\\\\lambda^2 $ - sld), sld (in  $nm^{-2}$) being the scattering length density. 
-";
-
-%feature("docstring")  WavelengthIndependentMaterial::print "void WavelengthIndependentMaterial::print(std::ostream &ostr) const override
-
-Prints object data. 
-";
-
-
 // File: classWavevectorInfo.xml
 %feature("docstring") WavevectorInfo "
 
@@ -14335,52 +14402,52 @@ C++ includes: ZLimits.h
 // File: namespace_0D249.xml
 
 
-// File: namespace_0D275.xml
+// File: namespace_0D276.xml
 
 
-// File: namespace_0D281.xml
+// File: namespace_0D280.xml
 
 
-// File: namespace_0D285.xml
+// File: namespace_0D282.xml
 
 
-// File: namespace_0D291.xml
+// File: namespace_0D284.xml
 
 
-// File: namespace_0D306.xml
+// File: namespace_0D292.xml
 
 
-// File: namespace_0D314.xml
+// File: namespace_0D307.xml
 
 
-// File: namespace_0D320.xml
+// File: namespace_0D315.xml
 
 
-// File: namespace_0D323.xml
+// File: namespace_0D321.xml
 
 
-// File: namespace_0D325.xml
+// File: namespace_0D324.xml
 
 
-// File: namespace_0D346.xml
+// File: namespace_0D326.xml
 
 
-// File: namespace_0D355.xml
+// File: namespace_0D347.xml
 
 
-// File: namespace_0D388.xml
+// File: namespace_0D356.xml
 
 
-// File: namespace_0D395.xml
+// File: namespace_0D389.xml
 
 
-// File: namespace_0D501.xml
+// File: namespace_0D396.xml
 
 
-// File: namespace_0D523.xml
+// File: namespace_0D502.xml
 
 
-// File: namespace_0D57.xml
+// File: namespace_0D524.xml
 
 
 // File: namespace_0D63.xml
@@ -15969,10 +16036,13 @@ make Swappable
 // File: Rectangle_8h.xml
 
 
-// File: BaseMaterialImpl_8cpp.xml
+// File: BaseMaterialImpl_8h.xml
 
 
-// File: BaseMaterialImpl_8h.xml
+// File: MagneticMaterialImpl_8cpp.xml
+
+
+// File: MagneticMaterialImpl_8h.xml
 
 
 // File: Material_8cpp.xml
@@ -15981,6 +16051,12 @@ make Swappable
 // File: Material_8h.xml
 
 
+// File: MaterialBySLDImpl_8cpp.xml
+
+
+// File: MaterialBySLDImpl_8h.xml
+
+
 // File: MaterialFactoryFuncs_8cpp.xml
 %feature("docstring")  HomogeneousMaterial "Material HomogeneousMaterial(const std::string &name, complex_t refractive_index, kvector_t magnetization)
 
@@ -16006,7 +16082,7 @@ magnetization (in A/m)
 
 %feature("docstring")  HomogeneousMaterial "Material HomogeneousMaterial()
 
-Constructs vacuum material based on refractive coefficients. Though in practice there is no difference between vacuum materials produced with  MaterialBySLD() and  HomogeneousMaterial(), they are not equal because of the difference in the type of underlying data 
+Constructs material with zero refractive coefficients and zero magnetization. 
 ";
 
 %feature("docstring")  MaterialBySLD "Material MaterialBySLD(const std::string &name, double sld, double abs_term, kvector_t magnetization)
@@ -16045,7 +16121,12 @@ magnetization (in A/m)
 
 %feature("docstring")  MaterialBySLD "Material MaterialBySLD()
 
-Constructs wavelength-independent vacuum material. Though in practice there is no difference between vacuum materials produced with  MaterialBySLD() and  HomogeneousMaterial(), they are not equal because of the difference in the type of underlying data 
+Constructs wavelength-independent material with zero sld and zero magnetization. 
+";
+
+%feature("docstring")  createAveragedMaterial "Material createAveragedMaterial(const Material &layer_mat, const std::vector< HomogeneousRegion > &regions)
+
+Creates averaged material. Square refractive index of returned material is arithmetic mean over  regions and  layer_mat. Magnetization (if present) is averaged linearly. 
 ";
 
 
@@ -16074,7 +16155,7 @@ magnetization (in A/m)
 
 %feature("docstring")  HomogeneousMaterial "BA_CORE_API_ Material HomogeneousMaterial()
 
-Constructs vacuum material based on refractive coefficients. Though in practice there is no difference between vacuum materials produced with  MaterialBySLD() and  HomogeneousMaterial(), they are not equal because of the difference in the type of underlying data 
+Constructs material with zero refractive coefficients and zero magnetization. 
 ";
 
 %feature("docstring")  MaterialBySLD "BA_CORE_API_ Material MaterialBySLD(const std::string &name, double sld, double abs_term, kvector_t magnetization=kvector_t())
@@ -16113,7 +16194,12 @@ magnetization (in A/m)
 
 %feature("docstring")  MaterialBySLD "BA_CORE_API_ Material MaterialBySLD()
 
-Constructs wavelength-independent vacuum material. Though in practice there is no difference between vacuum materials produced with  MaterialBySLD() and  HomogeneousMaterial(), they are not equal because of the difference in the type of underlying data 
+Constructs wavelength-independent material with zero sld and zero magnetization. 
+";
+
+%feature("docstring")  createAveragedMaterial "BA_CORE_API_ Material createAveragedMaterial(const Material &layer_mat, const std::vector< HomogeneousRegion > &regions)
+
+Creates averaged material. Square refractive index of returned material is arithmetic mean over  regions and  layer_mat. Magnetization (if present) is averaged linearly. 
 ";
 
 
@@ -16125,16 +16211,10 @@ Constructs wavelength-independent vacuum material. Though in practice there is n
 // File: MaterialUtils_8h.xml
 
 
-// File: RefractiveCoefMaterial_8cpp.xml
-
-
-// File: RefractiveCoefMaterial_8h.xml
-
-
-// File: WavelengthIndependentMaterial_8cpp.xml
+// File: RefractiveMaterialImpl_8cpp.xml
 
 
-// File: WavelengthIndependentMaterial_8h.xml
+// File: RefractiveMaterialImpl_8h.xml
 
 
 // File: DecouplingApproximationStrategy_8cpp.xml
diff --git a/auto/Wrap/libBornAgainCore.py b/auto/Wrap/libBornAgainCore.py
index 043c43a7f1c1d4f87521259e92159d0f0bbe6cce..9a6f89fcf3b852377da5a47ee34f7676ba0113c9 100644
--- a/auto/Wrap/libBornAgainCore.py
+++ b/auto/Wrap/libBornAgainCore.py
@@ -22713,6 +22713,18 @@ class Material(_object):
         return _libBornAgainCore.Material_isEmpty(self)
 
 
+    def isDefaultMaterial(self):
+        """
+        isDefaultMaterial(Material self) -> bool
+
+        bool Material::isDefaultMaterial() const
+
+        Returns true if material has refractive index of (1.0, 0.0) and zero magnetization. 
+
+        """
+        return _libBornAgainCore.Material_isDefaultMaterial(self)
+
+
     def scalarSubtrSLD(self, wavevectors):
         """
         scalarSubtrSLD(Material self, WavevectorInfo wavevectors) -> complex_t
@@ -22750,7 +22762,7 @@ def HomogeneousMaterial(*args):
 
     BA_CORE_API_ Material HomogeneousMaterial()
 
-    Constructs vacuum material based on refractive coefficients. Though in practice there is no difference between vacuum materials produced with  MaterialBySLD() and  HomogeneousMaterial(), they are not equal because of the difference in the type of underlying data 
+    Constructs material with zero refractive coefficients and zero magnetization. 
 
     """
     return _libBornAgainCore.HomogeneousMaterial(*args)
@@ -22787,10 +22799,21 @@ def MaterialBySLD(*args):
 
     BA_CORE_API_ Material MaterialBySLD()
 
-    Constructs wavelength-independent vacuum material. Though in practice there is no difference between vacuum materials produced with  MaterialBySLD() and  HomogeneousMaterial(), they are not equal because of the difference in the type of underlying data 
+    Constructs wavelength-independent material with zero sld and zero magnetization. 
 
     """
     return _libBornAgainCore.MaterialBySLD(*args)
+
+def createAveragedMaterial(layer_mat, regions):
+    """
+    createAveragedMaterial(Material layer_mat, std::vector< HomogeneousRegion,std::allocator< HomogeneousRegion > > const & regions) -> Material
+
+    BA_CORE_API_ Material createAveragedMaterial(const Material &layer_mat, const std::vector< HomogeneousRegion > &regions)
+
+    Creates averaged material. Square refractive index of returned material is arithmetic mean over  regions and  layer_mat. Magnetization (if present) is averaged linearly. 
+
+    """
+    return _libBornAgainCore.createAveragedMaterial(layer_mat, regions)
 class MesoCrystal(IParticle):
     """
 
diff --git a/auto/Wrap/libBornAgainCore_wrap.cpp b/auto/Wrap/libBornAgainCore_wrap.cpp
index 9665a82990e7260e7ae5c2cb6f80748e34d747e6..ca8ee24a4c2dc8630e1a9e65870705fd9a3e58df 100644
--- a/auto/Wrap/libBornAgainCore_wrap.cpp
+++ b/auto/Wrap/libBornAgainCore_wrap.cpp
@@ -5990,7 +5990,7 @@ SWIG_AsVal_std_complex_Sl_double_Sg_  (PyObject *o, std::complex<double>* val)
 
 
 SWIGINTERNINLINE PyObject*
-SWIG_From_std_complex_Sl_double_Sg_  (/*@SWIG:/home/pospelov/software/local/share/swig/3.0.8/typemaps/swigmacros.swg,104,%ifcplusplus@*/
+SWIG_From_std_complex_Sl_double_Sg_  (/*@SWIG:/usr/local/share/swig/3.0.8/typemaps/swigmacros.swg,104,%ifcplusplus@*/
 
 const std::complex<double>&
 
@@ -95795,6 +95795,28 @@ fail:
 }
 
 
+SWIGINTERN PyObject *_wrap_Material_isDefaultMaterial(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+  PyObject *resultobj = 0;
+  Material *arg1 = (Material *) 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  PyObject * obj0 = 0 ;
+  bool result;
+  
+  if (!PyArg_ParseTuple(args,(char *)"O:Material_isDefaultMaterial",&obj0)) SWIG_fail;
+  res1 = SWIG_ConvertPtr(obj0, &argp1,SWIGTYPE_p_Material, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "Material_isDefaultMaterial" "', argument " "1"" of type '" "Material const *""'"); 
+  }
+  arg1 = reinterpret_cast< Material * >(argp1);
+  result = (bool)((Material const *)arg1)->isDefaultMaterial();
+  resultobj = SWIG_From_bool(static_cast< bool >(result));
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
 SWIGINTERN PyObject *_wrap_Material_scalarSubtrSLD(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
   Material *arg1 = (Material *) 0 ;
@@ -96583,6 +96605,43 @@ fail:
 }
 
 
+SWIGINTERN PyObject *_wrap_createAveragedMaterial(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+  PyObject *resultobj = 0;
+  Material *arg1 = 0 ;
+  std::vector< HomogeneousRegion,std::allocator< HomogeneousRegion > > *arg2 = 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  void *argp2 = 0 ;
+  int res2 = 0 ;
+  PyObject * obj0 = 0 ;
+  PyObject * obj1 = 0 ;
+  SwigValueWrapper< Material > result;
+  
+  if (!PyArg_ParseTuple(args,(char *)"OO:createAveragedMaterial",&obj0,&obj1)) SWIG_fail;
+  res1 = SWIG_ConvertPtr(obj0, &argp1, SWIGTYPE_p_Material,  0  | 0);
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "createAveragedMaterial" "', argument " "1"" of type '" "Material const &""'"); 
+  }
+  if (!argp1) {
+    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "createAveragedMaterial" "', argument " "1"" of type '" "Material const &""'"); 
+  }
+  arg1 = reinterpret_cast< Material * >(argp1);
+  res2 = SWIG_ConvertPtr(obj1, &argp2, SWIGTYPE_p_std__vectorT_HomogeneousRegion_std__allocatorT_HomogeneousRegion_t_t,  0  | 0);
+  if (!SWIG_IsOK(res2)) {
+    SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "createAveragedMaterial" "', argument " "2"" of type '" "std::vector< HomogeneousRegion,std::allocator< HomogeneousRegion > > const &""'"); 
+  }
+  if (!argp2) {
+    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "createAveragedMaterial" "', argument " "2"" of type '" "std::vector< HomogeneousRegion,std::allocator< HomogeneousRegion > > const &""'"); 
+  }
+  arg2 = reinterpret_cast< std::vector< HomogeneousRegion,std::allocator< HomogeneousRegion > > * >(argp2);
+  result = createAveragedMaterial((Material const &)*arg1,(std::vector< HomogeneousRegion,std::allocator< HomogeneousRegion > > const &)*arg2);
+  resultobj = SWIG_NewPointerObj((new Material(static_cast< const Material& >(result))), SWIGTYPE_p_Material, SWIG_POINTER_OWN |  0 );
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
 SWIGINTERN PyObject *_wrap_new_MesoCrystal(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
   IClusteredParticles *arg1 = 0 ;
@@ -122237,6 +122296,14 @@ static PyMethodDef SwigMethods[] = {
 		"Returns true if material underlying data is nullptr. \n"
 		"\n"
 		""},
+	 { (char *)"Material_isDefaultMaterial", _wrap_Material_isDefaultMaterial, METH_VARARGS, (char *)"\n"
+		"Material_isDefaultMaterial(Material self) -> bool\n"
+		"\n"
+		"bool Material::isDefaultMaterial() const\n"
+		"\n"
+		"Returns true if material has refractive index of (1.0, 0.0) and zero magnetization. \n"
+		"\n"
+		""},
 	 { (char *)"Material_scalarSubtrSLD", _wrap_Material_scalarSubtrSLD, METH_VARARGS, (char *)"\n"
 		"Material_scalarSubtrSLD(Material self, WavevectorInfo wavevectors) -> complex_t\n"
 		"\n"
@@ -122262,7 +122329,7 @@ static PyMethodDef SwigMethods[] = {
 		"\n"
 		"BA_CORE_API_ Material HomogeneousMaterial()\n"
 		"\n"
-		"Constructs vacuum material based on refractive coefficients. Though in practice there is no difference between vacuum materials produced with  MaterialBySLD() and  HomogeneousMaterial(), they are not equal because of the difference in the type of underlying data \n"
+		"Constructs material with zero refractive coefficients and zero magnetization. \n"
 		"\n"
 		""},
 	 { (char *)"MaterialByAbsCX", _wrap_MaterialByAbsCX, METH_VARARGS, (char *)"\n"
@@ -122293,7 +122360,15 @@ static PyMethodDef SwigMethods[] = {
 		"\n"
 		"BA_CORE_API_ Material MaterialBySLD()\n"
 		"\n"
-		"Constructs wavelength-independent vacuum material. Though in practice there is no difference between vacuum materials produced with  MaterialBySLD() and  HomogeneousMaterial(), they are not equal because of the difference in the type of underlying data \n"
+		"Constructs wavelength-independent material with zero sld and zero magnetization. \n"
+		"\n"
+		""},
+	 { (char *)"createAveragedMaterial", _wrap_createAveragedMaterial, METH_VARARGS, (char *)"\n"
+		"createAveragedMaterial(Material layer_mat, std::vector< HomogeneousRegion,std::allocator< HomogeneousRegion > > const & regions) -> Material\n"
+		"\n"
+		"BA_CORE_API_ Material createAveragedMaterial(const Material &layer_mat, const std::vector< HomogeneousRegion > &regions)\n"
+		"\n"
+		"Creates averaged material. Square refractive index of returned material is arithmetic mean over  regions and  layer_mat. Magnetization (if present) is averaged linearly. \n"
 		"\n"
 		""},
 	 { (char *)"new_MesoCrystal", _wrap_new_MesoCrystal, METH_VARARGS, (char *)"\n"