From d324c2a56f6c47d7496802e9ec586a980c149d67 Mon Sep 17 00:00:00 2001
From: Gennady Pospelov <g.pospelov@fz-juelich.de>
Date: Tue, 13 Jan 2015 17:24:54 +0100
Subject: [PATCH] New unit test for MultiLayer, refactoring in
 SpecularSimulation, corresponding unit tests.

---
 Core/Algorithms/inc/SpecularSimulation.h      |  34 ++-
 Core/Algorithms/src/SpecularSimulation.cpp    | 222 ++++++++++--------
 Core/PythonAPI/src/PythonModule.cpp           | 212 ++++++++---------
 .../PythonAPI/src/SpecularSimulation.pypp.cpp |  44 +++-
 Core/Samples/inc/MultiLayer.h                 |   3 +
 Core/Samples/src/MultiLayer.cpp               |   8 +
 .../ex05_BeamAndDetector/EvanescentWaveMap.py | 100 ++++++++
 .../SpecularSimulation.py                     |   2 -
 Fit/PythonAPI/src/PythonModule.cpp            |  46 ++--
 Tests/UnitTests/TestCore/MultiLayerTest.h     |  46 ++++
 .../TestCore/SpecularSimulationTest.h         | 133 +++++++++++
 Tests/UnitTests/TestCore/main.cpp             |   1 +
 dev-tools/python-bindings/settings_core.py    |   5 +
 13 files changed, 607 insertions(+), 249 deletions(-)
 create mode 100644 Examples/python/simulation/ex05_BeamAndDetector/EvanescentWaveMap.py
 create mode 100644 Tests/UnitTests/TestCore/SpecularSimulationTest.h

diff --git a/Core/Algorithms/inc/SpecularSimulation.h b/Core/Algorithms/inc/SpecularSimulation.h
index 2ad1a897049..e161d736b62 100644
--- a/Core/Algorithms/inc/SpecularSimulation.h
+++ b/Core/Algorithms/inc/SpecularSimulation.h
@@ -16,7 +16,6 @@
 #ifndef SPECULARSIMULATION_H
 #define SPECULARSIMULATION_H
 
-#include "Simulation.h"
 //#include "MultiLayerRTCoefficients.h"
 #include "OutputData.h"
 
@@ -39,13 +38,10 @@ public:
     SpecularSimulation();
     SpecularSimulation(const ISample& sample);
     SpecularSimulation(SampleBuilder_t sample_builder);
-    ~SpecularSimulation();
+    virtual ~SpecularSimulation();
 
     SpecularSimulation *clone() const;
 
-    //! Put into a clean state for running a simulation
-    void prepareSimulation();
-
     //! Run a simulation with the current parameter settings
     void runSimulation();
 
@@ -65,19 +61,28 @@ public:
     void setBeamParameters(double lambda, const IAxis &alpha_axis);
     void setBeamParameters(double lambda, int nbins, double alpha_i_min, double alpha_i_max);
 
+    //! set axis for evanescent wave axis
+    void setEvanescentWaveAxis(const IAxis &z_axis);
+    void setEvanescentWaveAxis(int nbins, double z_min, double z_max);
+
     //! returns alpha_i axis
     const IAxis *getAlphaAxis() const;
 
     //! returns vector containing reflection coefficients for all alpha_i angles for given layer index
-    std::vector<complex_t > getScalarR(int i_layer = 0) const;
+    std::vector<complex_t > getScalarR(size_t i_layer) const;
 
     //! returns vector containing transmission coefficients for all alpha_i angles for given layer index
-    std::vector<complex_t > getScalarT(int i_layer = 0) const;
+    std::vector<complex_t > getScalarT(size_t i_layer) const;
 
     //! returns vector containing Kz coefficients for all alpha_i angles for given layer index
-    std::vector<complex_t > getScalarKz(int i_layer = 0) const;
+    std::vector<complex_t > getScalarKz(size_t i_layer) const;
 
-    LayerRTCoefficients_t getLayerRTCoefficients(int i_alpha, int i_layer) const;
+    LayerRTCoefficients_t getLayerRTCoefficients(size_t i_alpha, size_t i_layer) const;
+
+    //! Put into a clean state for running a simulation
+    void prepareSimulation();
+
+    OutputData<double>* getEvanescentWaveIntensity() const;
 
 protected:
     SpecularSimulation(const SpecularSimulation& other);
@@ -94,15 +99,26 @@ protected:
     //! calculates RT coefficients for multilayer with magnetic materials
     void collectRTCoefficientsMatrix(const MultiLayer *multilayer);
 
+    //! calculates the intensity of evanescent wave
+    void calculateEvanescentWaveIntensity();
+
+    //! check if simulation was run already and has valid coefficients
+    void checkCoefficients(size_t i_layer) const;
+
+    //! update data axes
+    void updateCoefficientDataAxes();
+
     ISample *m_sample;
     SampleBuilder_t m_sample_builder;
     IAxis *m_alpha_i_axis;
+    IAxis *m_z_axis;
     double m_lambda;
 
 #ifndef GCCXML_SKIP_THIS
 //    OutputData<SpecularMatrix::MultiLayerCoeff_t> *m_scalar_data;
 
     OutputData<MultiLayerRTCoefficients_t> m_data;
+    OutputData<double> m_ewave_intensity;
 #endif
 };
 
diff --git a/Core/Algorithms/src/SpecularSimulation.cpp b/Core/Algorithms/src/SpecularSimulation.cpp
index 03f5f8e9ed1..65858123711 100644
--- a/Core/Algorithms/src/SpecularSimulation.cpp
+++ b/Core/Algorithms/src/SpecularSimulation.cpp
@@ -16,14 +16,15 @@
 #include "SpecularSimulation.h"
 #include "MultiLayer.h"
 #include "SpecularMatrix.h"
+#include <iostream>
 
 
 SpecularSimulation::SpecularSimulation()
     : IParameterized("SpecularSimulation")
     , m_sample(0)
     , m_alpha_i_axis(0)
+    , m_z_axis(0)
     , m_lambda(0.0)
-//    , m_scalar_data(0)
 {
     init_parameters();
 }
@@ -33,8 +34,8 @@ SpecularSimulation::SpecularSimulation(const ISample& sample)
     : IParameterized("SpecularSimulation")
     , m_sample(sample.clone())
     , m_alpha_i_axis(0)
+    , m_z_axis(0)
     , m_lambda(0.0)
-//    , m_scalar_data(0)
 {
     init_parameters();
 }
@@ -45,24 +46,25 @@ SpecularSimulation::SpecularSimulation(SampleBuilder_t sample_builder)
     , m_sample(0)
     , m_sample_builder(sample_builder)
     , m_alpha_i_axis(0)
+    , m_z_axis(0)
     , m_lambda(0.0)
-//    , m_scalar_data(0)
 {
     init_parameters();
 }
 
-
 SpecularSimulation::SpecularSimulation(const SpecularSimulation& other)
     : ICloneable(), IParameterized(other)
     , m_sample(0)
     , m_sample_builder(other.m_sample_builder)
+    , m_alpha_i_axis(0)
+    , m_z_axis(0)
     , m_lambda(other.m_lambda)
 {
     if(other.m_sample) m_sample = other.m_sample->clone();
     if(other.m_alpha_i_axis) m_alpha_i_axis = other.m_alpha_i_axis->clone();
-//    if(other.m_scalar_data) m_scalar_data = other.m_scalar_data->clone();
+    if(other.m_z_axis) m_z_axis = other.m_z_axis->clone();
     m_data.copyFrom(other.m_data);
-
+    m_ewave_intensity.copyFrom(other.m_ewave_intensity);
     init_parameters();
 }
 
@@ -71,7 +73,7 @@ SpecularSimulation::~SpecularSimulation()
 {
     delete m_sample;
     delete m_alpha_i_axis;
-//    delete m_scalar_data;
+    delete m_z_axis;
 }
 
 
@@ -115,17 +117,32 @@ SampleBuilder_t SpecularSimulation::getSampleBuilder() const
 
 void SpecularSimulation::prepareSimulation()
 {
+    updateSample();
+
     if (!m_alpha_i_axis || m_alpha_i_axis->getSize()<1) {
         throw ClassInitializationException(
-                "SpecularSimulation::prepareSimulation() "
+                "SpecularSimulation::checkSimulation() "
                 "-> Error. Incoming alpha range not configured.");
     }
     if (m_lambda<=0.0) {
         throw ClassInitializationException(
-                "SpecularSimulation::prepareSimulation() "
-                "-> Error. Incoming wavelength < 0.");
+                "SpecularSimulation::checkSimulation() "
+                "-> Error. Incoming wavelength <= 0.");
     }
-    updateSample();
+
+    if(!m_sample)
+        throw ClassInitializationException("SpecularSimulation::checkSimulation() -> Error. No sample set");
+
+    updateCoefficientDataAxes();
+}
+
+
+OutputData<double> *SpecularSimulation::getEvanescentWaveIntensity() const
+{
+    if(m_ewave_intensity.getAllocatedSize() == 1)
+        throw ClassInitializationException("SpecularSimulation::getEvanescentWaveIntensity() -> Error. No evanescent wave calculations have been performed. Set corresponding axis with Simulation::setEvanescentWaveAxis.");
+
+    return m_ewave_intensity.clone();
 }
 
 
@@ -135,12 +152,13 @@ void SpecularSimulation::runSimulation()
 
     MultiLayer *multilayer = dynamic_cast<MultiLayer *>(m_sample);
     if(!multilayer)
-        throw RuntimeErrorException("SpecularSimulation::runSimulation() -> Error. Not a multilayer");
+        throw NullPointerException("SpecularSimulation::runSimulation() -> Error. Not a MultiLayer");
 
     if(multilayer->requiresMatrixRTCoefficients()) {
         collectRTCoefficientsMatrix(multilayer);
     } else {
         collectRTCoefficientsScalar(multilayer);
+        calculateEvanescentWaveIntensity();
     }
 }
 
@@ -149,75 +167,40 @@ void SpecularSimulation::setBeamParameters(double lambda, const IAxis &alpha_axi
 {
     delete m_alpha_i_axis;
     m_alpha_i_axis = alpha_axis.clone();
-    if (alpha_axis.getSize()<1) {
-        throw ClassInitializationException(
-                "SpecularSimulation::prepareSimulation() "
-                "-> Error. Incoming alpha range size < 1.");
-    }
     m_lambda = lambda;
 }
 
 
 void SpecularSimulation::setBeamParameters(double lambda, int nbins, double alpha_i_min, double alpha_i_max)
 {
-    delete m_alpha_i_axis;
-    m_alpha_i_axis = new FixedBinAxis("alpha_i", nbins, alpha_i_min, alpha_i_max);
-    m_lambda = lambda;
+    FixedBinAxis axis("alpha_i", nbins, alpha_i_min, alpha_i_max);
+    setBeamParameters(lambda, axis);
 }
 
-
-const IAxis *SpecularSimulation::getAlphaAxis() const
+void SpecularSimulation::setEvanescentWaveAxis(const IAxis &z_axis)
 {
-    return m_alpha_i_axis;
+    delete m_z_axis;
+    m_z_axis = z_axis.clone();
 }
 
+void SpecularSimulation::setEvanescentWaveAxis(int nbins, double z_min, double z_max)
+{
+    FixedBinAxis axis("z_axis", nbins, z_min, z_max);
+    setEvanescentWaveAxis(axis);
+}
 
-//std::vector<complex_t > SpecularSimulation::getScalarR(int i_layer) const
-//{
-//    if(!m_scalar_data)
-//        throw RuntimeErrorException("SpecularSimulation::getScalarR() -> Error. No scalar coefficients.");
-
-//    std::vector<complex_t > result;
-//    result.resize(m_alpha_i_axis->getSize());
-//    for(size_t i=0; i<m_scalar_data->getAllocatedSize(); ++i) {
-//        result[i] = (*m_scalar_data)[i][i_layer].getScalarR();
-//    }
-//    return result;
-//}
-
-
-//std::vector<complex_t > SpecularSimulation::getScalarT(int i_layer) const
-//{
-//    if(!m_scalar_data)
-//        throw RuntimeErrorException("SpecularSimulation::getScalarR() -> Error. No scalar coefficients.");
-
-//    std::vector<complex_t > result;
-//    result.resize(m_alpha_i_axis->getSize());
-//    for(size_t i=0; i<m_scalar_data->getAllocatedSize(); ++i) {
-//        result[i] = (*m_scalar_data)[i][i_layer].getScalarT();
-//    }
-//    return result;
-//}
-
-
-//std::vector<complex_t> SpecularSimulation::getScalarKz(int i_layer) const
-//{
-//    if(!m_scalar_data)
-//        throw RuntimeErrorException("SpecularSimulation::getScalarR() -> Error. No scalar coefficients.");
 
-//    std::vector<complex_t > result;
-//    result.resize(m_alpha_i_axis->getSize());
-//    for(size_t i=0; i<m_scalar_data->getAllocatedSize(); ++i) {
-//        result[i] = (*m_scalar_data)[i][i_layer].getScalarKz();
-//    }
-//    return result;
-//}
+const IAxis *SpecularSimulation::getAlphaAxis() const
+{
+    return m_alpha_i_axis;
+}
 
 
-std::vector<complex_t > SpecularSimulation::getScalarR(int i_layer) const
+std::vector<complex_t > SpecularSimulation::getScalarR(size_t i_layer) const
 {
+    checkCoefficients(i_layer);
     std::vector<complex_t > result;
-    result.resize(m_alpha_i_axis->getSize());
+    result.resize(m_data.getAllocatedSize());
     for(size_t i=0; i<m_data.getAllocatedSize(); ++i) {
         result[i] = m_data[i][i_layer]->getScalarR();
     }
@@ -225,10 +208,11 @@ std::vector<complex_t > SpecularSimulation::getScalarR(int i_layer) const
 }
 
 
-std::vector<complex_t > SpecularSimulation::getScalarT(int i_layer) const
+std::vector<complex_t > SpecularSimulation::getScalarT(size_t i_layer) const
 {
+    checkCoefficients(i_layer);
     std::vector<complex_t > result;
-    result.resize(m_alpha_i_axis->getSize());
+    result.resize(m_data.getAllocatedSize());
     for(size_t i=0; i<m_data.getAllocatedSize(); ++i) {
         result[i] = m_data[i][i_layer]->getScalarT();
     }
@@ -236,10 +220,11 @@ std::vector<complex_t > SpecularSimulation::getScalarT(int i_layer) const
 }
 
 
-std::vector<complex_t> SpecularSimulation::getScalarKz(int i_layer) const
+std::vector<complex_t> SpecularSimulation::getScalarKz(size_t i_layer) const
 {
+    checkCoefficients(i_layer);
     std::vector<complex_t > result;
-    result.resize(m_alpha_i_axis->getSize());
+    result.resize(m_data.getAllocatedSize());
     for(size_t i=0; i<m_data.getAllocatedSize(); ++i) {
         result[i] = m_data[i][i_layer]->getScalarKz();
     }
@@ -247,13 +232,13 @@ std::vector<complex_t> SpecularSimulation::getScalarKz(int i_layer) const
 }
 
 
-SpecularSimulation::LayerRTCoefficients_t SpecularSimulation::getLayerRTCoefficients(int i_alpha, int i_layer) const
+SpecularSimulation::LayerRTCoefficients_t SpecularSimulation::getLayerRTCoefficients(size_t i_alpha, size_t i_layer) const
 {
-    if((size_t)i_alpha >= m_data.getAllocatedSize())
+    if(i_alpha >= m_data.getAllocatedSize())
         throw RuntimeErrorException("SpecularSimulation::getLayerRTCoefficients() -> Error. Wrong i_alpha.");
 
-    if((size_t)i_layer >= m_data[i_alpha].size())
-        throw RuntimeErrorException("SpecularSimulation::getLayerRTCoefficients() -> Error. Wrong i_layer.");
+    if(i_layer >= m_data[i_alpha].size())
+        throw RuntimeErrorException("SpecularSimulation::getLayerRTCoefficients() -> Error. Wrong layer index.");
 
     return m_data[i_alpha][i_layer];
 }
@@ -279,41 +264,11 @@ void SpecularSimulation::updateSample()
             m_sample = new_sample;
         }
     }
-
-    if(!m_sample)
-        throw RuntimeErrorException("SpecularSimulation::updateSample() -> Error. Sample is not set.");
 }
 
 
-//void SpecularSimulation::collectRTCoefficientsScalar(const MultiLayer *multilayer)
-//{
-//    delete m_scalar_data;
-//    m_scalar_data = new OutputData<SpecularMatrix::MultiLayerCoeff_t>();
-
-//    m_scalar_data->addAxis(*m_alpha_i_axis);
-
-//    OutputData<SpecularMatrix::MultiLayerCoeff_t>::iterator it = m_scalar_data->begin();
-//    while (it != m_scalar_data->end()) {
-//        double alpha_i = m_scalar_data->getValueOfAxis(0,it.getIndex());
-//        kvector_t kvec;
-//        kvec.setLambdaAlphaPhi(m_lambda, -alpha_i, 0.0);
-
-//        SpecularMatrix::MultiLayerCoeff_t coeffs;
-//        SpecularMatrix matrixCalculator;
-//        matrixCalculator.execute(*multilayer, kvec, coeffs);
-
-//        *it = coeffs;
-//        ++it;
-
-//    } // alpha_i
-//}
-
-
 void SpecularSimulation::collectRTCoefficientsScalar(const MultiLayer *multilayer)
 {
-    m_data.clear();
-    m_data.addAxis(*m_alpha_i_axis);
-
     OutputData<MultiLayerRTCoefficients_t>::iterator it = m_data.begin();
     while (it != m_data.end()) {
         double alpha_i = m_data.getValueOfAxis(0,it.getIndex());
@@ -341,3 +296,66 @@ void SpecularSimulation::collectRTCoefficientsMatrix(const MultiLayer * /*multil
 {
     throw NotImplementedException("SpecularSimulation::collectRTCoefficientsMatrix() -> Error. Not implemented.");
 }
+
+
+void SpecularSimulation::calculateEvanescentWaveIntensity()
+{
+    if(!m_z_axis) return;
+
+    MultiLayer *multilayer = dynamic_cast<MultiLayer *>(m_sample);
+    std::vector<double> z_coordinates;
+    for(size_t i=0; i<multilayer->getNumberOfLayers(); ++i) {
+        z_coordinates.push_back(multilayer->getLayerBottomZ(i));
+    }
+    for(size_t i=0; i<z_coordinates.size(); ++i) {
+        double z = z_coordinates[i] -1.0;
+        std::cout << i << " " << z_coordinates[i] << " z:" << z << " " << multilayer->zToLayerIndex(z) << std::endl;
+    }
+
+    const IAxis *alpha_axis = m_ewave_intensity.getAxis(0);
+    const IAxis *z_axis = m_ewave_intensity.getAxis(1);
+
+    OutputData<double>::iterator it = m_ewave_intensity.begin();
+    while (it != m_ewave_intensity.end()) {
+        std::vector<int> indices =
+                m_ewave_intensity.toCoordinates(it.getIndex());
+
+        double alpha_axis_value = (*alpha_axis)[indices[0]];
+        double z_axis_value = (*z_axis)[indices[1]];
+
+
+        ++it;
+    }
+
+}
+
+
+void SpecularSimulation::checkCoefficients(size_t i_layer) const
+{
+    if(m_data.getAllocatedSize() == 1 || m_data[0].size()==0)
+        throw ClassInitializationException(
+                "SpecularSimulation::checkCoefficients() "
+                "-> Error. No coefficients found, check that (1) you have set beam parameters (2) you have run your simulation.");
+
+    if(i_layer >= m_data[0].size()) {
+        std::ostringstream message;
+        message << "SpecularSimulation::checkCoefficients() -> Error. Requested layer index "
+                << i_layer
+                << " is large than or equal to the total number of layers "
+                << m_data[0].size()
+                << std::endl;
+        throw OutOfBoundsException(message.str());
+    }
+}
+
+void SpecularSimulation::updateCoefficientDataAxes()
+{
+    m_data.clear();
+    m_data.addAxis(*m_alpha_i_axis);
+
+    if(m_z_axis) {
+        m_ewave_intensity.clear();
+        m_ewave_intensity.addAxis(*m_alpha_i_axis);
+        m_ewave_intensity.addAxis(*m_z_axis);
+    }
+}
diff --git a/Core/PythonAPI/src/PythonModule.cpp b/Core/PythonAPI/src/PythonModule.cpp
index 43c2f51e748..58ec0651a51 100644
--- a/Core/PythonAPI/src/PythonModule.cpp
+++ b/Core/PythonAPI/src/PythonModule.cpp
@@ -10,129 +10,129 @@ GCC_DIAG_ON(missing-field-initializers)
 #define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
 #define PY_ARRAY_UNIQUE_SYMBOL BORNAGAIN_PYTHONAPI_ARRAY
 #include "numpy/arrayobject.h"
-#include "FormFactorFullSpheroid.pypp.h"
-#include "DistributionGate.pypp.h"
-#include "SimpleSelectionRule.pypp.h"
-#include "RealParameterWrapper.pypp.h"
-#include "vdouble1d_t.pypp.h"
+#include "FormFactorFullSphere.pypp.h"
+#include "Lattice.pypp.h"
+#include "FormFactorInfLongBox.pypp.h"
+#include "Layer.pypp.h"
+#include "FormFactorHemiEllipsoid.pypp.h"
+#include "FormFactorPrism3.pypp.h"
 #include "SimulationParameters.pypp.h"
-#include "Transform3D.pypp.h"
-#include "ThreadInfo.pypp.h"
-#include "InterferenceFunction2DLattice.pypp.h"
-#include "LayerInterface.pypp.h"
-#include "ILayout.pypp.h"
-#include "FormFactorCone6.pypp.h"
-#include "FormFactorTetrahedron.pypp.h"
-#include "FTDistribution1DCosine.pypp.h"
-#include "FTDistribution1DTriangle.pypp.h"
-#include "FormFactorWeighted.pypp.h"
-#include "DistributionGaussian.pypp.h"
-#include "IDetectorResolution.pypp.h"
-#include "FormFactorCylinder.pypp.h"
-#include "Crystal.pypp.h"
-#include "FTDistribution1DCauchy.pypp.h"
-#include "IFormFactorBorn.pypp.h"
-#include "FormFactorEllipsoidalCylinder.pypp.h"
-#include "InterferenceFunctionNone.pypp.h"
-#include "FTDistribution2DGate.pypp.h"
+#include "FormFactorInfLongRipple1.pypp.h"
+#include "FTDistribution1DVoigt.pypp.h"
+#include "FormFactorRipple1.pypp.h"
+#include "ParticleInfo.pypp.h"
+#include "ICompositeSample.pypp.h"
+#include "OffSpecSimulation.pypp.h"
+#include "IResolutionFunction2D.pypp.h"
 #include "vector_kvector_t.pypp.h"
-#include "FormFactorTruncatedSpheroid.pypp.h"
-#include "Particle.pypp.h"
-#include "FormFactorTrivial.pypp.h"
-#include "ConstKBinAxis.pypp.h"
-#include "FTDistribution2DCauchy.pypp.h"
-#include "FormFactorCrystal.pypp.h"
-#include "ParticleDistribution.pypp.h"
-#include "vector_longinteger_t.pypp.h"
-#include "ResolutionFunction2DGaussian.pypp.h"
-#include "FTDistribution1DGauss.pypp.h"
 #include "FTDistribution1DGate.pypp.h"
-#include "FormFactorAnisoPyramid.pypp.h"
-#include "FixedBinAxis.pypp.h"
-#include "MultiLayer.pypp.h"
-#include "IFormFactor.pypp.h"
-#include "kvector_t.pypp.h"
-#include "FormFactorSphereUniformRadius.pypp.h"
-#include "OffSpecSimulation.pypp.h"
-#include "FormFactorRipple1.pypp.h"
+#include "LatticeBasis.pypp.h"
+#include "FormFactorRipple2.pypp.h"
+#include "FTDistribution2DVoigt.pypp.h"
+#include "IFormFactorBorn.pypp.h"
+#include "FormFactorPyramid.pypp.h"
+#include "ISample.pypp.h"
+#include "ILayout.pypp.h"
+#include "FormFactorTetrahedron.pypp.h"
 #include "InterferenceFunction1DParaCrystal.pypp.h"
 #include "Simulation.pypp.h"
+#include "ParticleLayout.pypp.h"
+#include "FTDistribution2DCone.pypp.h"
+#include "FormFactorEllipsoidalCylinder.pypp.h"
+#include "vector_IFormFactorPtr_t.pypp.h"
+#include "FTDistribution1DGauss.pypp.h"
+#include "HomogeneousMaterial.pypp.h"
+#include "cvector_t.pypp.h"
+#include "FTDistribution2DGauss.pypp.h"
 #include "IObservable.pypp.h"
+#include "FormFactorBox.pypp.h"
+#include "DistributionLogNormal.pypp.h"
+#include "FormFactorCylinder.pypp.h"
+#include "Detector.pypp.h"
+#include "IObserver.pypp.h"
+#include "ResolutionFunction2DGaussian.pypp.h"
+#include "IParameterized.pypp.h"
+#include "FormFactorGauss.pypp.h"
+#include "FormFactorDecoratorDebyeWaller.pypp.h"
+#include "IDetectorResolution.pypp.h"
+#include "FormFactorWeighted.pypp.h"
+#include "InterferenceFunction2DLattice.pypp.h"
+#include "RealParameterWrapper.pypp.h"
+#include "LayerInterface.pypp.h"
+#include "ParticleCoreShell.pypp.h"
+#include "InterferenceFunction1DLattice.pypp.h"
+#include "FormFactorCrystal.pypp.h"
+#include "Instrument.pypp.h"
+#include "PythonInterface_global_variables.pypp.h"
+#include "ParticleDistribution.pypp.h"
+#include "IFTDistribution1D.pypp.h"
+#include "FormFactorTrivial.pypp.h"
+#include "DistributionGaussian.pypp.h"
+#include "Lattice1DIFParameters.pypp.h"
+#include "DistributionCosine.pypp.h"
+#include "IClusteredParticles.pypp.h"
 #include "FormFactorLorentz.pypp.h"
-#include "SpecularSimulation.pypp.h"
-#include "ISelectionRule.pypp.h"
-#include "FormFactorRipple2.pypp.h"
-#include "LayerRoughness.pypp.h"
 #include "Bin1DCVector.pypp.h"
-#include "FormFactorSphereGaussianRadius.pypp.h"
-#include "ParameterPool.pypp.h"
-#include "FormFactorPrism3.pypp.h"
+#include "InterferenceFunctionNone.pypp.h"
+#include "IFormFactorDecorator.pypp.h"
+#include "Bin1D.pypp.h"
+#include "ISampleBuilder.pypp.h"
+#include "IntensityDataIOFactory.pypp.h"
+#include "FTDistribution2DGate.pypp.h"
+#include "FormFactorInfLongRipple2.pypp.h"
 #include "IMaterial.pypp.h"
-#include "FTDistribution1DVoigt.pypp.h"
-#include "IntensityDataFunctions.pypp.h"
-#include "FormFactorPrism6.pypp.h"
-#include "IClusteredParticles.pypp.h"
-#include "VariableBinAxis.pypp.h"
+#include "InterferenceFunction2DParaCrystal.pypp.h"
+#include "Particle.pypp.h"
+#include "FormFactorTruncatedSpheroid.pypp.h"
+#include "MesoCrystal.pypp.h"
+#include "FTDistribution2DCauchy.pypp.h"
+#include "IInterferenceFunction.pypp.h"
+#include "FTDistribution1DCosine.pypp.h"
 #include "IParticle.pypp.h"
-#include "DistributionCosine.pypp.h"
-#include "FormFactorHemiEllipsoid.pypp.h"
+#include "FTDistribution1DTriangle.pypp.h"
+#include "FormFactorPrism6.pypp.h"
+#include "HomogeneousMagneticMaterial.pypp.h"
+#include "FormFactorFullSpheroid.pypp.h"
+#include "Transform3D.pypp.h"
+#include "IntensityData.pypp.h"
+#include "DistributionGate.pypp.h"
+#include "ThreadInfo.pypp.h"
 #include "IAxis.pypp.h"
-#include "vector_integer_t.pypp.h"
-#include "IntensityDataIOFactory.pypp.h"
-#include "ParameterDistribution.pypp.h"
-#include "Layer.pypp.h"
-#include "FormFactorPyramid.pypp.h"
 #include "CustomBinAxis.pypp.h"
-#include "FTDistribution2DCone.pypp.h"
-#include "IFTDistribution1D.pypp.h"
-#include "DistributionLorentz.pypp.h"
-#include "IDistribution1D.pypp.h"
-#include "HomogeneousMagneticMaterial.pypp.h"
-#include "FormFactorCuboctahedron.pypp.h"
-#include "cvector_t.pypp.h"
+#include "FormFactorCone6.pypp.h"
+#include "ICloneable.pypp.h"
 #include "PythonInterface_free_functions.pypp.h"
+#include "vector_longinteger_t.pypp.h"
 #include "FormFactorSphereLogNormalRadius.pypp.h"
-#include "FormFactorInfLongRipple1.pypp.h"
-#include "IResolutionFunction2D.pypp.h"
-#include "vector_IFormFactorPtr_t.pypp.h"
-#include "FormFactorFullSphere.pypp.h"
-#include "ParticleLayout.pypp.h"
-#include "FormFactorBox.pypp.h"
-#include "IParameterized.pypp.h"
-#include "Lattice2DIFParameters.pypp.h"
-#include "IFormFactorDecorator.pypp.h"
-#include "InterferenceFunction1DLattice.pypp.h"
-#include "ISample.pypp.h"
-#include "ISampleBuilder.pypp.h"
-#include "PythonInterface_global_variables.pypp.h"
-#include "Beam.pypp.h"
-#include "HomogeneousMaterial.pypp.h"
-#include "ICloneable.pypp.h"
-#include "ParticleCoreShell.pypp.h"
-#include "FormFactorDecoratorDebyeWaller.pypp.h"
-#include "MesoCrystal.pypp.h"
-#include "Lattice1DIFParameters.pypp.h"
-#include "IObserver.pypp.h"
-#include "IntensityData.pypp.h"
-#include "Lattice.pypp.h"
-#include "IInterferenceFunction.pypp.h"
-#include "ParticleInfo.pypp.h"
-#include "Instrument.pypp.h"
-#include "FormFactorInfLongBox.pypp.h"
-#include "FormFactorCone.pypp.h"
-#include "FTDistribution2DGauss.pypp.h"
+#include "FormFactorSphereGaussianRadius.pypp.h"
 #include "FormFactorTruncatedSphere.pypp.h"
-#include "FTDistribution2DVoigt.pypp.h"
-#include "FormFactorGauss.pypp.h"
-#include "InterferenceFunction2DParaCrystal.pypp.h"
-#include "Detector.pypp.h"
-#include "FormFactorInfLongRipple2.pypp.h"
-#include "LatticeBasis.pypp.h"
-#include "ICompositeSample.pypp.h"
-#include "Bin1D.pypp.h"
+#include "IntensityDataFunctions.pypp.h"
+#include "DistributionLorentz.pypp.h"
+#include "ConstKBinAxis.pypp.h"
+#include "ParameterDistribution.pypp.h"
+#include "ParameterPool.pypp.h"
 #include "vector_complex_t.pypp.h"
-#include "DistributionLogNormal.pypp.h"
+#include "SpecularSimulation.pypp.h"
+#include "FormFactorAnisoPyramid.pypp.h"
+#include "FormFactorCuboctahedron.pypp.h"
+#include "MultiLayer.pypp.h"
+#include "FormFactorCone.pypp.h"
+#include "IDistribution1D.pypp.h"
+#include "LayerRoughness.pypp.h"
+#include "VariableBinAxis.pypp.h"
+#include "SimpleSelectionRule.pypp.h"
+#include "FixedBinAxis.pypp.h"
+#include "Lattice2DIFParameters.pypp.h"
+#include "IFormFactor.pypp.h"
+#include "vdouble1d_t.pypp.h"
 #include "IFTDistribution2D.pypp.h"
+#include "Beam.pypp.h"
+#include "FormFactorSphereUniformRadius.pypp.h"
+#include "FTDistribution1DCauchy.pypp.h"
+#include "kvector_t.pypp.h"
+#include "Crystal.pypp.h"
+#include "ISelectionRule.pypp.h"
+#include "vector_integer_t.pypp.h"
 #include "__call_policies.pypp.hpp"
 #include "__convenience.pypp.hpp"
 #include "__call_policies.pypp.hpp"
diff --git a/Core/PythonAPI/src/SpecularSimulation.pypp.cpp b/Core/PythonAPI/src/SpecularSimulation.pypp.cpp
index a3c7cfb2e6a..3d54432181e 100644
--- a/Core/PythonAPI/src/SpecularSimulation.pypp.cpp
+++ b/Core/PythonAPI/src/SpecularSimulation.pypp.cpp
@@ -199,10 +199,20 @@ void register_SpecularSimulation_class(){
                 , getAlphaAxis_function_type( &::SpecularSimulation::getAlphaAxis )
                 , bp::return_value_policy< bp::reference_existing_object >() );
         
+        }
+        { //::SpecularSimulation::getEvanescentWaveIntensity
+        
+            typedef ::OutputData< double > * ( ::SpecularSimulation::*getEvanescentWaveIntensity_function_type)(  ) const;
+            
+            SpecularSimulation_exposer.def( 
+                "getEvanescentWaveIntensity"
+                , getEvanescentWaveIntensity_function_type( &::SpecularSimulation::getEvanescentWaveIntensity )
+                , bp::return_value_policy< bp::manage_new_object >() );
+        
         }
         { //::SpecularSimulation::getLayerRTCoefficients
         
-            typedef ::boost::shared_ptr< const ILayerRTCoefficients > ( ::SpecularSimulation::*getLayerRTCoefficients_function_type)( int,int ) const;
+            typedef ::boost::shared_ptr< const ILayerRTCoefficients > ( ::SpecularSimulation::*getLayerRTCoefficients_function_type)( ::std::size_t,::std::size_t ) const;
             
             SpecularSimulation_exposer.def( 
                 "getLayerRTCoefficients"
@@ -231,32 +241,32 @@ void register_SpecularSimulation_class(){
         }
         { //::SpecularSimulation::getScalarKz
         
-            typedef ::std::vector< std::complex<double> > ( ::SpecularSimulation::*getScalarKz_function_type)( int ) const;
+            typedef ::std::vector< std::complex<double> > ( ::SpecularSimulation::*getScalarKz_function_type)( ::std::size_t ) const;
             
             SpecularSimulation_exposer.def( 
                 "getScalarKz"
                 , getScalarKz_function_type( &::SpecularSimulation::getScalarKz )
-                , ( bp::arg("i_layer")=(int)(0) ) );
+                , ( bp::arg("i_layer") ) );
         
         }
         { //::SpecularSimulation::getScalarR
         
-            typedef ::std::vector< std::complex<double> > ( ::SpecularSimulation::*getScalarR_function_type)( int ) const;
+            typedef ::std::vector< std::complex<double> > ( ::SpecularSimulation::*getScalarR_function_type)( ::std::size_t ) const;
             
             SpecularSimulation_exposer.def( 
                 "getScalarR"
                 , getScalarR_function_type( &::SpecularSimulation::getScalarR )
-                , ( bp::arg("i_layer")=(int)(0) ) );
+                , ( bp::arg("i_layer") ) );
         
         }
         { //::SpecularSimulation::getScalarT
         
-            typedef ::std::vector< std::complex<double> > ( ::SpecularSimulation::*getScalarT_function_type)( int ) const;
+            typedef ::std::vector< std::complex<double> > ( ::SpecularSimulation::*getScalarT_function_type)( ::std::size_t ) const;
             
             SpecularSimulation_exposer.def( 
                 "getScalarT"
                 , getScalarT_function_type( &::SpecularSimulation::getScalarT )
-                , ( bp::arg("i_layer")=(int)(0) ) );
+                , ( bp::arg("i_layer") ) );
         
         }
         { //::SpecularSimulation::prepareSimulation
@@ -296,6 +306,26 @@ void register_SpecularSimulation_class(){
                 , setBeamParameters_function_type( &::SpecularSimulation::setBeamParameters )
                 , ( bp::arg("lambda"), bp::arg("nbins"), bp::arg("alpha_i_min"), bp::arg("alpha_i_max") ) );
         
+        }
+        { //::SpecularSimulation::setEvanescentWaveAxis
+        
+            typedef void ( ::SpecularSimulation::*setEvanescentWaveAxis_function_type)( ::IAxis const & ) ;
+            
+            SpecularSimulation_exposer.def( 
+                "setEvanescentWaveAxis"
+                , setEvanescentWaveAxis_function_type( &::SpecularSimulation::setEvanescentWaveAxis )
+                , ( bp::arg("z_axis") ) );
+        
+        }
+        { //::SpecularSimulation::setEvanescentWaveAxis
+        
+            typedef void ( ::SpecularSimulation::*setEvanescentWaveAxis_function_type)( int,double,double ) ;
+            
+            SpecularSimulation_exposer.def( 
+                "setEvanescentWaveAxis"
+                , setEvanescentWaveAxis_function_type( &::SpecularSimulation::setEvanescentWaveAxis )
+                , ( bp::arg("nbins"), bp::arg("z_min"), bp::arg("z_max") ) );
+        
         }
         { //::SpecularSimulation::setSample
         
diff --git a/Core/Samples/inc/MultiLayer.h b/Core/Samples/inc/MultiLayer.h
index d83bd86794f..f2dc47644d4 100644
--- a/Core/Samples/inc/MultiLayer.h
+++ b/Core/Samples/inc/MultiLayer.h
@@ -128,6 +128,9 @@ public:
     //! returns true if contains magnetic materials and matrix calculations are required
     bool requiresMatrixRTCoefficients() const;
 
+    //! returns layer index corresponding to given global z coordinate
+    size_t zToLayerIndex(double z_value);
+
 protected:
     //! Registers some class members for later access via parameter pool
     virtual void init_parameters();
diff --git a/Core/Samples/src/MultiLayer.cpp b/Core/Samples/src/MultiLayer.cpp
index a8990d019e3..0354b95f4f8 100644
--- a/Core/Samples/src/MultiLayer.cpp
+++ b/Core/Samples/src/MultiLayer.cpp
@@ -279,3 +279,11 @@ bool MultiLayer::requiresMatrixRTCoefficients() const
     return false;
 }
 
+size_t MultiLayer::zToLayerIndex(double z_value)
+{
+    if(z_value < m_layers_z.back()) return m_layers_z.size()-1;
+    std::vector<double>::reverse_iterator top_limit = std::upper_bound(m_layers_z.rbegin(), m_layers_z.rend(), z_value);
+    size_t nbin = m_layers_z.rend() - top_limit;
+    return nbin;
+}
+
diff --git a/Examples/python/simulation/ex05_BeamAndDetector/EvanescentWaveMap.py b/Examples/python/simulation/ex05_BeamAndDetector/EvanescentWaveMap.py
new file mode 100644
index 00000000000..ef1a2760ecd
--- /dev/null
+++ b/Examples/python/simulation/ex05_BeamAndDetector/EvanescentWaveMap.py
@@ -0,0 +1,100 @@
+"""
+Demonstrates how to plot evanescent wave intensity as a function of z-coordinate for multilayer sample.
+Specular simulation geometry.
+"""
+import numpy
+import pylab
+from bornagain import *
+
+
+alpha_i_min, alpha_i_max = 0.0, 2.0  # incoming beam
+
+
+def get_sample():
+    """
+    Build and return the sample representing the layers with correlated roughness.
+    """
+    m_ambience = HomogeneousMaterial("ambience", 0.0, 0.0)
+    m_part_a = HomogeneousMaterial("PartA", 5e-6, 0.0)
+    m_part_b = HomogeneousMaterial("PartB", 10e-6, 0.0)
+    m_substrate = HomogeneousMaterial("substrate", 15e-6, 0.0)
+
+    l_ambience = Layer(m_ambience)
+    l_part_a = Layer(m_part_a, 5.0*nanometer)
+    l_part_b = Layer(m_part_b, 10.0*nanometer)
+    l_substrate = Layer(m_substrate)
+
+    roughness = LayerRoughness()
+    roughness.setSigma(1.0*nanometer)
+    roughness.setHurstParameter(0.3)
+    roughness.setLatteralCorrLength(500.0*nanometer)
+
+    my_sample = MultiLayer()
+
+    # adding layers
+    my_sample.addLayer(l_ambience)
+
+    n_repetitions = 10
+    for i in range(n_repetitions):
+        my_sample.addLayerWithTopRoughness(l_part_a, roughness)
+        my_sample.addLayerWithTopRoughness(l_part_b, roughness)
+
+    my_sample.addLayerWithTopRoughness(l_substrate, roughness)
+    # my_sample.setCrossCorrLength(1e-4)
+
+    return my_sample
+
+
+def get_simulation():
+    """
+    Create and return off-specular simulation with beam and detector defined
+    """
+    simulation = SpecularSimulation()
+    simulation.setBeamParameters(1.54*angstrom, 10, alpha_i_min*degree, alpha_i_max*degree)
+    simulation.setEvanescentWaveAxis(10, -200., 0.0)
+    return simulation
+
+
+def run_simulation():
+    """
+    Run simulation and plot results
+    """
+    sample = get_sample()
+    simulation = get_simulation()
+    simulation.setSample(sample)
+    simulation.runSimulation()
+
+    # plotting results for several selected layers
+    # selected_layers = [0, 1, 20, 21]
+    # alpha_angles = simulation.getAlphaAxis().getBinCenters()
+
+
+    # dpi = 72.
+    # xinch = 1024 / dpi
+    # yinch = 768 / dpi
+    # fig = pylab.figure(figsize=(xinch, yinch))
+    #
+    # nplot = 1
+    # for layer_index in selected_layers:
+    #
+    #     R = []
+    #     for coeff in simulation.getScalarR(layer_index):
+    #         R.append(numpy.abs(coeff))
+    #
+    #     T = []
+    #     for coeff in simulation.getScalarT(layer_index):
+    #         T.append(numpy.abs(coeff))
+    #
+    #     pylab.subplot(2, 2, nplot)
+    #     pylab.ylim(ymax=50.0, ymin=1e-06)
+    #     pylab.semilogy(alpha_angles, R)
+    #     pylab.semilogy(alpha_angles, T)
+    #     pylab.legend(['|R| layer #'+str(layer_index), '|T| layer #'+str(layer_index)], loc='upper right')
+    #     nplot = nplot + 1
+    #
+    #
+    # pylab.show()
+
+
+if __name__ == '__main__':
+    run_simulation()
diff --git a/Examples/python/simulation/ex05_BeamAndDetector/SpecularSimulation.py b/Examples/python/simulation/ex05_BeamAndDetector/SpecularSimulation.py
index c57e3a40d33..d1489277504 100644
--- a/Examples/python/simulation/ex05_BeamAndDetector/SpecularSimulation.py
+++ b/Examples/python/simulation/ex05_BeamAndDetector/SpecularSimulation.py
@@ -2,9 +2,7 @@
 R and T coefficients in multilayer, Specular simulation.
 """
 import numpy
-import matplotlib
 import pylab
-import cmath
 from bornagain import *
 
 
diff --git a/Fit/PythonAPI/src/PythonModule.cpp b/Fit/PythonAPI/src/PythonModule.cpp
index ef8312b19c0..9e316599b4c 100644
--- a/Fit/PythonAPI/src/PythonModule.cpp
+++ b/Fit/PythonAPI/src/PythonModule.cpp
@@ -5,38 +5,38 @@ GCC_DIAG_OFF(missing-field-initializers)
 #include "boost/python.hpp"
 GCC_DIAG_ON(unused-parameter)
 GCC_DIAG_ON(missing-field-initializers)
-#include "IntensityFunctionSqrt.pypp.h"
-#include "MinimizerFactory.pypp.h"
-#include "IMinimizer.pypp.h"
-#include "vector_string_t.pypp.h"
-#include "SquaredFunctionSystematicError.pypp.h"
-#include "IntensityNormalizer.pypp.h"
-#include "IIntensityFunction.pypp.h"
-#include "INamed.pypp.h"
+#include "FitObject.pypp.h"
 #include "IntensityFunctionLog.pypp.h"
-#include "FitSuiteParameters.pypp.h"
-#include "AttFitting.pypp.h"
-#include "FitParameter.pypp.h"
-#include "IntensityScaleAndShiftNormalizer.pypp.h"
+#include "IntensityFunctionSqrt.pypp.h"
+#include "FitStrategyDefault.pypp.h"
 #include "IChiSquaredModule.pypp.h"
+#include "AttFitting.pypp.h"
+#include "FitStrategyAdjustParameters.pypp.h"
 #include "FitStrategyAdjustMinimizer.pypp.h"
-#include "IFitStrategy.pypp.h"
-#include "FitStrategyFixParameters.pypp.h"
+#include "ISquaredFunction.pypp.h"
 #include "SquaredFunctionGaussianError.pypp.h"
+#include "IIntensityFunction.pypp.h"
+#include "IntensityScaleAndShiftNormalizer.pypp.h"
+#include "SquaredFunctionMeanSquaredError.pypp.h"
+#include "FitSuiteParameters.pypp.h"
+#include "IMinimizer.pypp.h"
 #include "IIntensityNormalizer.pypp.h"
+#include "FitParameter.pypp.h"
+#include "INamed.pypp.h"
+#include "FitStrategyReleaseParameters.pypp.h"
 #include "FitSuite.pypp.h"
-#include "FitStrategyAdjustParameters.pypp.h"
+#include "IntensityNormalizer.pypp.h"
+#include "FitStrategyFixParameters.pypp.h"
+#include "FitSuiteObjects.pypp.h"
+#include "SquaredFunctionSystematicError.pypp.h"
+#include "MinimizerFactory.pypp.h"
 #include "ChiSquaredModule.pypp.h"
+#include "SquaredFunctionSimError.pypp.h"
 #include "MinimizerOptions.pypp.h"
-#include "SquaredFunctionDefault.pypp.h"
-#include "SquaredFunctionMeanSquaredError.pypp.h"
-#include "ISquaredFunction.pypp.h"
-#include "FitStrategyDefault.pypp.h"
 #include "AttLimits.pypp.h"
-#include "FitObject.pypp.h"
-#include "FitSuiteObjects.pypp.h"
-#include "SquaredFunctionSimError.pypp.h"
-#include "FitStrategyReleaseParameters.pypp.h"
+#include "SquaredFunctionDefault.pypp.h"
+#include "vector_string_t.pypp.h"
+#include "IFitStrategy.pypp.h"
 
 BOOST_PYTHON_MODULE(libBornAgainFit){
     boost::python::docstring_options doc_options(true, true, false);
diff --git a/Tests/UnitTests/TestCore/MultiLayerTest.h b/Tests/UnitTests/TestCore/MultiLayerTest.h
index f7d7bb6e983..cb02b3b84c1 100644
--- a/Tests/UnitTests/TestCore/MultiLayerTest.h
+++ b/Tests/UnitTests/TestCore/MultiLayerTest.h
@@ -652,5 +652,51 @@ TEST_F(MultiLayerTest, MultiLayerCompositeTest)
 
 }
 
+/*
+
+index_0
+-----------------  0.0
+index_1
+-----------------  -10.0
+index_2
+
+-----------------  -30.0
+index_3
+
+
+-----------------  -60.0
+index_4
+
+*/
+
+TEST_F(MultiLayerTest, MultiLayerZtoIndex)
+{
+    MultiLayer multilayer;
+    HomogeneousMaterial air("air",0,1.0);
+
+    Layer layer0(air, 0*Units::nanometer);
+    Layer layer1(air, 10*Units::nanometer);
+    Layer layer2(air, 20*Units::nanometer);
+    Layer layer3(air, 30*Units::nanometer);
+    Layer layer4(air, 0*Units::nanometer);
+    multilayer.addLayer(layer0);
+    multilayer.addLayer(layer1);
+    multilayer.addLayer(layer2);
+    multilayer.addLayer(layer3);
+    multilayer.addLayer(layer4);
+
+    EXPECT_EQ(0, multilayer.zToLayerIndex(1.0));
+    EXPECT_EQ(0, multilayer.zToLayerIndex(0.0));
+    EXPECT_EQ(1, multilayer.zToLayerIndex(-1.0));
+    EXPECT_EQ(1, multilayer.zToLayerIndex(-9.0));
+    EXPECT_EQ(1, multilayer.zToLayerIndex(-10.0));
+    EXPECT_EQ(2, multilayer.zToLayerIndex(-11.0));
+    EXPECT_EQ(2, multilayer.zToLayerIndex(-30.0));
+    EXPECT_EQ(3, multilayer.zToLayerIndex(-31.0));
+    EXPECT_EQ(3, multilayer.zToLayerIndex(-60.0));
+    EXPECT_EQ(4, multilayer.zToLayerIndex(-61.0));
+
+}
+
 
 #endif
diff --git a/Tests/UnitTests/TestCore/SpecularSimulationTest.h b/Tests/UnitTests/TestCore/SpecularSimulationTest.h
new file mode 100644
index 00000000000..d2777495676
--- /dev/null
+++ b/Tests/UnitTests/TestCore/SpecularSimulationTest.h
@@ -0,0 +1,133 @@
+#ifndef SPECULARSIMULATIONTEST_H
+#define SPECULARSIMULATIONTEST_H
+
+#include "gtest/gtest.h"
+#include "SpecularSimulation.h"
+#include "Exceptions.h"
+#include "FixedBinAxis.h"
+#include "Units.h"
+#include <iostream>
+
+class SpecularSimulationTest : public ::testing::Test
+{
+ protected:
+    SpecularSimulationTest();
+
+    class SampleBuilder : public ISampleBuilder
+    {
+    public:
+        virtual ISample *buildSample() const { return new Layer(); }
+    };
+
+    SampleBuilder_t sample_builder;
+    MultiLayer multilayer;
+};
+
+SpecularSimulationTest::SpecularSimulationTest()
+    : sample_builder(new SampleBuilder)
+{
+    HomogeneousMaterial mat0("ambience", 0.0, 0.0);
+    HomogeneousMaterial mat1("PartA", 5e-6, 0.0);
+    HomogeneousMaterial mat2("substrate", 15e-6, 0.0);
+
+    Layer layer0(mat0);
+    Layer layer1(mat1, 10*Units::nanometer);
+    Layer layer2(mat2);
+
+    multilayer.addLayer(layer0);
+    multilayer.addLayer(layer1);
+    multilayer.addLayer(layer2);
+}
+
+
+TEST_F(SpecularSimulationTest, InitialState)
+{
+    SpecularSimulation sim;
+    ASSERT_THROW( sim.runSimulation(), Exceptions::ClassInitializationException);
+    EXPECT_EQ(NULL, sim.getAlphaAxis());
+    EXPECT_EQ(NULL, sim.getSample());
+    EXPECT_EQ(NULL, sim.getSampleBuilder().get());
+    ASSERT_THROW( sim.getScalarR(0), Exceptions::ClassInitializationException);
+    ASSERT_THROW( sim.getScalarT(0), Exceptions::ClassInitializationException);
+    ASSERT_THROW( sim.getScalarKz(0), Exceptions::ClassInitializationException);
+}
+
+TEST_F(SpecularSimulationTest, CloneOfEmpty)
+{
+    SpecularSimulation sim;
+
+    SpecularSimulation *clone = sim.clone();
+    ASSERT_THROW( clone->runSimulation(), Exceptions::ClassInitializationException);
+    EXPECT_EQ(NULL, clone->getAlphaAxis());
+    EXPECT_EQ(NULL, clone->getSample());
+    EXPECT_EQ(NULL, clone->getSampleBuilder().get());
+    ASSERT_THROW( clone->getScalarR(0), Exceptions::ClassInitializationException);
+    ASSERT_THROW( clone->getScalarT(0), Exceptions::ClassInitializationException);
+    ASSERT_THROW( clone->getScalarKz(0), Exceptions::ClassInitializationException);
+    delete clone;
+}
+
+TEST_F(SpecularSimulationTest, SetBeamParameters)
+{
+    SpecularSimulation sim;
+
+    sim.setBeamParameters(1.0, 10, -2.0, 3.0);
+    EXPECT_EQ(10, sim.getAlphaAxis()->getSize());
+    EXPECT_EQ(-2.0, sim.getAlphaAxis()->getMin());
+    EXPECT_EQ(3.0, sim.getAlphaAxis()->getMax());
+
+    FixedBinAxis axis("axis",2, -1.0, 2.0);
+    sim.setBeamParameters(1.0, axis);
+    EXPECT_EQ(2, sim.getAlphaAxis()->getSize());
+    EXPECT_EQ(-1.0, sim.getAlphaAxis()->getMin());
+    EXPECT_EQ(2.0, sim.getAlphaAxis()->getMax());
+}
+
+TEST_F(SpecularSimulationTest, ConstructSimulation)
+{
+    SpecularSimulation sim;
+    sim.setBeamParameters(1.0, 10, 0.0*Units::degree, 2.0*Units::degree);
+    sim.setSample(multilayer);
+    EXPECT_EQ( size_t(3), dynamic_cast<MultiLayer *>(sim.getSample())->getNumberOfLayers());
+
+    ASSERT_THROW( sim.getScalarR(0), Exceptions::ClassInitializationException);
+    ASSERT_THROW( sim.getScalarT(0), Exceptions::ClassInitializationException);
+    ASSERT_THROW( sim.getScalarKz(0), Exceptions::ClassInitializationException);
+
+    sim.runSimulation();
+    EXPECT_EQ( 10, sim.getScalarR(0).size());
+    EXPECT_EQ( 10, sim.getScalarT(0).size());
+    EXPECT_EQ( 10, sim.getScalarKz(0).size());
+
+    ASSERT_THROW( sim.getScalarR(3), Exceptions::OutOfBoundsException);
+    ASSERT_THROW( sim.getScalarT(3), Exceptions::OutOfBoundsException);
+    ASSERT_THROW( sim.getScalarKz(3), Exceptions::OutOfBoundsException);
+
+}
+
+TEST_F(SpecularSimulationTest, SimulationClone)
+{
+    SpecularSimulation sim;
+    sim.setBeamParameters(1.0, 10, 0.0*Units::degree, 2.0*Units::degree);
+    sim.setSample(multilayer);
+
+    SpecularSimulation *clone = sim.clone();
+
+    EXPECT_EQ( size_t(3), dynamic_cast<MultiLayer *>(clone->getSample())->getNumberOfLayers());
+
+    ASSERT_THROW( clone->getScalarR(0), Exceptions::ClassInitializationException);
+    ASSERT_THROW( clone->getScalarT(0), Exceptions::ClassInitializationException);
+    ASSERT_THROW( clone->getScalarKz(0), Exceptions::ClassInitializationException);
+    delete clone;
+
+    sim.runSimulation();
+
+    clone = sim.clone();
+    EXPECT_EQ( 10, clone->getScalarR(0).size());
+    EXPECT_EQ( 10, clone->getScalarT(0).size());
+    EXPECT_EQ( 10, clone->getScalarKz(0).size());
+
+}
+
+
+#endif
diff --git a/Tests/UnitTests/TestCore/main.cpp b/Tests/UnitTests/TestCore/main.cpp
index 464dd83fd7f..67556643907 100644
--- a/Tests/UnitTests/TestCore/main.cpp
+++ b/Tests/UnitTests/TestCore/main.cpp
@@ -41,6 +41,7 @@
 #include "ConstKBinAxisTest.h"
 #include "CustomBinAxisTest.h"
 #include "IntensityDataFunctionsTest.h"
+#include "SpecularSimulationTest.h"
 
 struct ErrorStreamRedirect {
     ErrorStreamRedirect( std::streambuf * new_buffer )
diff --git a/dev-tools/python-bindings/settings_core.py b/dev-tools/python-bindings/settings_core.py
index 7bb03804977..3d925045524 100644
--- a/dev-tools/python-bindings/settings_core.py
+++ b/dev-tools/python-bindings/settings_core.py
@@ -317,6 +317,11 @@ def ManualClassTunings(mb):
     cl.member_function("getPolarizedIntensityData").call_policies = \
         call_policies.return_value_policy(call_policies.manage_new_object)
     #
+    cl = mb.class_("SpecularSimulation")
+    cl.member_function("setSampleBuilder").include()
+    cl.member_function("getEvanescentWaveIntensity").call_policies = \
+    call_policies.return_value_policy(call_policies.manage_new_object)
+    #
     cl = mb.class_("OffSpecSimulation")
     cl.member_function("setSampleBuilder").include()
     cl.member_function("getOutputData").exclude()
-- 
GitLab