Added polarized functionality to LayerStrategyBuilder and DecouplingApproximationStrategy

Core/Algorithms/inc/IInterferenceFunctionStrategy.h

@@ -21,6 +21,7 @@
 #include "Bin.h"
 #include "SafePointerVector.h"
 #include "LayerStrategyBuilder.h"
+#include "FormFactorDWBAPol.h"
 
 #include <vector>
 

Core/Algorithms/inc/LayerStrategyBuilder.h

@@ -73,6 +73,11 @@ private:
         const ParticleInfo *p_particle_info,
         const IMaterial *p_ambient_material,
         complex_t factor) const;
+    //! Creates formfactor info for single particle in presence of polarization
+    FormFactorInfo *createFormFactorInfoPol(
+        const ParticleInfo *p_particle_info,
+        const IMaterial *p_ambient_material,
+        complex_t factor) const;
 
     SafePointerVector<FormFactorInfo> m_ff_infos;
     SafePointerVector<IInterferenceFunction> m_ifs;

Core/Algorithms/src/DecouplingApproximationStrategy.cpp

@@ -15,6 +15,8 @@
 
 #include "DecouplingApproximationStrategy.h"
 #include "Exceptions.h"
+#include "MathFunctions.h"
+
 #include <cassert>
 #include <iostream>
 
@@ -63,15 +65,26 @@ Eigen::Matrix2d DecouplingApproximationStrategy::evaluatePol(
         const Bin1DCVector& k_f2_bin, double alpha_i, double alpha_f,
         double phi_f) const
 {
-    (void)k_i;
-    (void)k_f1_bin;
-    (void)k_f2_bin;
-    (void)alpha_i;
-    (void)alpha_f;
-    (void)phi_f;
-    throw Exceptions::NotImplementedException(
-            "DecouplingApproximationStrategy::evaluatePol: "
-            "this strategy is not implemented for magnetic systems");
+    Eigen::Matrix2d intensity = Eigen::Matrix2d::Zero();
+    Eigen::Matrix2cd amplitude = Eigen::Matrix2cd::Zero();
+    for (size_t i=0; i<m_ff_infos.size(); ++i) {
+        FormFactorDWBAPol *p_ff_pol = dynamic_cast<FormFactorDWBAPol *>(
+                m_ff_infos[i]->mp_ff);
+        if (!p_ff_pol) {
+            throw Exceptions::ClassInitializationException(
+                    "DecouplingApproximationStrategy::evaluatePol: "
+                    "unpolarized form factor encountered");
+        }
+        Eigen::Matrix2cd  ff = p_ff_pol->evaluatePol(k_i, k_f1_bin, k_f2_bin,
+                alpha_i, alpha_f, phi_f);
+
+        double fraction = m_ff_infos[i]->m_abundance;
+        amplitude += fraction*ff;
+        intensity += fraction*(MathFunctions::Norm(ff));
+    }
+    Eigen::Matrix2d amplitude_norm = MathFunctions::Norm(amplitude);
+    double itf_function = m_ifs[0]->evaluate(k_i-k_f1_bin.getMidPoint());
+    return intensity + amplitude_norm*(itf_function-1.0);
 }
 
 bool DecouplingApproximationStrategy::checkVectorSizes() const

Core/Algorithms/src/LayerStrategyBuilder.cpp

@@ -25,6 +25,7 @@
 #include "PositionParticleInfo.h"
 
 #include <cmath>
+#include <boost/scoped_ptr.hpp>
 
 LayerStrategyBuilder::LayerStrategyBuilder(
         const Layer& decorated_layer, const Simulation& simulation,
@@ -126,9 +127,15 @@ void LayerStrategyBuilder::collectFormFactorInfos()
              0; particle_index<number_of_particles; ++particle_index) {
         const ParticleInfo *p_particle_info =
             p_decoration->getParticleInfo(particle_index);
-        FormFactorInfo *p_ff_info =
-            createFormFactorInfo(p_particle_info, p_layer_material,
-                                 wavevector_scattering_factor);
+        FormFactorInfo *p_ff_info;
+        if (mp_magnetic_coeff_map) {
+            p_ff_info = createFormFactorInfoPol(p_particle_info,
+                    p_layer_material, wavevector_scattering_factor);
+        }
+        else {
+            p_ff_info = createFormFactorInfo(p_particle_info, p_layer_material,
+                    wavevector_scattering_factor);
+        }
         p_ff_info->m_abundance =
             p_decoration->getAbundanceFractionOfParticle(particle_index);
         m_ff_infos.push_back(p_ff_info);
@@ -159,25 +166,21 @@ FormFactorInfo *LayerStrategyBuilder::createFormFactorInfo(
         complex_t factor) const
 {
     FormFactorInfo *p_result = new FormFactorInfo;
-    Particle *p_particle_clone = p_particle_info->getParticle()->clone();
+    boost::scoped_ptr<Particle> P_particle_clone(p_particle_info->getParticle()->clone());
     const Geometry::PTransform3D transform =
         p_particle_info->getPTransform3D();
 
     // formfactor
-    p_particle_clone->setAmbientMaterial(p_ambient_material);
-    IFormFactor *ff_particle = p_particle_clone->createFormFactor();
-    delete p_particle_clone;
-    IFormFactor *ff_transformed(0);
+    P_particle_clone->setAmbientMaterial(p_ambient_material);
+    IFormFactor *ff_particle = P_particle_clone->createFormFactor();
+    IFormFactor *ff_transformed(ff_particle);
     if(transform) {
         ff_transformed = new FormFactorDecoratorTransformation(ff_particle, transform);
-    } else{
-        ff_transformed = ff_particle;
     }
-    IFormFactor *p_ff_framework(0);
+    IFormFactor *p_ff_framework(ff_transformed);
     switch (m_sim_params.me_framework)
     {
     case SimulationParameters::BA:    // Born Approximation
-        p_ff_framework = ff_transformed;
         break;
     case SimulationParameters::DWBA:  // Distorted Wave Born Approximation
     {
@@ -195,8 +198,10 @@ FormFactorInfo *LayerStrategyBuilder::createFormFactorInfo(
     default:
         throw Exceptions::RuntimeErrorException("Framework must be BA or DWBA");
     }
-    FormFactorDecoratorFactor *p_ff =
-        new FormFactorDecoratorFactor(p_ff_framework, factor);
+    IFormFactor *p_ff(p_ff_framework);
+    if ( factor != complex_t(1.0, 0.0) ) {
+        p_ff = new FormFactorDecoratorFactor(p_ff_framework, factor);
+    }
     p_result->mp_ff = p_ff;
     // Other info (position and abundance
     const PositionParticleInfo *p_pos_particle_info =
@@ -210,6 +215,62 @@ FormFactorInfo *LayerStrategyBuilder::createFormFactorInfo(
     return p_result;
 }
 
+FormFactorInfo* LayerStrategyBuilder::createFormFactorInfoPol(
+        const ParticleInfo* p_particle_info,
+        const IMaterial* p_ambient_material, complex_t factor) const
+{
+    FormFactorInfo *p_result = new FormFactorInfo;
+    boost::scoped_ptr<Particle> P_particle_clone(p_particle_info->
+            getParticle()->clone());
+    const Geometry::PTransform3D transform = p_particle_info->getPTransform3D();
+    const IMaterial *p_material = P_particle_clone->getMaterial();
+
+    // formfactor
+    IFormFactor *ff_particle = P_particle_clone->getSimpleFormFactor()->clone();
+    IFormFactor *ff_particle_factor(ff_particle);
+    if ( factor!=complex_t(1.0,0.0) ) {
+        ff_particle_factor = new FormFactorDecoratorFactor(ff_particle, factor);
+    }
+    IFormFactor *ff_transformed(ff_particle_factor);
+    if(transform) {
+        ff_transformed = new FormFactorDecoratorTransformation(
+                ff_particle_factor, transform);
+    }
+    IFormFactor *p_ff_framework(ff_transformed);
+    switch (m_sim_params.me_framework)
+    {
+    case SimulationParameters::BA:    // Born Approximation
+        break;
+    case SimulationParameters::DWBA:  // Distorted Wave Born Approximation
+    {
+        if (!mp_magnetic_coeff_map) {
+            throw Exceptions::ClassInitializationException(
+                    "Magnetic coefficients are necessary for DWBA");
+        }
+        FormFactorDWBAPol *p_dwba_ff_pol =
+            new FormFactorDWBAPol(ff_transformed);
+        p_dwba_ff_pol->setRTInfo(*mp_magnetic_coeff_map);
+        p_dwba_ff_pol->setMaterial(p_material);
+        p_dwba_ff_pol->setAmbientMaterial(p_ambient_material);
+        p_ff_framework = p_dwba_ff_pol;
+        break;
+    }
+    default:
+        throw Exceptions::RuntimeErrorException("Framework must be BA or DWBA");
+    }
+    p_result->mp_ff = p_ff_framework;
+    // Other info (position and abundance)
+    const PositionParticleInfo *p_pos_particle_info =
+        dynamic_cast<const PositionParticleInfo *>(p_particle_info);
+    if (p_pos_particle_info) {
+        kvector_t position = p_pos_particle_info->getPosition();
+        p_result->m_pos_x = position.x();
+        p_result->m_pos_y = position.y();
+    }
+    p_result->m_abundance = p_particle_info->getAbundance();
+    return p_result;
+}
+
 // =============================================================================
 // Implementation of FormFactorInfo
 // =============================================================================
@@ -226,4 +287,3 @@ FormFactorInfo* FormFactorInfo::clone() const
     return p_result;
 }
 
-

Core/FormFactors/inc/FormFactors.h

@@ -28,6 +28,7 @@
 #include "FormFactorDecoratorTransformation.h"
 #include "FormFactorDWBA.h"
 #include "FormFactorDWBAConstZ.h"
+#include "FormFactorDWBAPol.h"
 #include "FormFactorEllipsoid.h"
 #include "FormFactorFullSphere.h"
 #include "FormFactorFullSpheroid.h"

Core/Tools/inc/MathFunctions.h

@@ -29,6 +29,8 @@
 #include "gsl/gsl_sf_expint.h"
 #include "gsl/gsl_integration.h"
 
+#include <Eigen/Core>
+
 //! Various mathematical functions.
 
 namespace MathFunctions
@@ -85,6 +87,8 @@ complex_t FastCos(const complex_t &x);
 //! simultaneous complex sine and cosine calculations
 void FastSinCos(const complex_t &x, complex_t &xsin, complex_t &xcos);
 
+Eigen::Matrix2d Norm(Eigen::Matrix2cd &M);
+
 } // Namespace MathFunctions
 
 inline double MathFunctions::GenerateNormalRandom(double average, double std_dev)
@@ -183,6 +187,15 @@ inline void MathFunctions::FastSinCos(const complex_t &x,
     xcos = complex_t( cosa*coshb, -sina*sinhb );
 }
 
+inline Eigen::Matrix2d MathFunctions::Norm(Eigen::Matrix2cd &M) {
+    Eigen::Matrix2d result;
+    result(0,0) = std::norm((complex_t)M(0,0));
+    result(0,1) = std::norm((complex_t)M(0,1));
+    result(1,0) = std::norm((complex_t)M(1,0));
+    result(1,1) = std::norm((complex_t)M(1,1));
+    return result;
+}
+
 #endif // MATHFUNCTIONS_H