From f82f5dd8bc7e8f63cc2cae88f033972e655f90d5 Mon Sep 17 00:00:00 2001
From: Walter Van Herck <w.van.herck@fz-juelich.de>
Date: Thu, 12 Sep 2013 17:05:59 +0200
Subject: [PATCH] Fixed calculation of polarized mesocrystals

---
 .cproject                                     |   3 -
 App/App.pro                                   |   2 +
 App/inc/TestPolarizedMeso.h                   |  57 ++++++
 App/src/FunctionalTestFactory.cpp             |   6 +
 App/src/TestPolarizedDWBA.cpp                 |   2 +-
 App/src/TestPolarizedMeso.cpp                 | 176 ++++++++++++++++++
 Core/Algorithms/src/LayerStrategyBuilder.cpp  |  33 ++--
 Core/FormFactors/inc/FormFactorCrystal.h      |   2 +-
 Core/FormFactors/src/FormFactorCrystal.cpp    |  18 +-
 .../src/FormFactorDecoratorMaterial.cpp       |   4 +-
 ...FormFactorDecoratorMultiPositionFactor.cpp |   4 +-
 Core/Samples/inc/Crystal.h                    |   2 +-
 Core/Samples/inc/IClusteredParticles.h        |   5 +
 Core/Samples/inc/ISampleVisitor.h             |  11 ++
 Core/Samples/inc/LatticeBasis.h               |   2 +-
 Core/Samples/inc/MesoCrystal.h                |   5 +-
 Core/Samples/src/LatticeBasis.cpp             |   2 +-
 Core/Samples/src/MesoCrystal.cpp              |   6 +
 Core/Tools/inc/FastVector.h                   |   6 +-
 Core/Tools/inc/SampleMaterialVisitor.h        |   3 +
 Core/Tools/inc/SamplePrintVisitor.h           |   3 +
 Core/Tools/src/SampleMaterialVisitor.cpp      |  30 +++
 Core/Tools/src/SamplePrintVisitor.cpp         |  52 +++++-
 .../SampleDesigner/ISampleToIView.h           |   4 +
 XCode_BornAgain.xcodeproj/project.pbxproj     |   6 +
 25 files changed, 399 insertions(+), 45 deletions(-)
 create mode 100644 App/inc/TestPolarizedMeso.h
 create mode 100644 App/src/TestPolarizedMeso.cpp

diff --git a/.cproject b/.cproject
index b511b7401b1..98b416ee0f8 100644
--- a/.cproject
+++ b/.cproject
@@ -71,9 +71,6 @@
 		</cconfiguration>
 		<cconfiguration id="cdt.managedbuild.toolchain.gnu.macosx.base.45089428.1343516600">
 			<storageModule buildSystemId="org.eclipse.cdt.managedbuilder.core.configurationDataProvider" id="cdt.managedbuild.toolchain.gnu.macosx.base.45089428.1343516600" moduleId="org.eclipse.cdt.core.settings" name="Debug">
-				<macros>
-					<stringMacro name="BORNAGAIN_DEBUG" type="VALUE_TEXT" value="yes"/>
-				</macros>
 				<externalSettings>
 					<externalSetting>
 						<entry flags="RESOLVED" kind="macro" name="BornAgain_DEBUG" value="yes"/>
diff --git a/App/App.pro b/App/App.pro
index d03b87353c1..9e88f607296 100644
--- a/App/App.pro
+++ b/App/App.pro
@@ -83,6 +83,7 @@ SOURCES += \
     src/TestPolarizedDWBA.cpp \
     src/TestPolarizedDWBATerms.cpp \
     src/TestPolarizedDWBAZeroMag.cpp \
+    src/TestPolarizedMeso.cpp \
     src/TestRootTree.cpp \
     src/TestRoughness.cpp \
     src/TestSpecularMagnetic.cpp \
@@ -147,6 +148,7 @@ HEADERS += \
     inc/TestPolarizedDWBA.h \
     inc/TestPolarizedDWBATerms.h \
     inc/TestPolarizedDWBAZeroMag.h \
+    inc/TestPolarizedMeso.h \
     inc/TestRootTree.h \
     inc/TestRoughness.h \
     inc/TestSpecularMagnetic.h \
diff --git a/App/inc/TestPolarizedMeso.h b/App/inc/TestPolarizedMeso.h
new file mode 100644
index 00000000000..8776205b474
--- /dev/null
+++ b/App/inc/TestPolarizedMeso.h
@@ -0,0 +1,57 @@
+// ************************************************************************** //
+//
+//  BornAgain: simulate and fit scattering at grazing incidence
+//
+//! @file      App/inc/TestPolarizedMeso.h
+//! @brief     Defines class TestPolarizedMeso.
+//!
+//! @homepage  http://apps.jcns.fz-juelich.de/BornAgain
+//! @license   GNU General Public License v3 or higher (see COPYING)
+//! @copyright Forschungszentrum Jülich GmbH 2013
+//! @authors   Scientific Computing Group at MLZ Garching
+//! @authors   C. Durniak, G. Pospelov, W. Van Herck, J. Wuttke
+//
+// ************************************************************************** //
+
+#ifndef TESTPOLARIZEDMESO_H_
+#define TESTPOLARIZEDMESO_H_
+
+#include "IFunctionalTest.h"
+#include "MultiLayer.h"
+
+class MesoCrystal;
+class Lattice;
+
+//! Functional test for polarized DWBA with mesocrystal.
+
+class TestPolarizedMeso : public IFunctionalTest
+{
+public:
+    TestPolarizedMeso();
+    virtual ~TestPolarizedMeso();
+
+    virtual void execute();
+
+private:
+    MultiLayer *createSample() const;
+    MesoCrystal* createMeso(double a, double c, const IMaterial *p_material,
+            double cube_size, const IFormFactor* p_meso_form_factor) const;
+    const Lattice *createLattice(double a, double c) const;
+    MultiLayer *mp_sample; //!< pointer to multilayer sample
+    double m_meso_width;
+    double m_surface_filling_ratio;
+    double m_meso_height;
+//    double m_sigma_meso_height;
+//    double m_sigma_meso_radius;
+    double m_lattice_length_a;
+    double m_lattice_length_c;
+    double m_nanoparticle_size;
+//    double m_sigma_nanoparticle_size;
+//    double m_sigma_lattice_length_a;
+//    double m_roughness;
+};
+
+
+
+
+#endif /* TESTPOLARIZEDMESO_H_ */
diff --git a/App/src/FunctionalTestFactory.cpp b/App/src/FunctionalTestFactory.cpp
index c2d75d83876..b9e4f95d82a 100644
--- a/App/src/FunctionalTestFactory.cpp
+++ b/App/src/FunctionalTestFactory.cpp
@@ -49,6 +49,7 @@
 #include "TestPolarizedDWBA.h"
 #include "TestPolarizedDWBATerms.h"
 #include "TestPolarizedDWBAZeroMag.h"
+#include "TestPolarizedMeso.h"
 #include "TestRootTree.h"
 #include "TestRoughness.h"
 #include "TestSpecularMagnetic.h"
@@ -324,5 +325,10 @@ void RegisterFunctionalTests(FunctionalTestFactory *p_test_factory)
          IFactoryCreateFunction<TestPolarizedDWBATerms, IFunctionalTest>,
          "functional test: compare different terms in DWBA between"
          " scalar and matrix computation");
+    p_test_factory->registerItem(
+         "polarizedMeso",
+         IFactoryCreateFunction<TestPolarizedMeso, IFunctionalTest>,
+         "functional test: polarized mesocrystals");
+
 }
 
diff --git a/App/src/TestPolarizedDWBA.cpp b/App/src/TestPolarizedDWBA.cpp
index 6357f331a00..49a992d5e43 100644
--- a/App/src/TestPolarizedDWBA.cpp
+++ b/App/src/TestPolarizedDWBA.cpp
@@ -35,7 +35,7 @@ void TestPolarizedDWBA::execute()
     mp_sample = dynamic_cast<MultiLayer *>(SampleFactory::createSample(
             "PolarizedDWBATestCase"));
 
-    // calculate scattered intesity from sample
+    // calculate scattered intensity from sample
     Simulation simulation(mp_options);
     simulation.setDetectorParameters(
         100, -1.0*Units::degree, 1.0*Units::degree, 100,
diff --git a/App/src/TestPolarizedMeso.cpp b/App/src/TestPolarizedMeso.cpp
new file mode 100644
index 00000000000..b00ea5e3bf3
--- /dev/null
+++ b/App/src/TestPolarizedMeso.cpp
@@ -0,0 +1,176 @@
+// ************************************************************************** //
+//
+//  BornAgain: simulate and fit scattering at grazing incidence
+//
+//! @file      App/src/TestPolarizedMeso.cpp
+//! @brief     Implements class TestPolarizedMeso.
+//!
+//! @homepage  http://apps.jcns.fz-juelich.de/BornAgain
+//! @license   GNU General Public License v3 or higher (see COPYING)
+//! @copyright Forschungszentrum Jülich GmbH 2013
+//! @authors   Scientific Computing Group at MLZ Garching
+//! @authors   C. Durniak, G. Pospelov, W. Van Herck, J. Wuttke
+//
+// ************************************************************************** //
+
+#include "TestPolarizedMeso.h"
+#include "Units.h"
+#include "IsGISAXSTools.h"
+#include "MaterialManager.h"
+#include "InterferenceFunctionNone.h"
+#include "Rotate3D.h"
+#include "FormFactors.h"
+#include "MesoCrystal.h"
+
+
+TestPolarizedMeso::TestPolarizedMeso()
+: mp_sample(0)
+, m_meso_width(3e+3*Units::nanometer)
+, m_surface_filling_ratio(0.3)
+, m_meso_height(3e+3*Units::nanometer)
+//, m_sigma_meso_height(10*Units::nanometer)
+//, m_sigma_meso_radius(10*Units::nanometer)
+, m_lattice_length_a(13.5*Units::nanometer)
+, m_lattice_length_c(23.8*Units::nanometer)
+, m_nanoparticle_size(9.6*Units::nanometer)
+//, m_sigma_nanoparticle_size(6.8688e-02*Units::nanometer)
+//, m_sigma_lattice_length_a(2.3596*Units::nanometer)
+//, m_roughness(0.6517*Units::nanometer)
+{
+    std::cout << "TestPolarizedMeso::TestPolarizedMeso() -> Info."
+            << std::endl;
+}
+
+TestPolarizedMeso::~TestPolarizedMeso()
+{
+}
+
+void TestPolarizedMeso::execute()
+{
+    std::cout << "TestPolarizedMeso::execute() -> Info." << std::endl;
+
+    // create sample
+    mp_sample = createSample();
+
+    // calculate scattered intensity from sample
+    Simulation simulation(mp_options);
+    simulation.setDetectorParameters(
+        100, -1.0*Units::degree, 1.0*Units::degree, 100,
+        0.0*Units::degree, 2.0*Units::degree);
+    simulation.setBeamParameters(
+        1.0*Units::angstrom, 0.2*Units::degree, 0.0*Units::degree);
+    simulation.setBeamIntensity(1e8);
+
+    // Run simulation
+    simulation.setSample(*mp_sample);
+    simulation.runSimulation();
+
+    simulation.normalize();
+
+    IsGISAXSTools::drawLogOutputDataPol(*simulation.getPolarizedOutputData(),
+            "c1_polMeso", "Polarized Mesocrystal", "CONT4 Z",
+            "Polarized Mesocrystal");
+//    IsGISAXSTools::drawLogOutputData(*simulation.getOutputData(), "c1_test_meso_crystal", "mesocrystal",
+//            "CONT4 Z", "mesocrystal");
+
+    delete mp_sample;
+
+}
+
+MultiLayer* TestPolarizedMeso::createSample() const
+{
+    // create mesocrystal
+    double surface_density =
+        m_surface_filling_ratio/m_meso_width/m_meso_width;
+    kvector_t magnetic_field(0.0, 0.0, 6.0);
+    const IMaterial *p_particle_material =
+            MaterialManager::getHomogeneousMagneticMaterial("nanoparticle",
+                    2.84e-5, 4.7e-7, magnetic_field);
+//    const IMaterial *p_particle_material =
+//            MaterialManager::getHomogeneousMaterial("nanoparticle",
+//                    2.84e-5, 4.7e-7);
+    FormFactorBox ff_box(m_meso_width, m_meso_width, m_meso_height);
+//    FormFactorDecoratorDebyeWaller
+//        ff_meso(ff_cyl.clone(),
+//                m_sigma_meso_height*m_sigma_meso_height/2.0,
+//                m_sigma_meso_radius*m_sigma_meso_radius/2.0);
+
+    // Create multilayer
+    MultiLayer *p_multi_layer = new MultiLayer();
+
+    const IMaterial *p_air_material =
+        MaterialManager::getHomogeneousMaterial("Air", 0.0, 0.0);
+    const IMaterial *p_substrate_material =
+        MaterialManager::getHomogeneousMaterial("Substrate", 7.57e-6, 1.73e-7);
+    Layer air_layer;
+    air_layer.setMaterial(p_air_material);
+    Layer substrate_layer;
+    substrate_layer.setMaterial(p_substrate_material);
+    IInterferenceFunction *p_interference_funtion =
+        new InterferenceFunctionNone();
+    ParticleDecoration particle_decoration;
+    size_t n_max_phi_rotation_steps = 1;
+    size_t n_alpha_rotation_steps = 1;
+
+    double phi_step = M_PI/4.0/n_max_phi_rotation_steps;
+    double phi_start = 0.0;
+    Geometry::PTransform3D trafo;
+    for (size_t i=0; i<n_max_phi_rotation_steps; ++i) {
+        for (size_t j=0; j<n_alpha_rotation_steps; ++j) {
+            trafo = Geometry::PTransform3D(new
+                Geometry::RotateZ_3D(phi_start + i*phi_step));
+            particle_decoration.addParticle(
+                createMeso(m_lattice_length_a, m_lattice_length_c,
+                        p_particle_material, m_nanoparticle_size, &ff_box),
+                trafo);
+        }
+    }
+
+    particle_decoration.setTotalParticleSurfaceDensity(surface_density);
+    particle_decoration.addInterferenceFunction(p_interference_funtion);
+
+    air_layer.setDecoration(particle_decoration);
+
+//    LayerRoughness roughness(m_roughness, 0.3, 500.0*Units::nanometer);
+
+    p_multi_layer->addLayer(air_layer);
+    p_multi_layer->addLayer(substrate_layer);
+
+    return p_multi_layer;
+}
+
+MesoCrystal* TestPolarizedMeso::createMeso(double a, double c,
+        const IMaterial *p_material, double cube_size,
+        const IFormFactor* p_meso_form_factor) const
+{
+    const Lattice *p_lat = createLattice(a, c);
+    kvector_t bas_a = p_lat->getBasisVectorA();
+    kvector_t bas_b = p_lat->getBasisVectorB();
+    kvector_t bas_c = p_lat->getBasisVectorC();
+
+
+    Particle particle(p_material, new FormFactorBox(cube_size, cube_size, cube_size));
+    kvector_t position_0 = kvector_t(0.0, 0.0, 0.0);
+    kvector_t position_1 = 1.0/2.0*(bas_a + bas_b + bas_c);
+    std::vector<kvector_t> pos_vector;
+    pos_vector.push_back(position_0);
+    pos_vector.push_back(position_1);
+    LatticeBasis basis(particle, pos_vector);
+
+    Crystal npc(basis, *p_lat);
+    delete p_lat;
+//    double dw_factor = m_sigma_lattice_length_a*m_sigma_lattice_length_a/6.0;
+//    npc.setDWFactor(dw_factor);
+    return new MesoCrystal(npc.clone(), p_meso_form_factor->clone());
+}
+
+const Lattice *TestPolarizedMeso::createLattice(double a, double c) const
+{
+    kvector_t a_vec(a, 0.0, 0.0);
+    kvector_t b_vec(0.0, a, 0.0);
+    kvector_t c_vec(0.0, 0.0, c);
+
+    Lattice *p_result = new Lattice(a_vec, b_vec, c_vec);
+    return p_result;
+}
+
diff --git a/Core/Algorithms/src/LayerStrategyBuilder.cpp b/Core/Algorithms/src/LayerStrategyBuilder.cpp
index 4149c1a1250..f25af31d227 100644
--- a/Core/Algorithms/src/LayerStrategyBuilder.cpp
+++ b/Core/Algorithms/src/LayerStrategyBuilder.cpp
@@ -186,22 +186,6 @@ FormFactorInfo *LayerStrategyBuilder::createFormFactorInfo(
     return p_result;
 }
 
-// =============================================================================
-// Implementation of FormFactorInfo
-// =============================================================================
-
-FormFactorInfo::~FormFactorInfo() { delete mp_ff; }
-
-FormFactorInfo* FormFactorInfo::clone() const
-{
-    FormFactorInfo *p_result = new FormFactorInfo();
-    p_result->m_abundance = m_abundance;
-    p_result->m_pos_x = m_pos_x;
-    p_result->m_pos_y = m_pos_y;
-    p_result->mp_ff = mp_ff->clone();
-    return p_result;
-}
-
 IFormFactor* LayerStrategyBuilder::createDWBAScalarFormFactor(
         IFormFactor* p_form_factor, double depth) const
 {
@@ -224,3 +208,20 @@ IFormFactor* LayerStrategyBuilder::createDWBAMatrixFormFactor(
     p_result->setSpecularInfo(*mp_specular_info);
     return p_result;
 }
+
+// =============================================================================
+// Implementation of FormFactorInfo
+// =============================================================================
+
+FormFactorInfo::~FormFactorInfo() { delete mp_ff; }
+
+FormFactorInfo* FormFactorInfo::clone() const
+{
+    FormFactorInfo *p_result = new FormFactorInfo();
+    p_result->m_abundance = m_abundance;
+    p_result->m_pos_x = m_pos_x;
+    p_result->m_pos_y = m_pos_y;
+    p_result->mp_ff = mp_ff->clone();
+    return p_result;
+}
+
diff --git a/Core/FormFactors/inc/FormFactorCrystal.h b/Core/FormFactors/inc/FormFactorCrystal.h
index ac3ce0dc5f5..b614ec314d3 100644
--- a/Core/FormFactors/inc/FormFactorCrystal.h
+++ b/Core/FormFactors/inc/FormFactorCrystal.h
@@ -53,7 +53,7 @@ private:
 
     Lattice m_lattice;
     complex_t m_wavevector_scattering_factor;
-    Particle *mp_particle;
+    LatticeBasis *mp_lattice_basis;
     IFormFactor *mp_basis_form_factor;
     IFormFactor *mp_meso_form_factor;
     const IMaterial *mp_ambient_material;
diff --git a/Core/FormFactors/src/FormFactorCrystal.cpp b/Core/FormFactors/src/FormFactorCrystal.cpp
index 606e2470053..16863f93586 100644
--- a/Core/FormFactors/src/FormFactorCrystal.cpp
+++ b/Core/FormFactors/src/FormFactorCrystal.cpp
@@ -25,8 +25,8 @@ FormFactorCrystal::FormFactorCrystal(
 , m_max_rec_length(0.0)
 {
     setName("FormFactorCrystal");
-    mp_particle = p_crystal.createBasis();
-    mp_basis_form_factor = mp_particle->createFormFactor(
+    mp_lattice_basis = p_crystal.createBasis();
+    mp_basis_form_factor = mp_lattice_basis->createFormFactor(
             m_wavevector_scattering_factor);
     mp_meso_form_factor = meso_crystal_form_factor.clone();
     setAmbientMaterial(mp_ambient_material);
@@ -35,14 +35,14 @@ FormFactorCrystal::FormFactorCrystal(
 
 FormFactorCrystal::~FormFactorCrystal()
 {
-    delete mp_particle;
+    delete mp_lattice_basis;
     delete mp_basis_form_factor;
     delete mp_meso_form_factor;
 }
 
 FormFactorCrystal* FormFactorCrystal::clone() const
 {
-    Crystal np_crystal(*mp_particle, m_lattice);
+    Crystal np_crystal(*mp_lattice_basis, m_lattice);
     FormFactorCrystal *result = new FormFactorCrystal(np_crystal,
             *mp_meso_form_factor, mp_ambient_material,
             m_wavevector_scattering_factor);
@@ -52,7 +52,7 @@ FormFactorCrystal* FormFactorCrystal::clone() const
 
 void FormFactorCrystal::setAmbientMaterial(const IMaterial *p_material)
 {
-    mp_particle->setAmbientMaterial(p_material);
+    mp_lattice_basis->setAmbientMaterial(p_material);
     mp_basis_form_factor->setAmbientMaterial(p_material);
 }
 
@@ -72,7 +72,7 @@ complex_t FormFactorCrystal::evaluate(const cvector_t& k_i,
     kvector_t q_real(q.x().real(), q.y().real(), q.z().real());
     cvector_t k_zero;
     // calculate the used radius in function of the reciprocal lattice scale
-    double radius = 1.1*m_max_rec_length;
+    double radius = 2.1*m_max_rec_length;
 
     // retrieve nearest reciprocal lattice vectors
     m_lattice.computeReciprocalLatticeVectorsWithinRadius(q_real, radius);
@@ -92,7 +92,7 @@ complex_t FormFactorCrystal::evaluate(const cvector_t& k_i,
         result += basis_factor*meso_factor;
     }
     // the transformed delta train gets a factor of (2pi)^3/V, but the (2pi)^3
-    // is cancelled by the convolution of Fourier transforms :
+    // is canceled by the convolution of Fourier transforms :
     double volume = m_lattice.getVolume();
     return result/volume;
 }
@@ -106,7 +106,7 @@ Eigen::Matrix2cd FormFactorCrystal::evaluatePol(const cvector_t& k_i,
     kvector_t q_real(q.x().real(), q.y().real(), q.z().real());
     cvector_t k_zero;
     // calculate the used radius in function of the reciprocal lattice scale
-    double radius = 1.1*m_max_rec_length;
+    double radius = 2.1*m_max_rec_length;
 
     // retrieve nearest reciprocal lattice vectors
     m_lattice.computeReciprocalLatticeVectorsWithinRadius(q_real, radius);
@@ -126,7 +126,7 @@ Eigen::Matrix2cd FormFactorCrystal::evaluatePol(const cvector_t& k_i,
         result += basis_factor*meso_factor;
     }
     // the transformed delta train gets a factor of (2pi)^3/V, but the (2pi)^3
-    // is cancelled by the convolution of Fourier transforms :
+    // is canceled by the convolution of Fourier transforms :
     double volume = m_lattice.getVolume();
     return result/volume;
 }
diff --git a/Core/FormFactors/src/FormFactorDecoratorMaterial.cpp b/Core/FormFactors/src/FormFactorDecoratorMaterial.cpp
index 90aed454164..1137719a7b6 100644
--- a/Core/FormFactors/src/FormFactorDecoratorMaterial.cpp
+++ b/Core/FormFactors/src/FormFactorDecoratorMaterial.cpp
@@ -66,9 +66,7 @@ Eigen::Matrix2cd FormFactorDecoratorMaterial::evaluatePol(const cvector_t& k_i,
     time_reverse_conj(0,1) = 1.0;
     time_reverse_conj(1,0) = -1.0;
     // the interaction and time reversal taken together:
-    double alpha_f = alpha_f_bin.getMidPoint();
-    double k_mag2 = k_f_bin.getMidPoint().magxy2().real()
-                  / std::cos(alpha_f) / std::cos(alpha_f);
+    double k_mag2 = 4.0 * M_PI * m_wavevector_scattering_factor.real();
     Eigen::Matrix2cd V_eff = m_wavevector_scattering_factor * time_reverse_conj
             * (mp_material->getScatteringMatrix(k_mag2) -
                mp_ambient_material->getScatteringMatrix(k_mag2));
diff --git a/Core/FormFactors/src/FormFactorDecoratorMultiPositionFactor.cpp b/Core/FormFactors/src/FormFactorDecoratorMultiPositionFactor.cpp
index b8f76861dc9..80fca1e5895 100644
--- a/Core/FormFactors/src/FormFactorDecoratorMultiPositionFactor.cpp
+++ b/Core/FormFactors/src/FormFactorDecoratorMultiPositionFactor.cpp
@@ -43,8 +43,8 @@ Eigen::Matrix2cd FormFactorDecoratorMultiPositionFactor::evaluatePol(
         Bin1D phi_f_bin) const
 {
     cvector_t q = k_i - k_f_bin.getMidPoint();
-    return getPositionsFactor(q)*mp_form_factor->
-               evaluatePol(k_i, k_f_bin, alpha_f_bin, phi_f_bin);
+    Eigen::Matrix2cd ff = mp_form_factor->evaluatePol(k_i, k_f_bin, alpha_f_bin, phi_f_bin);
+    return getPositionsFactor(q)*ff;
 }
 
 complex_t FormFactorDecoratorMultiPositionFactor::getPositionsFactor(
diff --git a/Core/Samples/inc/Crystal.h b/Core/Samples/inc/Crystal.h
index 1349557746c..48358d0c7a2 100644
--- a/Core/Samples/inc/Crystal.h
+++ b/Core/Samples/inc/Crystal.h
@@ -46,7 +46,7 @@ public:
         complex_t wavevector_scattering_factor) const;
 
     Lattice getLattice() const { return m_lattice; }
-    Particle *createBasis() const { return mp_lattice_basis->clone(); }
+    LatticeBasis *createBasis() const { return mp_lattice_basis->clone(); }
 
     const LatticeBasis *getLatticeBasis() const { return mp_lattice_basis; }
 
diff --git a/Core/Samples/inc/IClusteredParticles.h b/Core/Samples/inc/IClusteredParticles.h
index cdbf0ca7775..77acaeb884e 100644
--- a/Core/Samples/inc/IClusteredParticles.h
+++ b/Core/Samples/inc/IClusteredParticles.h
@@ -41,6 +41,11 @@ public:
                 "Error! Not implemented exception");
     }
 
+    //! Calls the ISampleVisitor's visit method
+    virtual void accept(ISampleVisitor *p_visitor) const {
+        p_visitor->visit(this);
+    }
+
     virtual void setAmbientMaterial(const IMaterial *p_ambient_material)=0;
 
     //! @brief create a total form factor for the mesocrystal with a specific
diff --git a/Core/Samples/inc/ISampleVisitor.h b/Core/Samples/inc/ISampleVisitor.h
index ee3fc33e736..a7fbece5110 100644
--- a/Core/Samples/inc/ISampleVisitor.h
+++ b/Core/Samples/inc/ISampleVisitor.h
@@ -26,6 +26,10 @@ class ParticleDecoration;
 class ParticleInfo;
 class Particle;
 class ParticleCoreShell;
+class MesoCrystal;
+class Crystal;
+class LatticeBasis;
+class IClusteredParticles;
 class IFormFactor;
 class IInterferenceFunction;
 class FormFactorFullSphere;
@@ -66,6 +70,13 @@ public:
         throw NotImplementedException(
             "ISampleVisitor::visit(ParticleCoreShell *)");
     }
+    virtual void visit(const MesoCrystal *) { throw NotImplementedException(
+            "ISampleVisitor::visit(MesoCrystal *)"); }
+    virtual void visit(const Crystal *) {
+        throw NotImplementedException(
+            "ISampleVisitor::visit(Crystal *)"); }
+    virtual void visit(const LatticeBasis *) { throw NotImplementedException(
+            "ISampleVisitor::visit(LatticeBasis *)"); }
     virtual void visit(const IFormFactor *) { throw NotImplementedException(
             "ISampleVisitor::visit(IFormFactor *)"); }
     virtual void visit(const FormFactorFullSphere *) {
diff --git a/Core/Samples/inc/LatticeBasis.h b/Core/Samples/inc/LatticeBasis.h
index 6d0419f94f6..be98c904f3b 100644
--- a/Core/Samples/inc/LatticeBasis.h
+++ b/Core/Samples/inc/LatticeBasis.h
@@ -24,7 +24,7 @@ class BA_CORE_API_ LatticeBasis : public Particle
 {
 public:
     LatticeBasis();
-    LatticeBasis(const Particle& particle);
+    explicit LatticeBasis(const Particle& particle);
     LatticeBasis(const Particle& particle, std::vector<kvector_t > positions);
     virtual ~LatticeBasis();
     virtual LatticeBasis *clone() const;
diff --git a/Core/Samples/inc/MesoCrystal.h b/Core/Samples/inc/MesoCrystal.h
index 43c2ca56ccf..abda6908d47 100644
--- a/Core/Samples/inc/MesoCrystal.h
+++ b/Core/Samples/inc/MesoCrystal.h
@@ -41,10 +41,7 @@ public:
         p_visitor->visit(this);
     }
 
-    virtual void setAmbientMaterial(const IMaterial *p_material)
-    {
-        mp_particle_structure->setAmbientMaterial(p_material);
-    }
+    virtual void setAmbientMaterial(const IMaterial *p_material);
 
     virtual IFormFactor* createFormFactor(
             complex_t wavevector_scattering_factor) const;
diff --git a/Core/Samples/src/LatticeBasis.cpp b/Core/Samples/src/LatticeBasis.cpp
index 6f1c289edd1..64695bc2df7 100644
--- a/Core/Samples/src/LatticeBasis.cpp
+++ b/Core/Samples/src/LatticeBasis.cpp
@@ -52,7 +52,7 @@ LatticeBasis* LatticeBasis::clone() const
         p_new->addParticle(*m_particles[index], m_positions_vector[index]);
     }
     p_new->setName(getName());
-    p_new->mp_ambient_material = this->mp_ambient_material;
+    p_new->setAmbientMaterial(this->mp_ambient_material);
     return p_new;
 }
 
diff --git a/Core/Samples/src/MesoCrystal.cpp b/Core/Samples/src/MesoCrystal.cpp
index e1e49552981..17db7110b0e 100644
--- a/Core/Samples/src/MesoCrystal.cpp
+++ b/Core/Samples/src/MesoCrystal.cpp
@@ -53,6 +53,12 @@ MesoCrystal* MesoCrystal::cloneInvertB() const
             mp_meso_form_factor->clone());
 }
 
+void MesoCrystal::setAmbientMaterial(const IMaterial* p_material)
+{
+    Particle::setAmbientMaterial(p_material);
+    mp_particle_structure->setAmbientMaterial(p_material);
+}
+
 IFormFactor* MesoCrystal::createFormFactor(
         complex_t wavevector_scattering_factor) const
 {
diff --git a/Core/Tools/inc/FastVector.h b/Core/Tools/inc/FastVector.h
index 7e17a44b99d..52ae7d5c88c 100644
--- a/Core/Tools/inc/FastVector.h
+++ b/Core/Tools/inc/FastVector.h
@@ -35,9 +35,11 @@ public:
 
     inline void push_back(const kvector_t& k) {
         if(m_current_position == m_max_buff_size) {
-            std::cout << "KVectorContainer::push_back() -> Info. Increasing size of the buffer from " << m_max_buff_size;
+//            std::cout << "KVectorContainer::push_back() -> "
+//                         "Info. Increasing size of the buffer from "
+//                      << m_max_buff_size;
             m_max_buff_size *=2;
-            std::cout << " to " << m_max_buff_size << std::endl;
+//            std::cout << " to " << m_max_buff_size << std::endl;
             m_buffer.resize(m_max_buff_size);
         }
         m_buffer[m_current_position][0] = k[0];
diff --git a/Core/Tools/inc/SampleMaterialVisitor.h b/Core/Tools/inc/SampleMaterialVisitor.h
index 2e633033b6f..098dd5dac69 100644
--- a/Core/Tools/inc/SampleMaterialVisitor.h
+++ b/Core/Tools/inc/SampleMaterialVisitor.h
@@ -34,6 +34,9 @@ public:
     void visit(const ParticleInfo *sample);
     void visit(const Particle *sample);
     void visit(const ParticleCoreShell *sample);
+    void visit(const MesoCrystal *sample);
+    void visit(const Crystal *sample);
+    void visit(const LatticeBasis *sample);
     void visit(const IFormFactor *sample);
 
     void visit(const FormFactorFullSphere *sample);
diff --git a/Core/Tools/inc/SamplePrintVisitor.h b/Core/Tools/inc/SamplePrintVisitor.h
index feb8def5262..052e9fa4367 100644
--- a/Core/Tools/inc/SamplePrintVisitor.h
+++ b/Core/Tools/inc/SamplePrintVisitor.h
@@ -18,6 +18,9 @@ public:
     void visit(const ParticleInfo *sample);
     void visit(const Particle *sample);
     void visit(const ParticleCoreShell *sample);
+    void visit(const MesoCrystal *sample);
+    void visit(const Crystal *sample);
+    void visit(const LatticeBasis *sample);
     void visit(const IFormFactor *sample);
 
     void visit(const FormFactorFullSphere *sample);
diff --git a/Core/Tools/src/SampleMaterialVisitor.cpp b/Core/Tools/src/SampleMaterialVisitor.cpp
index bf4de578153..9d294c68953 100644
--- a/Core/Tools/src/SampleMaterialVisitor.cpp
+++ b/Core/Tools/src/SampleMaterialVisitor.cpp
@@ -22,6 +22,7 @@
 #include "Particle.h"
 #include "ParticleInfo.h"
 #include "ParticleCoreShell.h"
+#include "MesoCrystal.h"
 #include "InterferenceFunction1DParaCrystal.h"
 #include "InterferenceFunction2DParaCrystal.h"
 
@@ -113,6 +114,35 @@ void SampleMaterialVisitor::visit(const ParticleCoreShell* sample)
     }
 }
 
+void SampleMaterialVisitor::visit(const MesoCrystal* sample)
+{
+    assert(sample);
+
+    const IClusteredParticles *p_clustered_particles =
+            sample->getClusteredParticles();
+    p_clustered_particles->accept(this);
+}
+
+void SampleMaterialVisitor::visit(const Crystal* sample)
+{
+    assert(sample);
+
+    const LatticeBasis *p_lattice_basis = sample->getLatticeBasis();
+    p_lattice_basis->accept(this);
+}
+
+void SampleMaterialVisitor::visit(const LatticeBasis* sample)
+{
+    assert(sample);
+
+    size_t nbr_particles = sample->getNbrParticles();
+    for (size_t i=0; i<nbr_particles; ++i)
+    {
+        const Particle *p_particle = sample->getParticle(i);
+        p_particle->accept(this);
+    }
+}
+
 void SampleMaterialVisitor::visit(const IFormFactor* sample)
 {
     (void)sample;
diff --git a/Core/Tools/src/SamplePrintVisitor.cpp b/Core/Tools/src/SamplePrintVisitor.cpp
index 3cb74f8ac75..d1b1a415bdf 100644
--- a/Core/Tools/src/SamplePrintVisitor.cpp
+++ b/Core/Tools/src/SamplePrintVisitor.cpp
@@ -9,6 +9,7 @@
 #include "InterferenceFunction2DParaCrystal.h"
 #include "ParticleInfo.h"
 #include <iostream>
+#include "MesoCrystal.h"
 
 
 void SamplePrintVisitor::visit(const ISample *sample)
@@ -153,6 +154,56 @@ void SamplePrintVisitor::visit(const ParticleCoreShell* sample)
     goBack();
 }
 
+void SamplePrintVisitor::visit(const MesoCrystal* sample)
+{
+    assert(sample);
+    std::cout << get_indent() << "PrintVisitor_MesoCrystal " <<
+            sample->getName() << std::endl;
+
+    goForward();
+
+    const IClusteredParticles *p_clustered_particles =
+            sample->getClusteredParticles();
+    p_clustered_particles->accept(this);
+
+    const IFormFactor *p_meso_ff = sample->getSimpleFormFactor();
+    p_meso_ff->accept(this);
+
+    goBack();
+}
+
+void SamplePrintVisitor::visit(const Crystal* sample)
+{
+    assert(sample);
+    std::cout << get_indent() << "PrintVisitor_Crystal " <<
+            sample->getName() << std::endl;
+
+    goForward();
+
+    const LatticeBasis *p_lattice_basis = sample->getLatticeBasis();
+    p_lattice_basis->accept(this);
+
+    goBack();
+}
+
+void SamplePrintVisitor::visit(const LatticeBasis* sample)
+{
+    assert(sample);
+    std::cout << get_indent() << "PrintVisitor_LatticeBasis " <<
+            sample->getName() << std::endl;
+
+    goForward();
+
+    size_t nbr_particles = sample->getNbrParticles();
+    for (size_t i=0; i<nbr_particles; ++i)
+    {
+        const Particle *p_particle = sample->getParticle(i);
+        p_particle->accept(this);
+    }
+
+    goBack();
+}
+
 void SamplePrintVisitor::visit(const IFormFactor *sample)
 {
     assert(sample);
@@ -216,7 +267,6 @@ void SamplePrintVisitor::visit(const InterferenceFunction1DParaCrystal *sample)
               << std::endl;
 }
 
-
 void SamplePrintVisitor::visit(const InterferenceFunction2DParaCrystal *sample)
 {
     assert(sample);
diff --git a/GUI/coregui/Views/Components/SampleDesigner/ISampleToIView.h b/GUI/coregui/Views/Components/SampleDesigner/ISampleToIView.h
index f98c451b7d1..0117b170646 100644
--- a/GUI/coregui/Views/Components/SampleDesigner/ISampleToIView.h
+++ b/GUI/coregui/Views/Components/SampleDesigner/ISampleToIView.h
@@ -23,6 +23,10 @@ public:
     void visit(const ParticleDecoration *sample);
     void visit(const ParticleInfo *sample);
     void visit(const Particle *sample);
+    void visit(const ParticleCoreShell *sample);
+    void visit(const MesoCrystal *sample);
+    void visit(const Crystal *sample);
+    void visit(const LatticeBasis *sample);
     void visit(const IFormFactor *sample);
 
     void visit(const FormFactorFullSphere *sample);
diff --git a/XCode_BornAgain.xcodeproj/project.pbxproj b/XCode_BornAgain.xcodeproj/project.pbxproj
index 17c05760e25..c36089ebef9 100644
--- a/XCode_BornAgain.xcodeproj/project.pbxproj
+++ b/XCode_BornAgain.xcodeproj/project.pbxproj
@@ -125,6 +125,7 @@
 		622F80DB17DE1ED20017FC52 /* FormFactorDecoratorMaterial.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 622F80DA17DE1ED20017FC52 /* FormFactorDecoratorMaterial.cpp */; };
 		6236DD0916CE708600ECED4F /* Instrument.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 6236DD0816CE708600ECED4F /* Instrument.cpp */; };
 		6236DD0D16CE9EC600ECED4F /* Simulation.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 6236DD0C16CE9EC600ECED4F /* Simulation.cpp */; };
+		6248BE2917E1C5000027D960 /* TestPolarizedMeso.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 6248BE2817E1C5000027D960 /* TestPolarizedMeso.cpp */; };
 		6254C2651666652E0098EE7E /* IFormFactorBorn.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 6254C2641666652E0098EE7E /* IFormFactorBorn.cpp */; };
 		625A174116BAAE77004943DB /* FormFactorCone.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 625A173E16BAAE77004943DB /* FormFactorCone.cpp */; };
 		625A174216BAAE77004943DB /* FormFactorFullSpheroid.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 625A173F16BAAE77004943DB /* FormFactorFullSpheroid.cpp */; };
@@ -1453,6 +1454,8 @@
 		6236DD0A16CE709400ECED4F /* MemberComplexFunctionIntegrator.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = MemberComplexFunctionIntegrator.h; sourceTree = "<group>"; };
 		6236DD0B16CE9EBD00ECED4F /* Simulation.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = Simulation.h; sourceTree = "<group>"; };
 		6236DD0C16CE9EC600ECED4F /* Simulation.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = Simulation.cpp; sourceTree = "<group>"; };
+		6248BE1A17E1C4F40027D960 /* TestPolarizedMeso.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = TestPolarizedMeso.h; sourceTree = "<group>"; };
+		6248BE2817E1C5000027D960 /* TestPolarizedMeso.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = TestPolarizedMeso.cpp; sourceTree = "<group>"; };
 		6254C2601666651F0098EE7E /* IFormFactorBorn.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = IFormFactorBorn.h; sourceTree = "<group>"; };
 		6254C2611666651F0098EE7E /* IFormFactorBornSeparable.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = IFormFactorBornSeparable.h; sourceTree = "<group>"; };
 		6254C2621666651F0098EE7E /* IFormFactorDecorator.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = IFormFactorDecorator.h; sourceTree = "<group>"; };
@@ -2471,6 +2474,7 @@
 		622221E9160CB745008205AC /* inc */ = {
 			isa = PBXGroup;
 			children = (
+				6248BE1A17E1C4F40027D960 /* TestPolarizedMeso.h */,
 				622127AB17D9FEB400393360 /* TestPolarizedDWBATerms.h */,
 				62ED4B4817D8B20700C8BAA0 /* TestPolarizedDWBAZeroMag.h */,
 				6221C37017D5CFDF00D43C2F /* TestBugs.h */,
@@ -2532,6 +2536,7 @@
 		62222224160CB745008205AC /* src */ = {
 			isa = PBXGroup;
 			children = (
+				6248BE2817E1C5000027D960 /* TestPolarizedMeso.cpp */,
 				622127B917D9FEC200393360 /* TestPolarizedDWBATerms.cpp */,
 				62ED4B4917D8B21200C8BAA0 /* TestPolarizedDWBAZeroMag.cpp */,
 				6221C37317D5CFF500D43C2F /* TestBugs.cpp */,
@@ -71097,6 +71102,7 @@
 				6221C37817D5CFF500D43C2F /* TestSpecularMagnetic.cpp in Sources */,
 				62ED4B4A17D8B21200C8BAA0 /* TestPolarizedDWBAZeroMag.cpp in Sources */,
 				622127BA17D9FEC200393360 /* TestPolarizedDWBATerms.cpp in Sources */,
+				6248BE2917E1C5000027D960 /* TestPolarizedMeso.cpp in Sources */,
 			);
 			runOnlyForDeploymentPostprocessing = 0;
 		};
-- 
GitLab