diff --git a/Core/Algorithms/inc/SpecularMagnetic.h b/Core/Algorithms/inc/SpecularMagnetic.h
new file mode 100644
index 0000000000000000000000000000000000000000..39e4684b5de98d12c107fe526ef8e43aebeb56e1
--- /dev/null
+++ b/Core/Algorithms/inc/SpecularMagnetic.h
@@ -0,0 +1,121 @@
+// ************************************************************************** //
+//
+// BornAgain: simulate and fit scattering at grazing incidence
+//
+//! @file Algorithms/inc/SpecularMagnetic.h
+//! @brief Defines class SpecularMagnetic.
+//!
+//! @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 SPECULARMAGNETIC_H_
+#define SPECULARMAGNETIC_H_
+
+#include "Types.h"
+#include "ISimulation.h"
+#include "MultiLayer.h"
+
+#include <Eigen/Core>
+
+//! Implements the matrix formalism for the calculation of wave amplitudes of
+//! the coherent wave solution in a multilayer with magnetization
+
+class SpecularMagnetic : public ISimulation
+{
+public:
+ SpecularMagnetic() {}
+
+ //! layer coefficients for matrix formalism
+ class LayerMatrixCoeff {
+ public:
+ LayerMatrixCoeff() {}
+ ~LayerMatrixCoeff() {}
+ Eigen::Vector2cd T1() const;
+ Eigen::Vector2cd R1() const;
+ Eigen::Vector2cd T2() const;
+ Eigen::Vector2cd R2() const;
+ // R, T - amplitudes of reflected and transmitted waves
+ Eigen::Vector2cd lambda; // positive eigenvalues of transfer matrix
+ Eigen::Vector2cd kz;
+ Eigen::Vector4cd phi_psi; // boundary values
+ Eigen::Matrix4cd l;
+ Eigen::Matrix4cd T1m;
+ Eigen::Matrix4cd R1m;
+ Eigen::Matrix4cd T2m;
+ Eigen::Matrix4cd R2m;
+ Eigen::Matrix2cd m_scatt_matrix; // scattering matrix
+ complex_t m_a; // polarization independent part
+ complex_t m_b_mag; // magnitude of polarization part
+ complex_t m_bz; // z-part of polarization scattering
+ void calculateTRMatrices();
+ private:
+ void calculateTRWithoutMagnetization();
+ };
+
+ //! multi layer coefficients for matrix formalism
+ class MultiLayerMatrixCoeff
+ {
+ public:
+ LayerMatrixCoeff& operator[](size_t i) { return m_data[i]; }
+ const LayerMatrixCoeff& operator[](size_t i) const { return m_data[i]; }
+ size_t size() const { return m_data.size(); }
+ void clear() { m_data.clear(); }
+ void resize(size_t size) { m_data.resize(size); }
+ private:
+ std::vector<LayerMatrixCoeff > m_data;
+ };
+
+ typedef MultiLayerMatrixCoeff MultiLayerCoeff_t; // set of layer coefficients for matrix formalism
+
+ //! Calculates Fresnel coefficients for given multi layer and kvector
+ void execute(const MultiLayer& sample, const kvector_t& k, MultiLayerCoeff_t& coeff);
+
+private:
+// std::vector<Eigen::Matrix2cd> m_roughness_pmatrices;
+
+ void calculateEigenvalues(const MultiLayer& sample, const kvector_t& k,
+ MultiLayerCoeff_t& coeff) const;
+ void calculateTransferAndBoundary(const MultiLayer& sample,
+ const kvector_t& k, MultiLayerCoeff_t& coeff) const;
+// Eigen::Matrix2cd calculatePMatrix(double sigma_eff,
+// complex_t lambda_lower, complex_t lambda_upper) const;
+ Eigen::Matrix2cd getUnitMatrix() const;
+// complex_t getPMatrixElement(complex_t sigma_lambda) const;
+ void setForNoTransmission(MultiLayerCoeff_t& coeff) const;
+ complex_t getImExponential(complex_t exponent) const;
+};
+
+inline Eigen::Vector2cd SpecularMagnetic::LayerMatrixCoeff::T1() const {
+ Eigen::Vector2cd result;
+ result(0) = T1m.row(2).dot(phi_psi);
+ result(1) = T1m.row(3).dot(phi_psi);
+ return result;
+}
+
+inline Eigen::Vector2cd SpecularMagnetic::LayerMatrixCoeff::R1() const {
+ Eigen::Vector2cd result;
+ result(0) = R1m.row(2).dot(phi_psi);
+ result(1) = R1m.row(3).dot(phi_psi);
+ return result;
+}
+
+inline Eigen::Vector2cd SpecularMagnetic::LayerMatrixCoeff::T2() const {
+ Eigen::Vector2cd result;
+ result(0) = T2m.row(2).dot(phi_psi);
+ result(1) = T2m.row(3).dot(phi_psi);
+ return result;
+}
+
+inline Eigen::Vector2cd SpecularMagnetic::LayerMatrixCoeff::R2() const {
+ Eigen::Vector2cd result;
+ result(0) = R2m.row(2).dot(phi_psi);
+ result(1) = R2m.row(3).dot(phi_psi);
+ return result;
+}
+
+#endif /* SPECULARMAGNETIC_H_ */
diff --git a/Core/Algorithms/src/SpecularMagnetic.cpp b/Core/Algorithms/src/SpecularMagnetic.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f7bdc6526fd20ffd1c345f74816c6ce16c8103ae
--- /dev/null
+++ b/Core/Algorithms/src/SpecularMagnetic.cpp
@@ -0,0 +1,218 @@
+// ************************************************************************** //
+//
+// BornAgain: simulate and fit scattering at grazing incidence
+//
+//! @file SpecularMagnetic.cpp
+//! @brief Implements class SpecularMagnetic.
+//!
+//! @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 "SpecularMagnetic.h"
+
+#include <Eigen/LU>
+
+void SpecularMagnetic::execute(const MultiLayer& sample, const kvector_t& k,
+ MultiLayerCoeff_t& coeff)
+{
+ coeff.clear();
+ coeff.resize(sample.getNumberOfLayers());
+
+ calculateEigenvalues(sample, k, coeff);
+
+ calculateTransferAndBoundary(sample, k, coeff);
+}
+
+void SpecularMagnetic::calculateEigenvalues(const MultiLayer& sample,
+ const kvector_t& k, MultiLayerCoeff_t& coeff) const
+{
+ double mag_k = k.mag();
+ for(size_t i=0; i<coeff.size(); ++i) {
+ coeff[i].m_scatt_matrix = sample.getLayer(i)->getMaterial()->
+ getScatteringMatrix(k);
+ coeff[i].m_a = coeff[i].m_scatt_matrix.trace()/2.0;
+ coeff[i].m_b_mag = std::sqrt(coeff[i].m_a*coeff[i].m_a -
+ coeff[i].m_scatt_matrix.determinant());
+ coeff[i].m_bz = ( coeff[i].m_scatt_matrix(0,0) -
+ coeff[i].m_scatt_matrix(1,1) )/2.0;
+ coeff[i].lambda(0) = std::sqrt(coeff[i].m_a - coeff[i].m_b_mag);
+ coeff[i].lambda(1) = std::sqrt(coeff[i].m_a + coeff[i].m_b_mag);
+ coeff[i].kz = mag_k*coeff[i].lambda;
+ }
+}
+
+void SpecularMagnetic::calculateTransferAndBoundary(const MultiLayer& sample,
+ const kvector_t& k, MultiLayerCoeff_t& coeff) const
+{
+ size_t N = coeff.size();
+ if (coeff[0].lambda == Eigen::Vector2cd::Zero() && N>1) {
+ setForNoTransmission(coeff);
+ return;
+ }
+
+ for (int i=1; i<(int)N; ++i) {
+ coeff[i].calculateTRMatrices();
+ double t = sample.getLayer(i)->getThickness();
+ coeff[i].l =
+ coeff[i].R1m * getImExponential((complex_t)(coeff[i].kz(0)*t)) +
+ coeff[i].T1m * getImExponential((complex_t)(-coeff[i].kz(0)*t)) +
+ coeff[i].R2m * getImExponential((complex_t)(coeff[i].kz(1)*t)) +
+ coeff[i].T2m * getImExponential((complex_t)(-coeff[i].kz(1)*t));
+ }
+}
+
+Eigen::Matrix2cd SpecularMagnetic::getUnitMatrix() const
+{
+ return Eigen::Matrix2cd::Identity();
+}
+
+void SpecularMagnetic::setForNoTransmission(MultiLayerCoeff_t& coeff) const
+{
+ size_t N = coeff.size();
+ for (size_t i=0; i<N; ++i) {
+ coeff[i].phi_psi.setZero();
+ coeff[i].l.setIdentity();
+ coeff[i].T1m = coeff[i].l/4.0;
+ coeff[i].R1m = coeff[i].T1m;
+ coeff[i].T2m = coeff[i].T1m;
+ coeff[i].R2m = coeff[i].T1m;
+ }
+}
+
+void SpecularMagnetic::LayerMatrixCoeff::calculateTRMatrices()
+{
+ if (m_b_mag == 0.0) {
+ calculateTRWithoutMagnetization();
+ return;
+ }
+
+ // T1m:
+ // row 0:
+ T1m(0,0) = (1.0 - m_bz/m_b_mag)/4.0;
+ T1m(0,1) = - m_scatt_matrix(0,1)/4.0/m_b_mag;
+ T1m(0,2) = lambda(0)*(m_bz/m_b_mag - 1.0)/4.0;
+ T1m(0,3) = m_scatt_matrix(0,1)*lambda(0)/4.0/m_b_mag;
+ // row 1:
+ T1m(1,0) = - m_scatt_matrix(1,0)/4.0/m_b_mag;
+ T1m(1,1) = (1.0 + m_bz/m_b_mag)/4.0;
+ T1m(1,2) = m_scatt_matrix(1,0)*lambda(0)/4.0/m_b_mag;
+ T1m(1,3) = - lambda(0)*(m_bz/m_b_mag + 1.0)/4.0;
+ // row 2:
+ T1m(2,0) = -(1.0 - m_bz/m_b_mag)/4.0/lambda(0);
+ T1m(2,1) = m_scatt_matrix(0,1)/4.0/m_b_mag/lambda(0);
+ T1m(2,2) = -(m_bz/m_b_mag - 1.0)/4.0;
+ T1m(2,3) = - m_scatt_matrix(0,1)/4.0/m_b_mag;
+ // row 3:
+ T1m(3,0) = m_scatt_matrix(1,0)/4.0/m_b_mag/lambda(0);
+ T1m(3,1) = -(1.0 + m_bz/m_b_mag)/4.0/lambda(0);
+ T1m(3,2) = - m_scatt_matrix(1,0)/4.0/m_b_mag;
+ T1m(3,3) = (m_bz/m_b_mag + 1.0)/4.0;
+
+ // R1m:
+ // row 0:
+ R1m(0,0) = T1m(0,0);
+ R1m(0,1) = T1m(0,1);
+ R1m(0,2) = -T1m(0,2);
+ R1m(0,3) = -T1m(0,3);
+ // row 1:
+ R1m(1,0) = T1m(1,0);
+ R1m(1,1) = T1m(1,1);
+ R1m(1,2) = -T1m(1,2);
+ R1m(1,3) = -T1m(1,3);
+ // row 2:
+ R1m(2,0) = -T1m(2,0);
+ R1m(2,1) = -T1m(2,1);
+ R1m(2,2) = T1m(2,2);
+ R1m(2,3) = T1m(2,3);
+ // row 3:
+ R1m(3,0) = -T1m(3,0);
+ R1m(3,1) = -T1m(3,1);
+ R1m(3,2) = T1m(3,2);
+ R1m(3,3) = T1m(3,3);
+
+ // T2m:
+ // row 0:
+ T2m(0,0) = (1.0 + m_bz/m_b_mag)/4.0;
+ T2m(0,1) = m_scatt_matrix(0,1)/4.0/m_b_mag;
+ T2m(0,2) = - lambda(1)*(m_bz/m_b_mag + 1.0)/4.0;
+ T2m(0,3) = - m_scatt_matrix(0,1)*lambda(1)/4.0/m_b_mag;
+ // row 1:
+ T2m(1,0) = m_scatt_matrix(1,0)/4.0/m_b_mag;
+ T2m(1,1) = (1.0 - m_bz/m_b_mag)/4.0;
+ T2m(1,2) = - m_scatt_matrix(1,0)*lambda(1)/4.0/m_b_mag;
+ T2m(1,3) = lambda(1)*(m_bz/m_b_mag - 1.0)/4.0;
+ // row 2:
+ T2m(2,0) = -(1.0 + m_bz/m_b_mag)/4.0/lambda(1);
+ T2m(2,1) = - m_scatt_matrix(0,1)/4.0/m_b_mag/lambda(1);
+ T2m(2,2) = (m_bz/m_b_mag + 1.0)/4.0;
+ T2m(2,3) = m_scatt_matrix(0,1)/4.0/m_b_mag;
+ // row 3:
+ T2m(3,0) = - m_scatt_matrix(1,0)/4.0/m_b_mag/lambda(1);
+ T2m(3,1) = -(1.0 - m_bz/m_b_mag)/4.0/lambda(1);
+ T2m(3,2) = m_scatt_matrix(1,0)/4.0/m_b_mag;
+ T2m(3,3) = (1.0 - m_bz/m_b_mag)/4.0;
+
+ // R2m:
+ // row 0:
+ R2m(0,0) = T2m(0,0);
+ R2m(0,1) = T2m(0,1);
+ R2m(0,2) = -T2m(0,2);
+ R2m(0,3) = -T2m(0,3);
+ // row 1:
+ R2m(1,0) = T2m(1,0);
+ R2m(1,1) = T2m(1,1);
+ R2m(1,2) = -T2m(1,2);
+ R2m(1,3) = -T2m(1,3);
+ // row 2:
+ R2m(2,0) = -T2m(2,0);
+ R2m(2,1) = -T2m(2,1);
+ R2m(2,2) = T2m(2,2);
+ R2m(2,3) = T2m(2,3);
+ // row 3:
+ R2m(3,0) = -T2m(3,0);
+ R2m(3,1) = -T2m(3,1);
+ R2m(3,2) = T2m(3,2);
+ R2m(3,3) = T2m(3,3);
+}
+
+void SpecularMagnetic::LayerMatrixCoeff::calculateTRWithoutMagnetization()
+{
+ // T1m:
+ T1m.setZero();
+ T1m(1,1) = 0.5;
+ T1m(1,3) = -std::sqrt(m_a)/2.0;
+ T1m(3,1) = -1.0/(2.0*std::sqrt(m_a));
+ T1m(3,3) = 0.5;
+
+ // R1m:
+ R1m.setZero();
+ R1m(0,0) = 0.5;
+ R1m(0,2) = -std::sqrt(m_a)/2.0;
+ R1m(2,0) = -1.0/(2.0*std::sqrt(m_a));
+ R1m(2,2) = 0.5;
+
+ // T2m:
+ T2m.setZero();
+ T2m(1,1) = 0.5;
+ T2m(1,3) = std::sqrt(m_a)/2.0;
+ T2m(3,1) = 1.0/(2.0*std::sqrt(m_a));
+ T2m(3,3) = 0.5;
+
+ // R2m:
+ R2m.setZero();
+ R2m(0,0) = 0.5;
+ R2m(0,2) = std::sqrt(m_a)/2.0;
+ R2m(2,0) = 1.0/(2.0*std::sqrt(m_a));
+ R2m(2,2) = 0.5;
+}
+
+complex_t SpecularMagnetic::getImExponential(complex_t exponent) const
+{
+ // TODO: add over- and underflow checks!
+ return std::exp(complex_t(0.0, 1.0)*exponent);
+}
diff --git a/Core/Algorithms/src/SpecularMatrix.cpp b/Core/Algorithms/src/SpecularMatrix.cpp
index de1e5eda41368e64536277a91d59ac2f8431b7b2..66cc9fbbac2b3db3c4a33ac32265d16fd1545b51 100644
--- a/Core/Algorithms/src/SpecularMatrix.cpp
+++ b/Core/Algorithms/src/SpecularMatrix.cpp
@@ -141,12 +141,7 @@ Eigen::Matrix2cd SpecularMatrix::calculatePMatrix(double sigma_eff,
Eigen::Matrix2cd SpecularMatrix::getUnitMatrix() const
{
- Eigen::Matrix2cd unit;
- unit(0,0) = 1.0;
- unit(0,1) = 0.0;
- unit(1,0) = 0.0;
- unit(1,1) = 1.0;
- return unit;
+ return Eigen::Matrix2cd::Identity();
}
complex_t SpecularMatrix::getPMatrixElement(complex_t sigma_lambda) const
diff --git a/Core/Core.pro b/Core/Core.pro
index 7a42368f8d04203d0106ef6f88885bc8bb41f475..03cbeed7ab32716f0ee55fb487bc719ee7735524 100644
--- a/Core/Core.pro
+++ b/Core/Core.pro
@@ -76,6 +76,7 @@ SOURCES += \
Algorithms/src/ResolutionFunction2DSimple.cpp \
Algorithms/src/Simulation.cpp \
Algorithms/src/SizeSpacingCorrelationApproximationStrategy.cpp \
+ Algorithms/src/SpecularMagnetic.cpp \
Algorithms/src/SpecularMatrix.cpp \
Algorithms/src/StrategyBuilder.cpp \
\
@@ -224,6 +225,7 @@ HEADERS += \
Algorithms/inc/Simulation.h \
Algorithms/inc/SimulationParameters.h \
Algorithms/inc/SizeSpacingCorrelationApproximationStrategy.h \
+ Algorithms/inc/SpecularMagnetic.h \
Algorithms/inc/SpecularMatrix.h \
Algorithms/inc/StrategyBuilder.h \
Algorithms/inc/ThreadInfo.h \
diff --git a/Core/Samples/inc/HomogeneousMagneticMaterial.h b/Core/Samples/inc/HomogeneousMagneticMaterial.h
index 7709384693629396a834f52216f60c4712b0f081..12a080449cc9ad6382b758c424abf515c5175f6b 100644
--- a/Core/Samples/inc/HomogeneousMagneticMaterial.h
+++ b/Core/Samples/inc/HomogeneousMagneticMaterial.h
@@ -52,7 +52,7 @@ public:
//! Get the scattering matrix from the refractive index, the
//! magnetic field and a given wavevector
- Eigen::Matrix2cd getScatteringMatrix(const kvector_t& k) const;
+ virtual Eigen::Matrix2cd getScatteringMatrix(const kvector_t& k) const;
protected:
kvector_t m_magnetic_field; //!< magnetic field in Tesla
private:
diff --git a/Core/Samples/inc/HomogeneousMaterial.h b/Core/Samples/inc/HomogeneousMaterial.h
index fef46754b910c15b4db9a5a9716de1f9ae326df6..db967ba946f7ccafce38e22d67a1dbd344f942f6 100644
--- a/Core/Samples/inc/HomogeneousMaterial.h
+++ b/Core/Samples/inc/HomogeneousMaterial.h
@@ -23,7 +23,7 @@
class HomogeneousMaterial : public IMaterial
{
- public:
+public:
//! Constructs a material with _name_ and _refractive_index_.
HomogeneousMaterial(const std::string& name,
const complex_t& refractive_index)
@@ -47,7 +47,10 @@ class HomogeneousMaterial : public IMaterial
void setRefractiveIndex(const complex_t &refractive_index)
{ m_refractive_index = refractive_index; }
- protected:
+ //! Get the scattering matrix from the refractive index
+ //! and a given wavevector
+ virtual Eigen::Matrix2cd getScatteringMatrix(const kvector_t& k) const;
+protected:
virtual void print(std::ostream& ostr) const
{
ostr << "HomMat:" << getName() << "<" << this << ">{ " <<
@@ -57,6 +60,16 @@ class HomogeneousMaterial : public IMaterial
complex_t m_refractive_index; //!< complex index of refraction
};
+inline Eigen::Matrix2cd HomogeneousMaterial::getScatteringMatrix(
+ const kvector_t& k) const
+{
+ Eigen::Matrix2cd result;
+ double xy_proj2 = k.magxy2()/k.mag2();
+ complex_t unit_factor = m_refractive_index*m_refractive_index - xy_proj2;
+ result = unit_factor*Eigen::Matrix2cd::Identity();
+ return result;
+}
+
#endif // HOMOGENEOUSMATERIAL_H
diff --git a/Core/Samples/inc/IMaterial.h b/Core/Samples/inc/IMaterial.h
index 88a1e84eb8379ce4f67967ef0202f35452aac90f..c30ac9254f4f449c14be4c3a182dfb4371a4fe3f 100644
--- a/Core/Samples/inc/IMaterial.h
+++ b/Core/Samples/inc/IMaterial.h
@@ -16,15 +16,18 @@
#ifndef IMATERIAL_H
#define IMATERIAL_H
+#include "INamed.h"
+#include "Types.h"
+
#include <string>
#include <iostream>
-#include "INamed.h"
+#include <Eigen/Core>
//! Interface to a named material.
class IMaterial : public INamed
{
- public:
+public:
//! Constructor that sets _name_.
explicit IMaterial(const std::string& name) : INamed(name) {}
@@ -39,7 +42,14 @@ class IMaterial : public INamed
friend std::ostream &operator<<(std::ostream &ostr, const IMaterial &m)
{ m.print(ostr); return ostr; }
- protected:
+ //! Get the scattering matrix from the refractive index
+ //! and a given wavevector
+ virtual Eigen::Matrix2cd getScatteringMatrix(const kvector_t& k) const {
+ (void)k;
+ return Eigen::Matrix2cd::Identity();
+ }
+
+protected:
virtual void print(std::ostream& ostr) const
{ ostr << "IMat:" << getName() << "<" << this << ">"; }
};
diff --git a/Core/Samples/src/HomogeneousMagneticMaterial.cpp b/Core/Samples/src/HomogeneousMagneticMaterial.cpp
index 14ae68748004ac47c42091a52c3d0a256e7cbd23..b2a119738d58e72d06f6523380e8c0a8b44f273a 100644
--- a/Core/Samples/src/HomogeneousMagneticMaterial.cpp
+++ b/Core/Samples/src/HomogeneousMagneticMaterial.cpp
@@ -41,8 +41,10 @@ Eigen::Matrix2cd HomogeneousMagneticMaterial::getScatteringMatrix(
{
Eigen::Matrix2cd result;
double kmag = k.mag();
+ double xy_proj2 = k.magxy2()/k.mag2();
double factor = m_magnetic_prefactor/kmag/kmag;
- result = m_refractive_index*m_refractive_index*m_unit_matrix
+ complex_t unit_factor = m_refractive_index*m_refractive_index - xy_proj2;
+ result = unit_factor*m_unit_matrix
+ factor*m_pauli_operator[0]*m_magnetic_field[0]
+ factor*m_pauli_operator[1]*m_magnetic_field[1]
+ factor*m_pauli_operator[2]*m_magnetic_field[2];