From 8b2e37d4eb4693c6e3c59f5d1e2b9e642246a67f Mon Sep 17 00:00:00 2001
From: Walter Van Herck <w.van.herck@fz-juelich.de>
Date: Fri, 9 Aug 2013 17:34:36 +0200
Subject: [PATCH] Added polarized functionality to LayerStrategyBuilder and
DecouplingApproximationStrategy
---
.../inc/IInterferenceFunctionStrategy.h | 1 +
Core/Algorithms/inc/LayerStrategyBuilder.h | 5 ++
.../src/DecouplingApproximationStrategy.cpp | 31 +++++--
Core/Algorithms/src/LayerStrategyBuilder.cpp | 90 +++++++++++++++----
Core/FormFactors/inc/FormFactors.h | 1 +
Core/Tools/inc/MathFunctions.h | 13 +++
6 files changed, 117 insertions(+), 24 deletions(-)
diff --git a/Core/Algorithms/inc/IInterferenceFunctionStrategy.h b/Core/Algorithms/inc/IInterferenceFunctionStrategy.h
index 9d9f0e54f1b..8fac942d3af 100644
--- a/Core/Algorithms/inc/IInterferenceFunctionStrategy.h
+++ b/Core/Algorithms/inc/IInterferenceFunctionStrategy.h
@@ -21,6 +21,7 @@
#include "Bin.h"
#include "SafePointerVector.h"
#include "LayerStrategyBuilder.h"
+#include "FormFactorDWBAPol.h"
#include <vector>
diff --git a/Core/Algorithms/inc/LayerStrategyBuilder.h b/Core/Algorithms/inc/LayerStrategyBuilder.h
index c23948bbc28..adc1f7e81b6 100644
--- a/Core/Algorithms/inc/LayerStrategyBuilder.h
+++ b/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;
diff --git a/Core/Algorithms/src/DecouplingApproximationStrategy.cpp b/Core/Algorithms/src/DecouplingApproximationStrategy.cpp
index c95771f9fdf..a94a538a6e5 100644
--- a/Core/Algorithms/src/DecouplingApproximationStrategy.cpp
+++ b/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
diff --git a/Core/Algorithms/src/LayerStrategyBuilder.cpp b/Core/Algorithms/src/LayerStrategyBuilder.cpp
index f90af0d98ea..510f9543197 100644
--- a/Core/Algorithms/src/LayerStrategyBuilder.cpp
+++ b/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;
}
-
diff --git a/Core/FormFactors/inc/FormFactors.h b/Core/FormFactors/inc/FormFactors.h
index 2308ab39deb..9c3e92910ab 100644
--- a/Core/FormFactors/inc/FormFactors.h
+++ b/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"
diff --git a/Core/Tools/inc/MathFunctions.h b/Core/Tools/inc/MathFunctions.h
index 86ff5ebe5ad..fa4ab8ed66d 100644
--- a/Core/Tools/inc/MathFunctions.h
+++ b/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
--
GitLab