Added polarization density matrix to Beam and fixed performance test

App/src/TestPerformance.cpp

@@ -22,6 +22,8 @@
 #include "MaterialManager.h"
 #include "OpticalFresnel.h"
 #include "SampleFactory.h"
+#include "SpecularMatrix.h"
+#include "SampleBuilderFactory.h"
 
 #include "TSystem.h"
 #include "TDatime.h"
@@ -32,7 +34,6 @@
 #include <sstream>
 #include <iomanip>
 #include <time.h>
-#include "SpecularMatrix.h"
 
 TestPerformance::TestPerformance()
 {
@@ -210,7 +211,8 @@ void PerfTest_Pyramid::initialise(ProgramOptions *p_options)
     IFunctionalTest::initialise(p_options);
     // sample
     if(m_sample) delete m_sample;
-    m_sample = dynamic_cast<MultiLayer *>(SampleFactory::createSample("IsGISAXS9_Pyramid"));
+    SampleBuilderFactory factory;
+    m_sample = dynamic_cast<MultiLayer *>(factory.createSample("isgisaxs09"));
 
     // simulation
     if(m_simulation) delete m_simulation;
@@ -234,7 +236,8 @@ void PerfTest_RotatedPyramid::initialise(ProgramOptions *p_options)
     IFunctionalTest::initialise(p_options);
     // sample
     if(m_sample) delete m_sample;
-    m_sample = dynamic_cast<MultiLayer *>(SampleFactory::createSample("IsGISAXS9_RotatedPyramid"));
+    SampleBuilderFactory factory;
+    m_sample = dynamic_cast<MultiLayer *>(factory.createSample("isgisaxs09_rotated"));
 
     // simulation
     if(m_simulation) delete m_simulation;

Core/Algorithms/inc/Beam.h

@@ -20,6 +20,8 @@
 #include "Types.h"
 #include "IParameterized.h"
 
+#include <Eigen/Core>
+
 //! Ideal collimated Beam defined by wavelength, direction and intensity.
 
 class BA_CORE_API_ Beam : public IParameterized
@@ -31,15 +33,35 @@ class BA_CORE_API_ Beam : public IParameterized
 
     virtual ~Beam() {}
 
+    //! Get the value of the wavevector
     cvector_t getCentralK() const { return m_central_k; }
 
+    //! Sets the value of the incoming wavevector
     void setCentralK(const cvector_t& k_i);
+
+    //! Sets the value of the incoming wavevector in terms of wavelength
+    //! and incoming angles
     void setCentralK(double lambda, double alpha_i, double phi_i);
 
+    //! Gets the beam intensity in neutrons/sec
     double getIntensity() const { return m_intensity; }
 
+    //! Sets the beam intensity in neutrons/sec
     void setIntensity(double intensity) { m_intensity = intensity; }
 
+    //! Gets the polarization density matrix (in spin basis along z-axis)
+    Eigen::Matrix2cd getPolarization() const { return m_polarization; }
+
+    //! Sets the polarization density matrix (in spin basis along z-axis)
+    void setPolarization(const Eigen::Matrix2cd &polarization);
+
+    //! Sets the polarization density matrix to a value representing
+    //! a mixed ensemble with the given fraction of positive z spin
+    void SetSpinUpFraction(double up_fraction);
+
+    //! Checks if the given matrix can represent a physical density matrix
+    bool checkPolarization(const Eigen::Matrix2cd &polarization) const;
+
  protected:
     //! Registers some class members for later access via parameter pool
     virtual void init_parameters();
@@ -48,8 +70,12 @@ class BA_CORE_API_ Beam : public IParameterized
     //! swap function
     void swapContent(Beam& other);
 
-    cvector_t m_central_k;
-    double m_intensity;
+    //! Initialize polarization (for constructors)
+    void initPolarization();
+
+    cvector_t m_central_k;  //!< incoming wavevector
+    double m_intensity;     //!< beam intensity (neutrons/sec)
+    Eigen::Matrix2cd m_polarization; //!< polarization density matrix
 };
 
 #endif /* BEAM_H_ */

Core/Algorithms/inc/SpecularMatrix.h

@@ -20,7 +20,7 @@
 #include "ISimulation.h"
 #include "MultiLayer.h"
 
-#include "Eigen/Core"
+#include <Eigen/Core>
 
 //! Implements the matrix formalism for the calculation of wave amplitudes of
 //! the coherent wave solution in a multilayer
@@ -42,7 +42,6 @@ public:
        complex_t kz;
        Eigen::Vector2cd phi_psi;
        Eigen::Matrix2cd l;
-       Eigen::Matrix2cd l_accumulated;
    };
 
    //! multi layer coefficients for matrix formalism
@@ -54,8 +53,6 @@ public:
        inline size_t size() const { return m_data.size(); }
        inline void clear() { m_data.clear(); }
        inline void resize(size_t size) { m_data.resize(size); }
-       complex_t R; // total reflection coefficient
-       Eigen::Matrix2cd L;
    private:
        std::vector<LayerMatrixCoeff > m_data;
    };

Core/Algorithms/src/Beam.cpp

@@ -14,14 +14,19 @@
 // ************************************************************************** //
 
 #include "Beam.h"
+#include "Exceptions.h"
+
+#include <Eigen/LU>
 
 Beam::Beam() : m_central_k(), m_intensity(1.0)
 {
     setName("Beam");
     init_parameters();
+    initPolarization();
 }
 
-Beam::Beam(const Beam& other) : IParameterized(), m_central_k(other.m_central_k), m_intensity(other.m_intensity)
+Beam::Beam(const Beam& other) : IParameterized(), m_central_k(other.m_central_k)
+  , m_intensity(other.m_intensity), m_polarization(other.m_polarization)
 {
     setName(other.getName());
     init_parameters();
@@ -37,12 +42,6 @@ Beam& Beam::operator=(const Beam& other)
     return *this;
 }
 
-void Beam::swapContent(Beam& other)
-{
-    std::swap(this->m_central_k, other.m_central_k);
-    std::swap(this->m_intensity, other.m_intensity);
-}
-
 void Beam::setCentralK(const cvector_t& k_i)
 {
     m_central_k = k_i;
@@ -55,10 +54,54 @@ void Beam::setCentralK(double lambda, double alpha_i, double phi_i)
     m_central_k = k_i;
 }
 
+void Beam::setPolarization(const Eigen::Matrix2cd &polarization)
+{
+    if (checkPolarization(polarization)) {
+        m_polarization = polarization;
+    }
+    else {
+        throw Exceptions::ClassInitializationException("The given polarization"
+                " matrix cannot represent a physical density matrix");
+    }
+}
+
+void Beam::SetSpinUpFraction(double up_fraction)
+{
+    if (up_fraction < 0.0 || up_fraction > 1.0) {
+        throw Exceptions::ClassInitializationException("The fraction of spin-up"
+                " states must be between 0.0 and 1.0");
+    }
+    m_polarization.setZero();
+    m_polarization(0,0) = up_fraction;
+    m_polarization(1,1) = 1.0 - up_fraction;
+}
+
+bool Beam::checkPolarization(const Eigen::Matrix2cd &polarization) const
+{
+    if (std::imag( (complex_t)polarization(0,0) ) != 0.0) return false;
+    if (polarization(0,0)+polarization(1,1) != 1.0) return false;
+    if (polarization(0,1) != std::conj( (complex_t)polarization(1,0))) return false;
+    if (std::abs( polarization.determinant() ) < 0.0) return false;
+
+    return true;
+}
+
 void Beam::init_parameters()
 {
     clearParameterPool();
     registerParameter("intensity", &m_intensity);
 }
 
+void Beam::swapContent(Beam& other)
+{
+    std::swap(this->m_central_k, other.m_central_k);
+    std::swap(this->m_intensity, other.m_intensity);
+    std::swap(this->m_polarization, other.m_polarization);
+}
 
+void Beam::initPolarization()
+{
+    m_polarization.setZero();
+    m_polarization(0,0) = 0.5;
+    m_polarization(1,1) = 0.5;
+}

Core/Samples/inc/DiffuseParticleInfo.h

@@ -29,7 +29,7 @@ class DiffuseParticleInfo: public ParticleInfo
         double depth=0, double abundance=0)
         : ParticleInfo(p_particle, transform, depth, abundance)
         , m_number_per_meso(0.0)
-    {}
+        , m_height_range(0.0) {}
 
     virtual ~DiffuseParticleInfo() {}
 

dev-tools/log/perf_history.txt

@@ -248,3 +248,4 @@
 2013-05-07 16:05:21 | jcnsopc74  | macosx64, 2800 MHz      | 294118      | 11.976      | 11.976      | 3.63636     | 1.05263e+06 | 
 2013-05-07 16:44:53 | jcnsopc74  | macosx64, 2800 MHz      | 298507      | 11.9048     | 11.8343     | 3.63636     | 1.05263e+06 | 
 2013-07-12 13:13:46 | jcnsopc74  | macosx64, 2800 MHz      | 298507      | 23.5294     | 17.2414     | 5.71429     | 800000      | 
+2013-07-17 17:19:14 | jcnsopc74  | macosx64, 2800 MHz      | 298507      | 23.8095     | 23.5294     | 6.06061     | 800000      |