diff --git a/Core/Export/SampleToPython.cpp b/Core/Export/SampleToPython.cpp
index 07301b1735112f5822230d9b284648c32cbee560..3a3b544cdf005206be7fb396940b6f1431e41e36 100644
--- a/Core/Export/SampleToPython.cpp
+++ b/Core/Export/SampleToPython.cpp
@@ -155,10 +155,10 @@ const std::map<MATERIAL_TYPES, std::string> factory_names{
std::string SampleToPython::defineMaterials() const {
const auto themap = m_materials->materialMap();
if (themap.empty())
- return "# No Materials.\n\n";
+ return "";
std::ostringstream result;
result << std::setprecision(12);
- result << indent() << "# Define Materials\n";
+ result << indent() << "# Define materials\n";
std::set<std::string> visitedMaterials;
for (auto it : themap) {
const std::string& key = it.first;
diff --git a/Core/Legacy/MatrixRTCoefficients_v1.cpp b/Core/Legacy/MatrixRTCoefficients_v1.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7f297e3ab5a23907a8d6b1f38748d87b9681448b
--- /dev/null
+++ b/Core/Legacy/MatrixRTCoefficients_v1.cpp
@@ -0,0 +1,303 @@
+// ************************************************************************************************
+//
+// BornAgain: simulate and fit scattering at grazing incidence
+//
+//! @file Sample/RT/MatrixRTCoefficients_v1.cpp
+//! @brief Implements class MatrixRTCoefficients_v1.
+//!
+//! @homepage http://www.bornagainproject.org
+//! @license GNU General Public License v3 or higher (see COPYING)
+//! @copyright Forschungszentrum Jülich GmbH 2018
+//! @authors Scientific Computing Group at MLZ (see CITATION, AUTHORS)
+//
+// ************************************************************************************************
+
+#include "Core/Legacy/MatrixRTCoefficients_v1.h"
+
+MatrixRTCoefficients_v1* MatrixRTCoefficients_v1::clone() const {
+ return new MatrixRTCoefficients_v1(*this);
+}
+
+Eigen::Vector2cd MatrixRTCoefficients_v1::T1plus() const {
+ Eigen::Vector2cd result;
+ Eigen::Matrix<complex_t, 4, 1> m = T1m * phi_psi_plus;
+ result(0) = m(2);
+ result(1) = m(3);
+ if (lambda(0) == 0.0 && result == Eigen::Vector2cd::Zero())
+ result(0) = 0.5;
+ return result;
+}
+
+Eigen::Vector2cd MatrixRTCoefficients_v1::R1plus() const {
+ Eigen::Vector2cd result;
+ Eigen::Matrix<complex_t, 4, 1> m = R1m * phi_psi_plus;
+ result(0) = m(2);
+ result(1) = m(3);
+ Eigen::Matrix<complex_t, 4, 1> mT = T1m * phi_psi_plus;
+ if (lambda(0) == 0.0 && mT(2) == 0.0 && mT(3) == 0.0)
+ result(0) = -0.5;
+ return result;
+}
+
+Eigen::Vector2cd MatrixRTCoefficients_v1::T2plus() const {
+ Eigen::Vector2cd result;
+ Eigen::Matrix<complex_t, 4, 1> m = T2m * phi_psi_plus;
+ result(0) = m(2);
+ result(1) = m(3);
+ if (lambda(1) == 0.0 && result == Eigen::Vector2cd::Zero())
+ result(0) = 0.5;
+ return result;
+}
+
+Eigen::Vector2cd MatrixRTCoefficients_v1::R2plus() const {
+ Eigen::Vector2cd result;
+ Eigen::Matrix<complex_t, 4, 1> m = R2m * phi_psi_plus;
+ result(0) = m(2);
+ result(1) = m(3);
+ Eigen::Matrix<complex_t, 4, 1> mT = T2m * phi_psi_plus;
+ if (lambda(1) == 0.0 && mT(2) == 0.0 && mT(3) == 0.0)
+ result(0) = -0.5;
+ return result;
+}
+
+Eigen::Vector2cd MatrixRTCoefficients_v1::T1min() const {
+ Eigen::Vector2cd result;
+ Eigen::Matrix<complex_t, 4, 1> m = T1m * phi_psi_min;
+ result(0) = m(2);
+ result(1) = m(3);
+ if (lambda(0) == 0.0 && result == Eigen::Vector2cd::Zero())
+ result(1) = 0.5;
+ return result;
+}
+
+Eigen::Vector2cd MatrixRTCoefficients_v1::R1min() const {
+ Eigen::Vector2cd result;
+ Eigen::Matrix<complex_t, 4, 1> m = R1m * phi_psi_min;
+ result(0) = m(2);
+ result(1) = m(3);
+ Eigen::Matrix<complex_t, 4, 1> mT = T1m * phi_psi_min;
+ if (lambda(0) == 0.0 && mT(2) == 0.0 && mT(3) == 0.0)
+ result(1) = -0.5;
+ return result;
+}
+
+Eigen::Vector2cd MatrixRTCoefficients_v1::T2min() const {
+ Eigen::Vector2cd result;
+ Eigen::Matrix<complex_t, 4, 1> m = T2m * phi_psi_min;
+ result(0) = m(2);
+ result(1) = m(3);
+ if (lambda(1) == 0.0 && result == Eigen::Vector2cd::Zero())
+ result(1) = 0.5;
+ return result;
+}
+
+Eigen::Vector2cd MatrixRTCoefficients_v1::R2min() const {
+ Eigen::Vector2cd result;
+ Eigen::Matrix<complex_t, 4, 1> m = R2m * phi_psi_min;
+ result(0) = m(2);
+ result(1) = m(3);
+ Eigen::Matrix<complex_t, 4, 1> mT = T2m * phi_psi_min;
+ if (lambda(1) == 0.0 && mT(2) == 0.0 && mT(3) == 0.0)
+ result(1) = -0.5;
+ return result;
+}
+
+void MatrixRTCoefficients_v1::calculateTRMatrices() {
+ if (m_b_mag == 0.0) {
+ calculateTRWithoutMagnetization();
+ return;
+ }
+
+ if (lambda(0) == 0.0) {
+ complex_t ikt = mul_I(m_kt);
+ // Lambda1 component contained only in T1m (R1m=0)
+ // row 0:
+ T1m(0, 0) = (1.0 - m_bz / m_b_mag) / 2.0;
+ T1m(0, 1) = -m_scatt_matrix(0, 1) / 2.0 / m_b_mag;
+ // row 1:
+ T1m(1, 0) = -m_scatt_matrix(1, 0) / 2.0 / m_b_mag;
+ T1m(1, 1) = (1.0 + m_bz / m_b_mag) / 2.0;
+ // row 2:
+ T1m(2, 0) = ikt * (1.0 - m_bz / m_b_mag) / 2.0;
+ T1m(2, 1) = -ikt * m_scatt_matrix(0, 1) / 2.0 / m_b_mag;
+ T1m(2, 2) = T1m(0, 0);
+ T1m(2, 3) = T1m(0, 1);
+ // row 3:
+ T1m(3, 0) = -ikt * m_scatt_matrix(1, 0) / 2.0 / m_b_mag;
+ T1m(3, 1) = ikt * (1.0 + m_bz / m_b_mag) / 2.0;
+ T1m(3, 2) = T1m(1, 0);
+ T1m(3, 3) = T1m(1, 1);
+ } else {
+ // 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);
+ }
+
+ if (lambda(1) == 0.0) {
+ complex_t ikt = mul_I(m_kt);
+ // Lambda2 component contained only in T2m (R2m=0)
+ // row 0:
+ T2m(0, 0) = (1.0 + m_bz / m_b_mag) / 2.0;
+ T2m(0, 1) = m_scatt_matrix(0, 1) / 2.0 / m_b_mag;
+ // row 1:
+ T2m(1, 0) = m_scatt_matrix(1, 0) / 2.0 / m_b_mag;
+ T2m(1, 1) = (1.0 - m_bz / m_b_mag) / 2.0;
+ // row 2:
+ T2m(2, 0) = ikt * (1.0 + m_bz / m_b_mag) / 2.0;
+ T2m(2, 1) = ikt * m_scatt_matrix(0, 1) / 2.0 / m_b_mag;
+ T2m(2, 2) = T2m(0, 0);
+ T2m(2, 3) = T2m(0, 1);
+ // row 3:
+ T2m(3, 0) = ikt * m_scatt_matrix(1, 0) / 2.0 / m_b_mag;
+ T2m(3, 1) = ikt * (1.0 - m_bz / m_b_mag) / 2.0;
+ T2m(3, 2) = T2m(1, 0);
+ T2m(3, 3) = T2m(1, 1);
+ } else {
+ // 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 MatrixRTCoefficients_v1::calculateTRWithoutMagnetization() {
+ T1m.setZero();
+ R1m.setZero();
+ T2m.setZero();
+ R2m.setZero();
+
+ if (m_a == 0.0) {
+ // Spin down component contained only in T1 (R1=0)
+ T1m(1, 1) = 1.0;
+ T1m(3, 1) = mul_I(m_kt);
+ T1m(3, 3) = 1.0;
+
+ // Spin up component contained only in T2 (R2=0)
+ T2m(0, 0) = 1.0;
+ T2m(2, 0) = mul_I(m_kt);
+ T2m(2, 2) = 1.0;
+ return;
+ }
+
+ // T1m:
+ 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(1, 1) = 0.5;
+ R1m(1, 3) = std::sqrt(m_a) / 2.0;
+ R1m(3, 1) = 1.0 / (2.0 * std::sqrt(m_a));
+ R1m(3, 3) = 0.5;
+
+ // T2m:
+ T2m(0, 0) = 0.5;
+ T2m(0, 2) = -std::sqrt(m_a) / 2.0;
+ T2m(2, 0) = -1.0 / (2.0 * std::sqrt(m_a));
+ T2m(2, 2) = 0.5;
+
+ // R2m:
+ 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;
+}
+
+void MatrixRTCoefficients_v1::initializeBottomLayerPhiPsi() {
+ if (m_b_mag == 0.0) {
+ phi_psi_min << 0.0, -std::sqrt(m_a), 0.0, 1.0;
+ phi_psi_plus << -std::sqrt(m_a), 0.0, 1.0, 0.0;
+ return;
+ }
+ // First basis vector that has no upward going wave amplitude
+ phi_psi_min(0) = m_scatt_matrix(0, 1) * (lambda(0) - lambda(1)) / 2.0 / m_b_mag;
+ phi_psi_min(1) = (m_bz * (lambda(1) - lambda(0)) / m_b_mag - lambda(1) - lambda(0)) / 2.0;
+ phi_psi_min(2) = 0.0;
+ phi_psi_min(3) = 1.0;
+
+ // Second basis vector that has no upward going wave amplitude
+ phi_psi_plus(0) = -(m_scatt_matrix(0, 0) + lambda(0) * lambda(1)) / (lambda(0) + lambda(1));
+ phi_psi_plus(1) = m_scatt_matrix(1, 0) * (lambda(0) - lambda(1)) / 2.0 / m_b_mag;
+ phi_psi_plus(2) = 1.0;
+ phi_psi_plus(3) = 0.0;
+}
diff --git a/Core/Legacy/MatrixRTCoefficients_v1.h b/Core/Legacy/MatrixRTCoefficients_v1.h
new file mode 100644
index 0000000000000000000000000000000000000000..b983370cc21c5a9af1b90af76132cf3de7db5f6e
--- /dev/null
+++ b/Core/Legacy/MatrixRTCoefficients_v1.h
@@ -0,0 +1,69 @@
+// ************************************************************************************************
+//
+// BornAgain: simulate and fit scattering at grazing incidence
+//
+//! @file Sample/RT/MatrixRTCoefficients_v1.h
+//! @brief Defines class MatrixRTCoefficients_v1.
+//!
+//! @homepage http://www.bornagainproject.org
+//! @license GNU General Public License v3 or higher (see COPYING)
+//! @copyright Forschungszentrum Jülich GmbH 2018
+//! @authors Scientific Computing Group at MLZ (see CITATION, AUTHORS)
+//
+// ************************************************************************************************
+
+#ifndef BORNAGAIN_CORE_LEGACY_MATRIXRTCOEFFICIENTS_V1_H
+#define BORNAGAIN_CORE_LEGACY_MATRIXRTCOEFFICIENTS_V1_H
+
+#include "Sample/RT/ILayerRTCoefficients.h"
+
+//! Specular reflection and transmission coefficients in a layer in case
+//! of 2x2 matrix interactions between the layers and the scattered particle.
+//! @ingroup algorithms_internal
+
+class MatrixRTCoefficients_v1 : public ILayerRTCoefficients {
+public:
+ MatrixRTCoefficients_v1() : m_kt(0.0) {}
+ virtual ~MatrixRTCoefficients_v1() {}
+
+ virtual MatrixRTCoefficients_v1* clone() const;
+
+ //! The following functions return the transmitted and reflected amplitudes
+ //! for different incoming beam polarizations and eigenmodes
+ virtual Eigen::Vector2cd T1plus() const;
+ virtual Eigen::Vector2cd R1plus() const;
+ virtual Eigen::Vector2cd T2plus() const;
+ virtual Eigen::Vector2cd R2plus() const;
+ virtual Eigen::Vector2cd T1min() const;
+ virtual Eigen::Vector2cd R1min() const;
+ virtual Eigen::Vector2cd T2min() const;
+ virtual Eigen::Vector2cd R2min() const;
+ //! Returns z-part of the two wavevector eigenmodes
+ virtual Eigen::Vector2cd getKz() const { return kz; }
+
+ void calculateTRMatrices();
+ void calculateTRWithoutMagnetization();
+ void initializeBottomLayerPhiPsi();
+
+ // NOTE: exceptionally, this class has member variables without prefix m_
+ Eigen::Vector2cd kz; //!< z-part of the two wavevector eigenmodes
+ Eigen::Vector2cd lambda; //!< positive eigenvalues of transfer matrix
+ Eigen::Vector4cd phi_psi_plus; //!< boundary values for up-polarization
+ Eigen::Vector4cd phi_psi_min; //!< boundary values for down-polarization
+ Eigen::Matrix4cd T1m; //!< matrix selecting the transmitted part of
+ //!< the first eigenmode
+ Eigen::Matrix4cd R1m; //!< matrix selecting the reflected part of
+ //!< the first eigenmode
+ Eigen::Matrix4cd T2m; //!< matrix selecting the transmitted part of
+ //!< the second eigenmode
+ Eigen::Matrix4cd R2m; //!< matrix selecting the reflected part of
+ //!< the second eigenmode
+ Eigen::Matrix2cd m_scatt_matrix; //!< scattering matrix
+ complex_t m_a; //!< polarization independent part of scattering matrix
+ complex_t m_b_mag; //!< magnitude of magnetic interaction term
+ complex_t m_bz; //!< z-part of magnetic interaction term
+ double m_kt; //!< wavevector length times thickness of layer for use when
+ //!< lambda=0
+};
+
+#endif // BORNAGAIN_CORE_LEGACY_MATRIXRTCOEFFICIENTS_V1_H
diff --git a/Sample/RT/MatrixRTCoefficients_v2.cpp b/Core/Legacy/MatrixRTCoefficients_v2.cpp
similarity index 98%
rename from Sample/RT/MatrixRTCoefficients_v2.cpp
rename to Core/Legacy/MatrixRTCoefficients_v2.cpp
index f976fcf69574ab1080f53797fa7b06d305a1a799..24f23b5b5a74384c7885b86197784022a164dd5d 100644
--- a/Sample/RT/MatrixRTCoefficients_v2.cpp
+++ b/Core/Legacy/MatrixRTCoefficients_v2.cpp
@@ -12,7 +12,7 @@
//
// ************************************************************************************************
-#include "Sample/RT/MatrixRTCoefficients_v2.h"
+#include "Core/Legacy/MatrixRTCoefficients_v2.h"
namespace {
Eigen::Vector2cd waveVector(const Eigen::Matrix4cd& frob_matrix,
diff --git a/Sample/RT/MatrixRTCoefficients_v2.h b/Core/Legacy/MatrixRTCoefficients_v2.h
similarity index 91%
rename from Sample/RT/MatrixRTCoefficients_v2.h
rename to Core/Legacy/MatrixRTCoefficients_v2.h
index dfdc242bcb1830d9ebb902242f15733711ae11dd..2146ea36eace9cd1286c4a082df35c47f4321412 100644
--- a/Sample/RT/MatrixRTCoefficients_v2.h
+++ b/Core/Legacy/MatrixRTCoefficients_v2.h
@@ -12,8 +12,8 @@
//
// ************************************************************************************************
-#ifndef BORNAGAIN_SAMPLE_RT_MATRIXRTCOEFFICIENTS_V2_H
-#define BORNAGAIN_SAMPLE_RT_MATRIXRTCOEFFICIENTS_V2_H
+#ifndef BORNAGAIN_CORE_LEGACY_MATRIXRTCOEFFICIENTS_V2_H
+#define BORNAGAIN_CORE_LEGACY_MATRIXRTCOEFFICIENTS_V2_H
#include "Base/Vector/Vectors3D.h"
#include "Sample/RT/ILayerRTCoefficients.h"
@@ -25,7 +25,8 @@
class MatrixRTCoefficients_v2 : public ILayerRTCoefficients {
public:
- friend class SpecularMagneticStrategy;
+ friend class SpecularMagneticStrategy_v2;
+ friend class SpecularMagneticOriginalStrategy;
MatrixRTCoefficients_v2(double kz_sign, Eigen::Vector2cd eigenvalues, kvector_t b);
MatrixRTCoefficients_v2(const MatrixRTCoefficients_v2& other);
@@ -66,4 +67,4 @@ private:
//!< the second eigenmode
};
-#endif // BORNAGAIN_SAMPLE_RT_MATRIXRTCOEFFICIENTS_V2_H
+#endif // BORNAGAIN_CORE_LEGACY_MATRIXRTCOEFFICIENTS_V2_H
diff --git a/Sample/Specular/SpecularMagneticOldStrategy.cpp b/Core/Legacy/SpecularMagneticStrategy_v1.cpp
similarity index 89%
rename from Sample/Specular/SpecularMagneticOldStrategy.cpp
rename to Core/Legacy/SpecularMagneticStrategy_v1.cpp
index 601711f3c4f788dfa5f7ecf4fd5f2ca662edd824..0171d494e2be2c5d7b975a08aaf21f2e530155ac 100644
--- a/Sample/Specular/SpecularMagneticOldStrategy.cpp
+++ b/Core/Legacy/SpecularMagneticStrategy_v1.cpp
@@ -2,8 +2,8 @@
//
// BornAgain: simulate and fit scattering at grazing incidence
//
-//! @file Sample/Specular/SpecularMagneticOldStrategy.cpp
-//! @brief Implements class SpecularMagneticOldStrategy.
+//! @file Sample/Specular/SpecularMagneticStrategy_v1.cpp
+//! @brief Implements class SpecularMagneticStrategy_v1.
//!
//! @homepage http://www.bornagainproject.org
//! @license GNU General Public License v3 or higher (see COPYING)
@@ -12,7 +12,7 @@
//
// ************************************************************************************************
-#include "Sample/Specular/SpecularMagneticOldStrategy.h"
+#include "SpecularMagneticStrategy_v1.h"
#include "Sample/Material/WavevectorInfo.h"
#include "Sample/Multilayer/Layer.h"
#include "Sample/Multilayer/MultiLayer.h"
@@ -22,35 +22,35 @@
namespace {
void CalculateEigenvalues(const std::vector<Slice>& slices, const kvector_t k,
- std::vector<MatrixRTCoefficients>& coeff);
+ std::vector<MatrixRTCoefficients_v1>& coeff);
void CalculateTransferAndBoundary(const std::vector<Slice>& slices,
- std::vector<MatrixRTCoefficients>& coeff);
-void SetForNoTransmission(std::vector<MatrixRTCoefficients>& coeff);
+ std::vector<MatrixRTCoefficients_v1>& coeff);
+void SetForNoTransmission(std::vector<MatrixRTCoefficients_v1>& coeff);
complex_t GetImExponential(complex_t exponent);
} // namespace
-ISpecularStrategy::coeffs_t SpecularMagneticOldStrategy::Execute(const std::vector<Slice>& slices,
+ISpecularStrategy::coeffs_t SpecularMagneticStrategy_v1::Execute(const std::vector<Slice>& slices,
const kvector_t& k) const {
- std::vector<MatrixRTCoefficients> result(slices.size());
+ std::vector<MatrixRTCoefficients_v1> result(slices.size());
CalculateEigenvalues(slices, k, result);
CalculateTransferAndBoundary(slices, result);
- coeffs_t resultConvert;
+ ISpecularStrategy::coeffs_t resultConvert;
for (auto& coeff : result)
- resultConvert.push_back(std::make_unique<MatrixRTCoefficients>(coeff));
+ resultConvert.push_back(std::make_unique<MatrixRTCoefficients_v1>(coeff));
return resultConvert;
}
ISpecularStrategy::coeffs_t
-SpecularMagneticOldStrategy::Execute(const std::vector<Slice>&,
+SpecularMagneticStrategy_v1::Execute(const std::vector<Slice>&,
const std::vector<complex_t>&) const {
throw std::runtime_error("Not implemented");
}
namespace {
void CalculateEigenvalues(const std::vector<Slice>& slices, const kvector_t k,
- std::vector<MatrixRTCoefficients>& coeff) {
+ std::vector<MatrixRTCoefficients_v1>& coeff) {
double mag_k = k.mag();
double n_ref = slices[0].material().refractiveIndex(2 * M_PI / mag_k).real();
double sign_kz = k.z() > 0.0 ? -1.0 : 1.0;
@@ -78,7 +78,7 @@ void CalculateEigenvalues(const std::vector<Slice>& slices, const kvector_t k,
// todo: avoid overflows (see SpecularMatrix.cpp)
void CalculateTransferAndBoundary(const std::vector<Slice>& slices,
- std::vector<MatrixRTCoefficients>& coeff) {
+ std::vector<MatrixRTCoefficients_v1>& coeff) {
size_t N = coeff.size();
if (coeff[0].lambda == Eigen::Vector2cd::Zero() && N > 1) {
SetForNoTransmission(coeff);
@@ -136,7 +136,7 @@ void CalculateTransferAndBoundary(const std::vector<Slice>& slices,
}
}
-void SetForNoTransmission(std::vector<MatrixRTCoefficients>& coeff) {
+void SetForNoTransmission(std::vector<MatrixRTCoefficients_v1>& coeff) {
size_t N = coeff.size();
for (size_t i = 0; i < N; ++i) {
coeff[i].phi_psi_plus.setZero();
diff --git a/Core/Legacy/SpecularMagneticStrategy_v1.h b/Core/Legacy/SpecularMagneticStrategy_v1.h
new file mode 100644
index 0000000000000000000000000000000000000000..ea99ca97e451e62fae64b4558bc576ea326ddc70
--- /dev/null
+++ b/Core/Legacy/SpecularMagneticStrategy_v1.h
@@ -0,0 +1,47 @@
+// ************************************************************************************************
+//
+// BornAgain: simulate and fit scattering at grazing incidence
+//
+//! @file Sample/Specular/SpecularMagneticStrategy_v1.h
+//! @brief Defines class SpecularMagneticStrategy_v1.
+//!
+//! @homepage http://www.bornagainproject.org
+//! @license GNU General Public License v3 or higher (see COPYING)
+//! @copyright Forschungszentrum Jülich GmbH 2018
+//! @authors Scientific Computing Group at MLZ (see CITATION, AUTHORS)
+//
+// ************************************************************************************************
+
+#ifndef BORNAGAIN_CORE_LEGACY_SPECULARMAGNETICSTRATEGY_V1_H
+#define BORNAGAIN_CORE_LEGACY_SPECULARMAGNETICSTRATEGY_V1_H
+
+#include "MatrixRTCoefficients_v1.h"
+#include "Sample/Specular/ISpecularStrategy.h"
+#include <memory>
+#include <vector>
+
+class Slice;
+
+//! Implements the matrix formalism for the calculation of wave amplitudes of
+//! the coherent wave solution in a multilayer with magnetization.
+//! @ingroup algorithms_internal
+
+class SpecularMagneticStrategy_v1 : public ISpecularStrategy {
+public:
+ // TODO remove once external test code is not needed anmyore
+ // for the moment i need them!
+ using coefficient_type = MatrixRTCoefficients_v1;
+ using coefficient_pointer_type = std::unique_ptr<const coefficient_type>;
+ using coeffs_t = std::vector<coefficient_pointer_type>;
+
+ //! Computes refraction angle reflection/transmission coefficients
+ //! for given sliced multilayer and wavevector k
+ ISpecularStrategy::coeffs_t Execute(const std::vector<Slice>& slices,
+ const kvector_t& k) const;
+
+ ISpecularStrategy::coeffs_t Execute(const std::vector<Slice>& slices,
+ const std::vector<complex_t>& kz) const;
+
+}; // class SpecularMagneticStrategy_v1
+
+#endif // BORNAGAIN_CORE_LEGACY_SPECULARMAGNETICSTRATEGY_V1_H
diff --git a/Core/Legacy/SpecularMagneticStrategy_v2.cpp b/Core/Legacy/SpecularMagneticStrategy_v2.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..58cc59474368d4f9095bd697be1d4b94c2181182
--- /dev/null
+++ b/Core/Legacy/SpecularMagneticStrategy_v2.cpp
@@ -0,0 +1,273 @@
+// ************************************************************************************************
+//
+// BornAgain: simulate and fit scattering at grazing incidence
+//
+//! @file Sample/Specular/SpecularMagneticStrategy_v2.cpp
+//! @brief Implements class SpecularMagneticStrategy_v2.
+//!
+//! @homepage http://www.bornagainproject.org
+//! @license GNU General Public License v3 or higher (see COPYING)
+//! @copyright Forschungszentrum Jülich GmbH 2018
+//! @authors Scientific Computing Group at MLZ (see CITATION, AUTHORS)
+//
+// ************************************************************************************************
+
+#include "SpecularMagneticStrategy_v2.h"
+#include "Base/Const/PhysicalConstants.h"
+#include "Sample/Slice/KzComputation.h"
+#include "Sample/Slice/Slice.h"
+
+namespace {
+kvector_t magneticImpact(kvector_t B_field);
+Eigen::Vector2cd eigenvalues(complex_t kz, double b_mag);
+Eigen::Vector2cd checkForUnderflow(const Eigen::Vector2cd& eigenvs);
+complex_t GetImExponential(complex_t exponent);
+
+// The factor 1e-18 is here to have unit: 1/T*nm^-2
+constexpr double magnetic_prefactor = PhysConsts::m_n * PhysConsts::g_factor_n * PhysConsts::mu_N
+ / PhysConsts::h_bar / PhysConsts::h_bar * 1e-18;
+} // namespace
+
+ISpecularStrategy::coeffs_t SpecularMagneticStrategy_v2::Execute(const std::vector<Slice>& slices,
+ const kvector_t& k) const {
+ return Execute(slices, KzComputation::computeReducedKz(slices, k));
+}
+
+ISpecularStrategy::coeffs_t
+SpecularMagneticStrategy_v2::Execute(const std::vector<Slice>& slices,
+ const std::vector<complex_t>& kz) const {
+ if (slices.size() != kz.size())
+ throw std::runtime_error("Number of slices does not match the size of the kz-vector");
+
+ ISpecularStrategy::coeffs_t result;
+ for (auto& coeff : computeTR(slices, kz))
+ result.push_back(std::make_unique<MatrixRTCoefficients_v2>(coeff));
+
+ return result;
+}
+
+std::vector<MatrixRTCoefficients_v2>
+SpecularMagneticStrategy_v2::computeTR(const std::vector<Slice>& slices,
+ const std::vector<complex_t>& kzs) {
+ if (kzs[0] == 0.)
+ throw std::runtime_error("Edge case k_z = 0 not implemented");
+
+ if (slices.size() != kzs.size())
+ throw std::runtime_error(
+ "Error in SpecularMagnetic_::execute: kz vector and slices size shall coinside.");
+ if (slices.empty())
+ return {};
+
+ std::vector<MatrixRTCoefficients_v2> result;
+ result.reserve(slices.size());
+
+ const double kz_sign = kzs.front().real() > 0.0 ? 1.0 : -1.0; // save sign to restore it later
+ const kvector_t b_0 = magneticImpact(slices.front().bField());
+ result.emplace_back(kz_sign, eigenvalues(kzs.front(), 0.0), kvector_t{0.0, 0.0, 0.0});
+ for (size_t i = 1, size = slices.size(); i < size; ++i) {
+ kvector_t b = magneticImpact(slices[i].bField()) - b_0;
+ result.emplace_back(kz_sign, checkForUnderflow(eigenvalues(kzs[i], b.mag())), b);
+ }
+
+ if (result.front().m_lambda == Eigen::Vector2cd::Zero()) {
+ std::for_each(result.begin(), result.end(), [](auto& coeff) { setNoTransmission(coeff); });
+ return result;
+ }
+
+ std::for_each(result.begin(), result.end(), [](auto& coeff) { calculateTR(coeff); });
+ nullifyBottomReflection(result.back());
+ propagateBackwardsForwards(result, slices);
+
+ return result;
+}
+
+void SpecularMagneticStrategy_v2::calculateTR(MatrixRTCoefficients_v2& coeff) {
+ const double b = coeff.m_b.mag();
+ if (b == 0.0) {
+ calculateZeroFieldTR(coeff);
+ return;
+ }
+
+ const double bpbz = b + coeff.m_b.z();
+ const double bmbz = b - coeff.m_b.z();
+ const complex_t bxmby = coeff.m_b.x() - I * coeff.m_b.y();
+ const complex_t bxpby = coeff.m_b.x() + I * coeff.m_b.y();
+ const complex_t l_1 = coeff.m_lambda(0);
+ const complex_t l_2 = coeff.m_lambda(1);
+
+ auto& T1 = coeff.T1;
+ T1 << bmbz, -bxmby, -bmbz * l_1, bxmby * l_1, -bxpby, bpbz, bxpby * l_1, -bpbz * l_1,
+ -bmbz / l_1, bxmby / l_1, bmbz, -bxmby, bxpby / l_1, -bpbz / l_1, -bxpby, bpbz;
+ T1 /= 4.0 * b;
+
+ auto& R1 = coeff.R1;
+ R1 << T1(0, 0), T1(0, 1), -T1(0, 2), -T1(0, 3), T1(1, 0), T1(1, 1), -T1(1, 2), -T1(1, 3),
+ -T1(2, 0), -T1(2, 1), T1(2, 2), T1(2, 3), -T1(3, 0), -T1(3, 1), T1(3, 2), T1(3, 3);
+
+ auto& T2 = coeff.T2;
+ T2 << bpbz, bxmby, -bpbz * l_2, -bxmby * l_2, bxpby, bmbz, -bxpby * l_2, -bmbz * l_2,
+ -bpbz / l_2, -bxmby / l_2, bpbz, bxmby, -bxpby / l_2, -bmbz / l_2, bxpby, bmbz;
+ T2 /= 4.0 * b;
+
+ auto& R2 = coeff.R2;
+ R2 << T2(0, 0), T2(0, 1), -T2(0, 2), -T2(0, 3), T2(1, 0), T2(1, 1), -T2(1, 2), -T2(1, 3),
+ -T2(2, 0), -T2(2, 1), T2(2, 2), T2(2, 3), -T2(3, 0), -T2(3, 1), T2(3, 2), T2(3, 3);
+}
+
+void SpecularMagneticStrategy_v2::calculateZeroFieldTR(MatrixRTCoefficients_v2& coeff) {
+ coeff.T1 = Eigen::Matrix4cd::Zero();
+ coeff.R1 = Eigen::Matrix4cd::Zero();
+ coeff.T2 = Eigen::Matrix4cd::Zero();
+ coeff.R2 = Eigen::Matrix4cd::Zero();
+
+ // lambda_1 == lambda_2, no difference which one to use
+ const complex_t eigen_value = coeff.m_lambda(0);
+
+ Eigen::Matrix3cd Tblock;
+ Tblock << 0.5, 0.0, -0.5 * eigen_value, 0.0, 0.0, 0.0, -0.5 / eigen_value, 0.0, 0.5;
+
+ Eigen::Matrix3cd Rblock;
+ Rblock << 0.5, 0.0, 0.5 * eigen_value, 0.0, 0.0, 0.0, 0.5 / eigen_value, 0.0, 0.5;
+
+ coeff.T1.block<3, 3>(1, 1) = Tblock;
+ coeff.R1.block<3, 3>(1, 1) = Rblock;
+ coeff.T2.block<3, 3>(0, 0) = Tblock;
+ coeff.R2.block<3, 3>(0, 0) = Rblock;
+}
+
+void SpecularMagneticStrategy_v2::setNoTransmission(MatrixRTCoefficients_v2& coeff) {
+ coeff.m_w_plus = Eigen::Vector4cd::Zero();
+ coeff.m_w_min = Eigen::Vector4cd::Zero();
+ coeff.T1 = Eigen::Matrix4cd::Identity() / 4.0;
+ coeff.R1 = coeff.T1;
+ coeff.T2 = coeff.T1;
+ coeff.R2 = coeff.T1;
+}
+
+void SpecularMagneticStrategy_v2::nullifyBottomReflection(MatrixRTCoefficients_v2& coeff) {
+ const complex_t l_1 = coeff.m_lambda(0);
+ const complex_t l_2 = coeff.m_lambda(1);
+ const double b_mag = coeff.m_b.mag();
+ const kvector_t& b = coeff.m_b;
+
+ if (b_mag == 0.0) {
+ // both eigenvalues are the same, no difference which one to take
+ coeff.m_w_plus << -l_1, 0.0, 1.0, 0.0;
+ coeff.m_w_min << 0.0, -l_1, 0.0, 1.0;
+ return;
+ }
+
+ // First basis vector that has no upward going wave amplitude
+ coeff.m_w_min(0) = (b.x() - I * b.y()) * (l_1 - l_2) / 2.0 / b_mag;
+ coeff.m_w_min(1) = b.z() * (l_2 - l_1) / 2.0 / b_mag - (l_1 + l_2) / 2.0;
+ coeff.m_w_min(2) = 0.0;
+ coeff.m_w_min(3) = 1.0;
+
+ // Second basis vector that has no upward going wave amplitude
+ coeff.m_w_plus(0) = -(l_1 + l_2) / 2.0 - b.z() / (l_1 + l_2);
+ coeff.m_w_plus(1) = (b.x() + I * b.y()) * (l_1 - l_2) / 2.0 / b_mag;
+ coeff.m_w_plus(2) = 1.0;
+ coeff.m_w_plus(3) = 0.0;
+}
+
+void SpecularMagneticStrategy_v2::propagateBackwardsForwards(
+ std::vector<MatrixRTCoefficients_v2>& coeff, const std::vector<Slice>& slices) {
+ const int size = static_cast<int>(coeff.size());
+ std::vector<Eigen::Matrix2cd> SMatrices(coeff.size());
+ std::vector<complex_t> Normalization(coeff.size());
+
+ for (int index = size - 2; index >= 0; --index) {
+ const size_t i = static_cast<size_t>(index);
+ const double t = slices[i].thickness();
+ const auto kz = coeff[i].getKz();
+ const Eigen::Matrix4cd l = coeff[i].R1 * GetImExponential(kz(0) * t)
+ + coeff[i].T1 * GetImExponential(-kz(0) * t)
+ + coeff[i].R2 * GetImExponential(kz(1) * t)
+ + coeff[i].T2 * GetImExponential(-kz(1) * t);
+ coeff[i].m_w_plus = l * coeff[i + 1].m_w_plus;
+ coeff[i].m_w_min = l * coeff[i + 1].m_w_min;
+
+ // rotate and normalize polarization
+ const auto Snorm = findNormalizationCoefficients(coeff[i]);
+ auto S = Snorm.first;
+ auto norm = Snorm.second;
+
+ SMatrices[i] = S;
+ Normalization[i] = norm;
+
+ const complex_t a_plus = S(0, 0) / norm;
+ const complex_t b_plus = S(1, 0) / norm;
+ const complex_t a_min = S(0, 1) / norm;
+ const complex_t b_min = S(1, 1) / norm;
+
+ const Eigen::Vector4cd w_plus = a_plus * coeff[i].m_w_plus + b_plus * coeff[i].m_w_min;
+ const Eigen::Vector4cd w_min = a_min * coeff[i].m_w_plus + b_min * coeff[i].m_w_min;
+
+ coeff[i].m_w_plus = std::move(w_plus);
+ coeff[i].m_w_min = std::move(w_min);
+ }
+
+ complex_t dumpingFactor = 1;
+ Eigen::Matrix2cd S = Eigen::Matrix2cd::Identity();
+ for (size_t i = 1; i < coeff.size(); ++i) {
+ dumpingFactor = dumpingFactor * Normalization[i - 1];
+ S = SMatrices[i - 1] * S;
+
+ if (std::isinf(std::norm(dumpingFactor))) {
+ // not entirely sure, whether this is the correct edge case
+ std::for_each(coeff.begin() + i, coeff.end(),
+ [](auto& coeff) { setNoTransmission(coeff); });
+ break;
+ }
+
+ const complex_t a_plus = S(0, 0) / dumpingFactor;
+ const complex_t b_plus = S(1, 0) / dumpingFactor;
+ const complex_t a_min = S(0, 1) / dumpingFactor;
+ const complex_t b_min = S(1, 1) / dumpingFactor;
+
+ Eigen::Vector4cd w_plus = a_plus * coeff[i].m_w_plus + b_plus * coeff[i].m_w_min;
+ Eigen::Vector4cd w_min = a_min * coeff[i].m_w_plus + b_min * coeff[i].m_w_min;
+
+ coeff[i].m_w_plus = std::move(w_plus);
+ coeff[i].m_w_min = std::move(w_min);
+ }
+}
+
+std::pair<Eigen::Matrix2cd, complex_t>
+SpecularMagneticStrategy_v2::findNormalizationCoefficients(const MatrixRTCoefficients_v2& coeff) {
+ const Eigen::Vector2cd Ta = coeff.T1plus() + coeff.T2plus();
+ const Eigen::Vector2cd Tb = coeff.T1min() + coeff.T2min();
+
+ Eigen::Matrix2cd S;
+ S << Ta(0), Tb(0), Ta(1), Tb(1);
+
+ Eigen::Matrix2cd SInverse;
+ SInverse << S(1, 1), -S(0, 1), -S(1, 0), S(0, 0);
+ const complex_t d1 = S(1, 1) - S(0, 1);
+ const complex_t d2 = S(1, 0) - S(0, 0);
+ const complex_t denominator = S(0, 0) * d1 - d2 * S(0, 1);
+
+ return {SInverse, denominator};
+}
+
+namespace {
+kvector_t magneticImpact(kvector_t B_field) {
+ return -magnetic_prefactor * B_field;
+}
+
+Eigen::Vector2cd eigenvalues(complex_t kz, double b_mag) {
+ const complex_t a = kz * kz;
+ return {I * std::sqrt(a + b_mag), I * std::sqrt(a - b_mag)};
+}
+
+Eigen::Vector2cd checkForUnderflow(const Eigen::Vector2cd& eigenvs) {
+ auto lambda = [](complex_t value) { return std::abs(value) < 1e-40 ? 1e-40 : value; };
+ return {lambda(eigenvs(0)), lambda(eigenvs(1))};
+}
+
+complex_t GetImExponential(complex_t exponent) {
+ if (exponent.imag() > -std::log(std::numeric_limits<double>::min()))
+ return 0.0;
+ return std::exp(I * exponent);
+}
+} // namespace
diff --git a/Core/Legacy/SpecularMagneticStrategy_v2.h b/Core/Legacy/SpecularMagneticStrategy_v2.h
new file mode 100644
index 0000000000000000000000000000000000000000..1de56608eabc4df319e9f01513d2c2bd5234e287
--- /dev/null
+++ b/Core/Legacy/SpecularMagneticStrategy_v2.h
@@ -0,0 +1,77 @@
+// ************************************************************************************************
+//
+// BornAgain: simulate and fit scattering at grazing incidence
+//
+//! @file Sample/Specular/SpecularMagneticStrategy_v2.h
+//! @brief Defines class SpecularMagneticStrategy_v2.
+//!
+//! @homepage http://www.bornagainproject.org
+//! @license GNU General Public License v3 or higher (see COPYING)
+//! @copyright Forschungszentrum Jülich GmbH 2018
+//! @authors Scientific Computing Group at MLZ (see CITATION, AUTHORS)
+//
+// ************************************************************************************************
+
+#ifndef BORNAGAIN_CORE_LEGACY_SPECULARMAGNETICSTRATEGY_V2_H
+#define BORNAGAIN_CORE_LEGACY_SPECULARMAGNETICSTRATEGY_V2_H
+
+#include "MatrixRTCoefficients_v2.h"
+#include "Sample/Specular/ISpecularStrategy.h"
+#include <memory>
+#include <vector>
+
+class Slice;
+
+//! Implements the magnetic Fresnel computation without roughness
+//!
+//! Implements the matrix formalism for the calculation of wave amplitudes of
+//! the coherent wave solution in a multilayer with magnetization.
+//! For a detailed description see internal document "Polarized Specular Reflectometry"
+//!
+//! @ingroup algorithms_internal
+class SpecularMagneticStrategy_v2 : public ISpecularStrategy {
+public:
+ // TODO remove once external test code is not needed anmyore
+ // for the moment i need them!
+ using coefficient_type = MatrixRTCoefficients_v2;
+ using coefficient_pointer_type = std::unique_ptr<const coefficient_type>;
+ using coeffs_t = std::vector<coefficient_pointer_type>;
+
+ //! Computes refraction angle reflection/transmission coefficients
+ //! for given sliced multilayer and wavevector k
+ ISpecularStrategy::coeffs_t Execute(const std::vector<Slice>& slices,
+ const kvector_t& k) const;
+
+ //! Computes refraction angle reflection/transmission coefficients
+ //! for given sliced multilayer and a set of kz projections corresponding to each slice
+ ISpecularStrategy::coeffs_t Execute(const std::vector<Slice>& slices,
+ const std::vector<complex_t>& kz) const;
+
+private:
+ static std::vector<MatrixRTCoefficients_v2> computeTR(const std::vector<Slice>& slices,
+ const std::vector<complex_t>& kzs);
+
+ //! Computes frobenius matrices for multilayer solution
+ static void calculateTR(MatrixRTCoefficients_v2& coeff);
+ static void calculateZeroFieldTR(MatrixRTCoefficients_v2& coeff);
+
+ static void setNoTransmission(MatrixRTCoefficients_v2& coeff);
+
+ //! initializes reflectionless bottom boundary condition.
+ static void nullifyBottomReflection(MatrixRTCoefficients_v2& coeff);
+
+ //! Propagates boundary conditions from the bottom to the top of the layer stack.
+ //! Used to compute boundary conditions from the bottom one (with nullified reflection)
+ //! simultaneously propagates amplitudes forward again
+ //! Due to the use of temporary objects this is combined into one function now
+ static void propagateBackwardsForwards(std::vector<MatrixRTCoefficients_v2>& coeff,
+ const std::vector<Slice>& slices);
+
+ //! finds linear coefficients for normalizing transmitted wave to unity.
+ //! The left column of the returned matrix corresponds to the coefficients for pure spin-up
+ //! wave, while the right column - to the coefficients for the spin-down one.
+ static std::pair<Eigen::Matrix2cd, complex_t>
+ findNormalizationCoefficients(const MatrixRTCoefficients_v2& coeff);
+};
+
+#endif // BORNAGAIN_CORE_LEGACY_SPECULARMAGNETICSTRATEGY_V2_H
diff --git a/Examples/Python/sim01_Particles/CylindersAndPrisms.py b/Examples/Python/sim01_Particles/CylindersAndPrisms.py
index 12b91c1c085e607d189d54d251f9a0e745e40171..cbab1b3b6d47c381c6f2adc6dc66d69f608dbf28 100644
--- a/Examples/Python/sim01_Particles/CylindersAndPrisms.py
+++ b/Examples/Python/sim01_Particles/CylindersAndPrisms.py
@@ -2,38 +2,49 @@
Mixture of cylinders and prisms without interference
"""
import bornagain as ba
-from bornagain import deg, angstrom, nm
+from bornagain import angstrom, deg, nm, nm2, kvector_t
def get_sample():
"""
Returns a sample with uncorrelated cylinders and prisms on a substrate.
"""
- # defining materials
- m_vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
- m_substrate = ba.HomogeneousMaterial("Substrate", 6e-6, 2e-8)
- m_particle = ba.HomogeneousMaterial("Particle", 6e-4, 2e-8)
-
- # collection of particles
- cylinder_ff = ba.FormFactorCylinder(5*nm, 5*nm)
- cylinder = ba.Particle(m_particle, cylinder_ff)
- prism_ff = ba.FormFactorPrism3(10*nm, 5*nm)
- prism = ba.Particle(m_particle, prism_ff)
- particle_layout = ba.ParticleLayout()
- particle_layout.addParticle(cylinder, 0.5)
- particle_layout.addParticle(prism, 0.5)
- interference = ba.InterferenceFunctionNone()
- particle_layout.setInterferenceFunction(interference)
-
- # vacuum layer with particles and substrate form multi layer
- vacuum_layer = ba.Layer(m_vacuum)
- vacuum_layer.addLayout(particle_layout)
- substrate_layer = ba.Layer(m_substrate)
- multi_layer = ba.MultiLayer()
- multi_layer.addLayer(vacuum_layer)
- multi_layer.addLayer(substrate_layer)
- print(multi_layer.treeToString())
- return multi_layer
+
+ # Define materials
+ material_Particle = ba.HomogeneousMaterial("Particle", 0.0006, 2e-08)
+ material_Substrate = ba.HomogeneousMaterial("Substrate", 6e-06, 2e-08)
+ material_Vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
+
+ # Define form factors
+ ff_1 = ba.FormFactorCylinder(5.0*nm, 5.0*nm)
+ ff_2 = ba.FormFactorPrism3(10.0*nm, 5.0*nm)
+
+ # Define particles
+ particle_1 = ba.Particle(material_Particle, ff_1)
+ particle_2 = ba.Particle(material_Particle, ff_2)
+
+ # Define interference functions
+ iff = ba.InterferenceFunctionNone()
+
+ # Define particle layouts
+ layout = ba.ParticleLayout()
+ layout.addParticle(particle_1, 0.5)
+ layout.addParticle(particle_2, 0.5)
+ layout.setInterferenceFunction(iff)
+ layout.setWeight(1)
+ layout.setTotalParticleSurfaceDensity(0.01)
+
+ # Define layers
+ layer_1 = ba.Layer(material_Vacuum)
+ layer_1.addLayout(layout)
+ layer_2 = ba.Layer(material_Substrate)
+
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer_1)
+ sample.addLayer(layer_2)
+
+ return sample
def get_simulation():
diff --git a/Examples/Python/sim01_Particles/CylindersInBA.py b/Examples/Python/sim01_Particles/CylindersInBA.py
index 69e2e08e03e4859478dd1189c4d2b0c96a0dc3ee..26855e44a93ba0f63acf6bb3e266b453afc3050d 100644
--- a/Examples/Python/sim01_Particles/CylindersInBA.py
+++ b/Examples/Python/sim01_Particles/CylindersInBA.py
@@ -2,7 +2,7 @@
Cylinder form factor in Born approximation
"""
import bornagain as ba
-from bornagain import deg, angstrom, nm
+from bornagain import angstrom, deg, nm, nm2, kvector_t
def get_sample():
@@ -10,22 +10,32 @@ def get_sample():
Returns a sample with cylinders in a homogeneous environment ("Vacuum"),
implying a simulation in plain Born approximation.
"""
- # defining materials
- m_vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
- m_particle = ba.HomogeneousMaterial("Particle", 6e-4, 2e-8)
-
- # collection of particles
- cylinder_ff = ba.FormFactorCylinder(5*nm, 5*nm)
- cylinder = ba.Particle(m_particle, cylinder_ff)
- particle_layout = ba.ParticleLayout()
- particle_layout.addParticle(cylinder, 1.0)
-
- vacuum_layer = ba.Layer(m_vacuum)
- vacuum_layer.addLayout(particle_layout)
-
- multi_layer = ba.MultiLayer()
- multi_layer.addLayer(vacuum_layer)
- return multi_layer
+
+ # Define materials
+ material_Particle = ba.HomogeneousMaterial("Particle", 0.0006, 2e-08)
+ material_Vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
+
+ # Define form factors
+ ff = ba.FormFactorCylinder(5.0*nm, 5.0*nm)
+
+ # Define particles
+ particle = ba.Particle(material_Particle, ff)
+
+ # Define particle layouts
+ layout = ba.ParticleLayout()
+ layout.addParticle(particle, 1.0)
+ layout.setWeight(1)
+ layout.setTotalParticleSurfaceDensity(0.01)
+
+ # Define layers
+ layer = ba.Layer(material_Vacuum)
+ layer.addLayout(layout)
+
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer)
+
+ return sample
def get_simulation():
diff --git a/Examples/Python/sim01_Particles/CylindersInDWBA.py b/Examples/Python/sim01_Particles/CylindersInDWBA.py
index 2d60be5dd1a91c96c073bec41a57a67a417c95fa..28d884bf62df1e399244fdd48aba86bc76b148d0 100644
--- a/Examples/Python/sim01_Particles/CylindersInDWBA.py
+++ b/Examples/Python/sim01_Particles/CylindersInDWBA.py
@@ -2,32 +2,42 @@
Cylinder form factor in DWBA
"""
import bornagain as ba
-from bornagain import deg, angstrom, nm
+from bornagain import angstrom, deg, nm, nm2, kvector_t
def get_sample():
"""
Returns a sample with cylinders on a substrate.
"""
- # defining materials
- m_vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
- m_substrate = ba.HomogeneousMaterial("Substrate", 6e-6, 2e-8)
- m_particle = ba.HomogeneousMaterial("Particle", 6e-4, 2e-8)
-
- # collection of particles
- cylinder_ff = ba.FormFactorCylinder(5*nm, 5*nm)
- cylinder = ba.Particle(m_particle, cylinder_ff)
- particle_layout = ba.ParticleLayout()
- particle_layout.addParticle(cylinder, 1.0)
-
- vacuum_layer = ba.Layer(m_vacuum)
- vacuum_layer.addLayout(particle_layout)
- substrate_layer = ba.Layer(m_substrate)
-
- multi_layer = ba.MultiLayer()
- multi_layer.addLayer(vacuum_layer)
- multi_layer.addLayer(substrate_layer)
- return multi_layer
+
+ # Define materials
+ material_Particle = ba.HomogeneousMaterial("Particle", 0.0006, 2e-08)
+ material_Substrate = ba.HomogeneousMaterial("Substrate", 6e-06, 2e-08)
+ material_Vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
+
+ # Define form factors
+ ff = ba.FormFactorCylinder(5.0*nm, 5.0*nm)
+
+ # Define particles
+ particle = ba.Particle(material_Particle, ff)
+
+ # Define particle layouts
+ layout = ba.ParticleLayout()
+ layout.addParticle(particle, 1.0)
+ layout.setWeight(1)
+ layout.setTotalParticleSurfaceDensity(0.01)
+
+ # Define layers
+ layer_1 = ba.Layer(material_Vacuum)
+ layer_1.addLayout(layout)
+ layer_2 = ba.Layer(material_Substrate)
+
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer_1)
+ sample.addLayer(layer_2)
+
+ return sample
def get_simulation():
diff --git a/Examples/Python/sim01_Particles/CylindersWithSizeDistribution.py b/Examples/Python/sim01_Particles/CylindersWithSizeDistribution.py
index a5dfd5d20c4d7568dee0a3cc38741b78b6bf6066..adc1fd1decbcae25800c687551080346903e08d2 100644
--- a/Examples/Python/sim01_Particles/CylindersWithSizeDistribution.py
+++ b/Examples/Python/sim01_Particles/CylindersWithSizeDistribution.py
@@ -2,7 +2,7 @@
Cylinders with size distribution
"""
import bornagain as ba
-from bornagain import deg, angstrom, nm
+from bornagain import angstrom, deg, nm, nm2, kvector_t
def get_sample():
@@ -10,38 +10,38 @@ def get_sample():
Return a sample with cylinders on a substrate.
The cylinders have a Gaussian size distribution.
"""
- m_vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
- m_particle = ba.HomogeneousMaterial("Particle", 6e-4, 2e-8)
-
- # cylindrical particle
- radius = 5*nm
- height = radius
- cylinder_ff = ba.FormFactorCylinder(radius, height)
- cylinder = ba.Particle(m_particle, cylinder_ff)
-
- # collection of particles with size distribution
- nparticles = 100
- sigma = 0.2*radius
-
- gauss_distr = ba.DistributionGaussian(radius, sigma)
-
- sigma_factor = 2.0
- par_distr = ba.ParameterDistribution("/Particle/Cylinder/Radius",
- gauss_distr, nparticles, sigma_factor)
- # by uncommenting the line below, the height of the cylinders
- # can be scaled proportionally to the radius:
- # par_distr.linkParameter("/Particle/Cylinder/Height")
- part_coll = ba.ParticleDistribution(cylinder, par_distr)
-
- # assembling the sample
- particle_layout = ba.ParticleLayout()
- particle_layout.addParticle(part_coll)
-
- vacuum_layer = ba.Layer(m_vacuum)
- vacuum_layer.addLayout(particle_layout)
- multi_layer = ba.MultiLayer()
- multi_layer.addLayer(vacuum_layer)
- return multi_layer
+
+ # Define materials
+ material_Particle = ba.HomogeneousMaterial("Particle", 0.0006, 2e-08)
+ material_Vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
+
+ # Define form factors
+ ff = ba.FormFactorCylinder(5.0*nm, 5.0*nm)
+
+ # Define particles
+ particle = ba.Particle(material_Particle, ff)
+
+ # Define particles with parameter following a distribution
+ distr_1 = ba.DistributionGaussian(5.0*nm, 1.0*nm)
+ par_distr_1 = ba.ParameterDistribution("/Particle/Cylinder/Radius", distr_1,
+ 100, 2.0)
+ particle_distrib = ba.ParticleDistribution(particle, par_distr_1)
+
+ # Define particle layouts
+ layout = ba.ParticleLayout()
+ layout.addParticle(particle_distrib, 1.0)
+ layout.setWeight(1)
+ layout.setTotalParticleSurfaceDensity(0.01)
+
+ # Define layers
+ layer = ba.Layer(material_Vacuum)
+ layer.addLayout(layout)
+
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer)
+
+ return sample
def get_simulation():
diff --git a/Examples/Python/sim01_Particles/RotatedPyramids.py b/Examples/Python/sim01_Particles/RotatedPyramids.py
index 4b8f260c3c685cf9d3e5eae9e7c3b7bd77dc4997..0077becb7bfe0fbb2de1ba8ed717227810b81960 100644
--- a/Examples/Python/sim01_Particles/RotatedPyramids.py
+++ b/Examples/Python/sim01_Particles/RotatedPyramids.py
@@ -2,34 +2,44 @@
Rotated pyramids on top of substrate
"""
import bornagain as ba
-from bornagain import deg, angstrom, nm
+from bornagain import angstrom, deg, nm, nm2, kvector_t
def get_sample():
"""
Returns a sample with rotated pyramids on top of a substrate.
"""
- # defining materials
- m_vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
- m_substrate = ba.HomogeneousMaterial("Substrate", 6e-6, 2e-8)
- m_particle = ba.HomogeneousMaterial("Particle", 6e-4, 2e-8)
-
- # collection of particles
- pyramid_ff = ba.FormFactorPyramid(10*nm, 5*nm, 54.73*deg)
- pyramid = ba.Particle(m_particle, pyramid_ff)
- transform = ba.RotationZ(45.*deg)
- particle_layout = ba.ParticleLayout()
- particle_layout.addParticle(pyramid, 1.0, ba.kvector_t(0.0, 0.0, 0.0),
- transform)
-
- vacuum_layer = ba.Layer(m_vacuum)
- vacuum_layer.addLayout(particle_layout)
- substrate_layer = ba.Layer(m_substrate)
-
- multi_layer = ba.MultiLayer()
- multi_layer.addLayer(vacuum_layer)
- multi_layer.addLayer(substrate_layer)
- return multi_layer
+
+ # Define materials
+ material_Particle = ba.HomogeneousMaterial("Particle", 0.0006, 2e-08)
+ material_Substrate = ba.HomogeneousMaterial("Substrate", 6e-06, 2e-08)
+ material_Vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
+
+ # Define form factors
+ ff = ba.FormFactorPyramid(10.0*nm, 5.0*nm, 54.73*deg)
+
+ # Define particles
+ particle = ba.Particle(material_Particle, ff)
+ particle_rotation = ba.RotationZ(45.0*deg)
+ particle.setRotation(particle_rotation)
+
+ # Define particle layouts
+ layout = ba.ParticleLayout()
+ layout.addParticle(particle, 1.0)
+ layout.setWeight(1)
+ layout.setTotalParticleSurfaceDensity(0.01)
+
+ # Define layers
+ layer_1 = ba.Layer(material_Vacuum)
+ layer_1.addLayout(layout)
+ layer_2 = ba.Layer(material_Substrate)
+
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer_1)
+ sample.addLayer(layer_2)
+
+ return sample
def get_simulation():
diff --git a/Examples/Python/sim01_Particles/TwoTypesOfCylindersWithSizeDistribution.py b/Examples/Python/sim01_Particles/TwoTypesOfCylindersWithSizeDistribution.py
index 064dd0b650f93f2fb5723c3ac4cf960435ccbf9e..4609d50824a4e035e5a8a960b09a00ec1564754b 100644
--- a/Examples/Python/sim01_Particles/TwoTypesOfCylindersWithSizeDistribution.py
+++ b/Examples/Python/sim01_Particles/TwoTypesOfCylindersWithSizeDistribution.py
@@ -2,7 +2,7 @@
Mixture cylinder particles with different size distribution
"""
import bornagain as ba
-from bornagain import deg, angstrom, nm
+from bornagain import angstrom, deg, nm, nm2, kvector_t
def get_sample():
@@ -10,57 +10,47 @@ def get_sample():
Returns a sample with cylinders in a homogeneous medium ("Vacuum").
The cylinders are a 95:5 mixture of two different size distributions.
"""
- # defining materials
- m_vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
- m_particle = ba.HomogeneousMaterial("Particle", 6e-4, 2e-8)
- # collection of particles #1
- radius1 = 5.0*nm
- height1 = radius1
- sigma1 = radius1*0.2
-
- cylinder_ff1 = ba.FormFactorCylinder(radius1, height1)
- cylinder1 = ba.Particle(m_particle, cylinder_ff1)
-
- gauss_distr1 = ba.DistributionGaussian(radius1, sigma1)
-
- nparticles = 150
- sigma_factor = 3.0
-
- # limits will assure, that generated Radius'es are >=0
- limits = ba.RealLimits.nonnegative()
-
- par_distr1 = ba.ParameterDistribution("/Particle/Cylinder/Radius",
- gauss_distr1, nparticles,
- sigma_factor, limits)
- part_coll1 = ba.ParticleDistribution(cylinder1, par_distr1)
-
- # collection of particles #2
- radius2 = 10.0*nm
- height2 = radius2
- sigma2 = radius2*0.02
-
- cylinder_ff2 = ba.FormFactorCylinder(radius2, height2)
- cylinder2 = ba.Particle(m_particle, cylinder_ff2)
-
- gauss_distr2 = ba.DistributionGaussian(radius2, sigma2)
-
- par_distr2 = ba.ParameterDistribution("/Particle/Cylinder/Radius",
- gauss_distr2, nparticles,
- sigma_factor, limits)
- part_coll2 = ba.ParticleDistribution(cylinder2, par_distr2)
-
- # assembling the sample
- particle_layout = ba.ParticleLayout()
- particle_layout.addParticle(part_coll1, 0.95)
- particle_layout.addParticle(part_coll2, 0.05)
-
- vacuum_layer = ba.Layer(m_vacuum)
- vacuum_layer.addLayout(particle_layout)
-
- multi_layer = ba.MultiLayer()
- multi_layer.addLayer(vacuum_layer)
- return multi_layer
+ # Define materials
+ material_Particle = ba.HomogeneousMaterial("Particle", 0.0006, 2e-08)
+ material_Vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
+
+ # Define form factors
+ ff_1 = ba.FormFactorCylinder(5.0*nm, 5.0*nm)
+ ff_2 = ba.FormFactorCylinder(10.0*nm, 10.0*nm)
+
+ # Define particles
+ particle_1 = ba.Particle(material_Particle, ff_1)
+ particle_2 = ba.Particle(material_Particle, ff_2)
+
+ # Define particles with parameter following a distribution
+ distr_1 = ba.DistributionGaussian(5.0*nm, 1.0*nm)
+ par_distr_1 = ba.ParameterDistribution("/Particle/Cylinder/Radius",
+ distr_1, 150, 3.0,
+ ba.RealLimits.nonnegative())
+ particle_distrib_1 = ba.ParticleDistribution(particle_1, par_distr_1)
+ distr_2 = ba.DistributionGaussian(10.0*nm, 0.2*nm)
+ par_distr_2 = ba.ParameterDistribution("/Particle/Cylinder/Radius",
+ distr_2, 150, 3.0,
+ ba.RealLimits.nonnegative())
+ particle_distrib_2 = ba.ParticleDistribution(particle_2, par_distr_2)
+
+ # Define particle layouts
+ layout = ba.ParticleLayout()
+ layout.addParticle(particle_distrib_1, 0.95)
+ layout.addParticle(particle_distrib_2, 0.05)
+ layout.setWeight(1)
+ layout.setTotalParticleSurfaceDensity(0.01)
+
+ # Define layers
+ layer = ba.Layer(material_Vacuum)
+ layer.addLayout(layout)
+
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer)
+
+ return sample
def get_simulation():
diff --git a/Examples/Python/sim02_Complexes/BiMaterialCylinders.py b/Examples/Python/sim02_Complexes/BiMaterialCylinders.py
index 959735afef08db72e1287c57a3e287d436c8ac55..d2c4b45a9bb0b4c96391021b1a34f63c5fec72c6 100644
--- a/Examples/Python/sim02_Complexes/BiMaterialCylinders.py
+++ b/Examples/Python/sim02_Complexes/BiMaterialCylinders.py
@@ -3,7 +3,7 @@ Cylindrical particle made from two materials.
Particle crosses air/substrate interface.
"""
import bornagain as ba
-from bornagain import deg, angstrom, nm
+from bornagain import angstrom, deg, nm, nm2, kvector_t
def get_composition(top_material,
@@ -35,30 +35,47 @@ def get_sample():
Particle shifted down to cross interface.
"""
- # defining materials
- m_vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
- m_substrate = ba.HomogeneousMaterial("Substrate", 3.212e-6, 3.244e-8)
- m_top_part = ba.HomogeneousMaterial("Ag", 1.245e-5, 5.419e-7)
- m_bottom_part = ba.HomogeneousMaterial("Teflon", 2.900e-6, 6.019e-9)
-
- # collection of particles
- composition = get_composition(m_top_part, m_bottom_part)
- shift = 10.0*nm
- composition.setPosition(0, 0, -shift) # will be shifted below interface
-
- particle_layout = ba.ParticleLayout()
- particle_layout.addParticle(composition)
- particle_layout.setTotalParticleSurfaceDensity(1)
-
- # vacuum layer with particles and substrate form multi layer
- vacuum_layer = ba.Layer(m_vacuum)
- vacuum_layer.addLayout(particle_layout)
- substrate_layer = ba.Layer(m_substrate)
- multi_layer = ba.MultiLayer()
- multi_layer.addLayer(vacuum_layer)
- multi_layer.addLayer(substrate_layer)
- print(multi_layer.treeToString())
- return multi_layer
+ # Define materials
+ material_Ag = ba.HomogeneousMaterial("Ag", 1.245e-05, 5.419e-07)
+ material_Substrate = ba.HomogeneousMaterial("Substrate", 3.212e-06,
+ 3.244e-08)
+ material_Teflon = ba.HomogeneousMaterial("Teflon", 2.9e-06, 6.019e-09)
+ material_Vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
+
+ # Define form factors
+ ff_1 = ba.FormFactorCylinder(10.0*nm, 4.0*nm)
+ ff_2 = ba.FormFactorCylinder(10.0*nm, 10.0*nm)
+
+ # Define particles
+ particle_1 = ba.Particle(material_Ag, ff_1)
+ particle_1_position = kvector_t(0.0*nm, 0.0*nm, 10.0*nm)
+ particle_1.setPosition(particle_1_position)
+ particle_2 = ba.Particle(material_Teflon, ff_2)
+
+ # Define composition of particles at specific positions
+ particle_3 = ba.ParticleComposition()
+ particle_3.addParticle(particle_1)
+ particle_3.addParticle(particle_2)
+ particle_3_position = kvector_t(0.0*nm, 0.0*nm, -10.0*nm)
+ particle_3.setPosition(particle_3_position)
+
+ # Define particle layouts
+ layout = ba.ParticleLayout()
+ layout.addParticle(particle_3, 1.0)
+ layout.setWeight(1)
+ layout.setTotalParticleSurfaceDensity(1)
+
+ # Define layers
+ layer_1 = ba.Layer(material_Vacuum)
+ layer_1.addLayout(layout)
+ layer_2 = ba.Layer(material_Substrate)
+
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer_1)
+ sample.addLayer(layer_2)
+
+ return sample
def get_simulation():
diff --git a/Examples/Python/sim02_Complexes/CoreShellNanoparticles.py b/Examples/Python/sim02_Complexes/CoreShellNanoparticles.py
index 9f26bd7df5b29521b0bdc22b13511e63b0251200..0ff5326cf1f6ba7d591fdfa20063ec2540d3fde2 100644
--- a/Examples/Python/sim02_Complexes/CoreShellNanoparticles.py
+++ b/Examples/Python/sim02_Complexes/CoreShellNanoparticles.py
@@ -2,39 +2,49 @@
Core shell nanoparticles
"""
import bornagain as ba
-from bornagain import deg, angstrom, nm
+from bornagain import angstrom, deg, nm, nm2, kvector_t
def get_sample():
"""
Returns a sample with box-shaped core-shell particles on a substrate.
"""
- # defining materials
- m_vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
- m_shell = ba.HomogeneousMaterial("Shell", 1e-4, 2e-8)
- m_core = ba.HomogeneousMaterial("Core", 6e-5, 2e-8)
-
- # collection of particles
- parallelepiped1_ff = ba.FormFactorBox(16*nm, 16*nm, 8*nm)
- parallelepiped2_ff = ba.FormFactorBox(12*nm, 12*nm, 7*nm)
- shell_particle = ba.Particle(m_shell, parallelepiped1_ff)
- core_particle = ba.Particle(m_core, parallelepiped2_ff)
- core_position = ba.kvector_t(0.0, 0.0, 0.0)
-
- particle = ba.ParticleCoreShell(shell_particle, core_particle,
- core_position)
- particle_layout = ba.ParticleLayout()
- particle_layout.addParticle(particle)
- interference = ba.InterferenceFunctionNone()
- particle_layout.setInterferenceFunction(interference)
-
- vacuum_layer = ba.Layer(m_vacuum)
- vacuum_layer.addLayout(particle_layout)
-
- multi_layer = ba.MultiLayer()
- multi_layer.addLayer(vacuum_layer)
-
- return multi_layer
+
+ # Define materials
+ material_Core = ba.HomogeneousMaterial("Core", 6e-05, 2e-08)
+ material_Shell = ba.HomogeneousMaterial("Shell", 0.0001, 2e-08)
+ material_Vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
+
+ # Define form factors
+ ff_1 = ba.FormFactorBox(12.0*nm, 12.0*nm, 7.0*nm)
+ ff_2 = ba.FormFactorBox(16.0*nm, 16.0*nm, 8.0*nm)
+
+ # Define particles
+ particle_1 = ba.Particle(material_Core, ff_1)
+ particle_2 = ba.Particle(material_Shell, ff_2)
+
+ # Define core shell particles
+ particle_3 = ba.ParticleCoreShell(particle_2, particle_1)
+
+ # Define interference functions
+ iff = ba.InterferenceFunctionNone()
+
+ # Define particle layouts
+ layout = ba.ParticleLayout()
+ layout.addParticle(particle_3, 1.0)
+ layout.setInterferenceFunction(iff)
+ layout.setWeight(1)
+ layout.setTotalParticleSurfaceDensity(0.01)
+
+ # Define layers
+ layer = ba.Layer(material_Vacuum)
+ layer.addLayout(layout)
+
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer)
+
+ return sample
def get_simulation():
diff --git a/Examples/Python/sim02_Complexes/HexagonalLatticesWithBasis.py b/Examples/Python/sim02_Complexes/HexagonalLatticesWithBasis.py
index 951758ed0b6e9422d2bcbbbd4002630943fdb26c..860732d97d96c21cdd8a05304edf2b4bc84b9855 100644
--- a/Examples/Python/sim02_Complexes/HexagonalLatticesWithBasis.py
+++ b/Examples/Python/sim02_Complexes/HexagonalLatticesWithBasis.py
@@ -3,7 +3,7 @@ Spheres on two hexagonal close packed layers
"""
import numpy
import bornagain as ba
-from bornagain import deg, angstrom, nm
+from bornagain import angstrom, deg, nm, nm2, kvector_t
def get_sample():
@@ -11,35 +11,53 @@ def get_sample():
Returns a sample with spheres on a substrate,
forming two hexagonal close packed layers.
"""
- m_vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
- m_substrate = ba.HomogeneousMaterial("Substrate", 6e-6, 2e-8)
- m_particle = ba.HomogeneousMaterial("Particle", 6e-4, 2e-8)
-
- radius = 10.0*nm
- sphere_ff = ba.FormFactorFullSphere(radius)
- sphere = ba.Particle(m_particle, sphere_ff)
- particle_layout = ba.ParticleLayout()
-
- pos0 = ba.kvector_t(0.0, 0.0, 0.0)
- pos1 = ba.kvector_t(radius, radius, numpy.sqrt(3.0)*radius)
- basis = ba.ParticleComposition()
- basis.addParticles(sphere, [pos0, pos1])
- particle_layout.addParticle(basis)
-
- interference = ba.InterferenceFunction2DLattice(
- ba.HexagonalLattice2D(radius*2.0, 0*deg))
- pdf = ba.FTDecayFunction2DCauchy(10*nm, 10*nm, 0)
- interference.setDecayFunction(pdf)
-
- particle_layout.setInterferenceFunction(interference)
-
- vacuum_layer = ba.Layer(m_vacuum)
- vacuum_layer.addLayout(particle_layout)
- substrate_layer = ba.Layer(m_substrate, 0)
- multi_layer = ba.MultiLayer()
- multi_layer.addLayer(vacuum_layer)
- multi_layer.addLayer(substrate_layer)
- return multi_layer
+
+ # Define materials
+ material_Particle = ba.HomogeneousMaterial("Particle", 0.0006, 2e-08)
+ material_Substrate = ba.HomogeneousMaterial("Substrate", 6e-06, 2e-08)
+ material_Vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
+
+ # Define form factors
+ ff_1 = ba.FormFactorFullSphere(10.0*nm)
+ ff_2 = ba.FormFactorFullSphere(10.0*nm)
+
+ # Define particles
+ particle_1 = ba.Particle(material_Particle, ff_1)
+ particle_2 = ba.Particle(material_Particle, ff_2)
+ particle_2_position = kvector_t(10.0*nm, 10.0*nm, 17.3205080757*nm)
+ particle_2.setPosition(particle_2_position)
+
+ # Define composition of particles at specific positions
+ particle_3 = ba.ParticleComposition()
+ particle_3.addParticle(particle_1)
+ particle_3.addParticle(particle_2)
+
+ # Define 2D lattices
+ lattice = ba.BasicLattice2D(20.0*nm, 20.0*nm, 120.0*deg, 0.0*deg)
+
+ # Define interference functions
+ iff = ba.InterferenceFunction2DLattice(lattice)
+ iff_pdf = ba.FTDecayFunction2DCauchy(10.0*nm, 10.0*nm, 0.0*deg)
+ iff.setDecayFunction(iff_pdf)
+
+ # Define particle layouts
+ layout = ba.ParticleLayout()
+ layout.addParticle(particle_3, 1.0)
+ layout.setInterferenceFunction(iff)
+ layout.setWeight(1)
+ layout.setTotalParticleSurfaceDensity(0.00288675134595)
+
+ # Define layers
+ layer_1 = ba.Layer(material_Vacuum)
+ layer_1.addLayout(layout)
+ layer_2 = ba.Layer(material_Substrate)
+
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer_1)
+ sample.addLayer(layer_2)
+
+ return sample
def get_simulation():
diff --git a/Examples/Python/sim02_Complexes/MesoCrystal.py b/Examples/Python/sim02_Complexes/MesoCrystal.py
index 84572be605e5ddaec54f1d863857048727e1c693..3bd203daadd04ab6d752e726b4b8dc616ddf2e31 100644
--- a/Examples/Python/sim02_Complexes/MesoCrystal.py
+++ b/Examples/Python/sim02_Complexes/MesoCrystal.py
@@ -2,46 +2,54 @@
Cylindrical mesocrystal on a substrate
"""
import bornagain as ba
-from bornagain import deg, angstrom, nm
+from bornagain import angstrom, deg, nm, nm2, kvector_t
def get_sample():
"""
Returns a sample with a cylindrically shaped mesocrystal on a substrate.
"""
- # defining materials
- m_vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
- m_substrate = ba.HomogeneousMaterial("Substrate", 6e-6, 2e-8)
- m_particle = ba.HomogeneousMaterial("Particle", 6e-4, 2e-8)
-
- # mesocrystal lattice
- lattice_basis_1 = ba.kvector_t(5.0, 0.0, 0.0)
- lattice_basis_2 = ba.kvector_t(0.0, 5.0, 0.0)
- lattice_basis_3 = ba.kvector_t(0.0, 0.0, 5.0)
- lattice = ba.Lattice3D(lattice_basis_1, lattice_basis_2, lattice_basis_3)
-
- # spherical particle that forms the base of the mesocrystal
- sphere_ff = ba.FormFactorFullSphere(2*nm)
- sphere = ba.Particle(m_particle, sphere_ff)
-
- # crystal structure
- crystal = ba.Crystal(sphere, lattice)
-
- # mesocrystal
- meso_ff = ba.FormFactorCylinder(20*nm, 50*nm)
- meso = ba.MesoCrystal(crystal, meso_ff)
-
- particle_layout = ba.ParticleLayout()
- particle_layout.addParticle(meso)
-
- vacuum_layer = ba.Layer(m_vacuum)
- vacuum_layer.addLayout(particle_layout)
- substrate_layer = ba.Layer(m_substrate)
-
- multi_layer = ba.MultiLayer()
- multi_layer.addLayer(vacuum_layer)
- multi_layer.addLayer(substrate_layer)
- return multi_layer
+
+ # Define materials
+ material_Particle = ba.HomogeneousMaterial("Particle", 0.0006, 2e-08)
+ material_Substrate = ba.HomogeneousMaterial("Substrate", 6e-06, 2e-08)
+ material_Vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
+
+ # Define form factors
+ ff_1 = ba.FormFactorFullSphere(2.0*nm)
+ ff_2 = ba.FormFactorCylinder(20.0*nm, 50.0*nm)
+
+ # Define particles
+ particle_1 = ba.Particle(material_Particle, ff_1)
+
+ # Define 3D lattices
+ lattice = ba.Lattice3D(ba.kvector_t(5.0*nm, 0.0*nm, 0.0*nm),
+ ba.kvector_t(0.0*nm, 5.0*nm, 0.0*nm),
+ ba.kvector_t(0.0*nm, 0.0*nm, 5.0*nm))
+
+ # Define crystals
+ crystal = ba.Crystal(particle_1, lattice)
+
+ # Define mesocrystals
+ particle_2 = ba.MesoCrystal(crystal, ff_2)
+
+ # Define particle layouts
+ layout = ba.ParticleLayout()
+ layout.addParticle(particle_2, 1.0)
+ layout.setWeight(1)
+ layout.setTotalParticleSurfaceDensity(0.01)
+
+ # Define layers
+ layer_1 = ba.Layer(material_Vacuum)
+ layer_1.addLayout(layout)
+ layer_2 = ba.Layer(material_Substrate)
+
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer_1)
+ sample.addLayer(layer_2)
+
+ return sample
def get_simulation():
diff --git a/Examples/Python/sim02_Complexes/ParticlesCrossingInterface.py b/Examples/Python/sim02_Complexes/ParticlesCrossingInterface.py
index 870133138cfce84342c7e9af11db3d01847e0565..4065dbd16b0bebe16f00d52bd210fcf489f372e0 100644
--- a/Examples/Python/sim02_Complexes/ParticlesCrossingInterface.py
+++ b/Examples/Python/sim02_Complexes/ParticlesCrossingInterface.py
@@ -17,7 +17,7 @@ For example, X or Y rotated particles can not yet cross interfaces (exception
will be thrown when trying to simulate such geometries).
"""
import bornagain as ba
-from bornagain import deg, angstrom, nm
+from bornagain import angstrom, deg, nm, nm2, kvector_t
phi_min, phi_max = -1.0, 1.0
alpha_min, alpha_max = 0.0, 2.0
@@ -27,37 +27,46 @@ def get_sample():
"""
Returns a sample with uncorrelated cylinders and prisms on a substrate.
"""
- # defining materials
- m_vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
- m_substrate = ba.HomogeneousMaterial("Substrate", 6e-6, 2e-8)
- m_particle = ba.HomogeneousMaterial("Particle", 6e-4, 2e-8)
-
- shift_down = 3*nm
-
- # collection of particles
- cylinder_ff = ba.FormFactorCylinder(5*nm, 5*nm)
- cylinder = ba.Particle(m_particle, cylinder_ff)
- cylinder.setPosition(0, 0, -shift_down)
-
- prism_ff = ba.FormFactorPrism3(10*nm, 5*nm)
- prism = ba.Particle(m_particle, prism_ff)
- prism.setPosition(0, 0, -shift_down)
-
- particle_layout = ba.ParticleLayout()
- particle_layout.addParticle(cylinder, 0.5)
- particle_layout.addParticle(prism, 0.5)
- interference = ba.InterferenceFunctionNone()
- particle_layout.setInterferenceFunction(interference)
-
- # vacuum layer with particles and substrate form multi layer
- vacuum_layer = ba.Layer(m_vacuum)
- vacuum_layer.addLayout(particle_layout)
- substrate_layer = ba.Layer(m_substrate)
- multi_layer = ba.MultiLayer()
- multi_layer.addLayer(vacuum_layer)
- multi_layer.addLayer(substrate_layer)
- print(multi_layer.treeToString())
- return multi_layer
+
+ # Define materials
+ material_Particle = ba.HomogeneousMaterial("Particle", 0.0006, 2e-08)
+ material_Substrate = ba.HomogeneousMaterial("Substrate", 6e-06, 2e-08)
+ material_Vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
+
+ # Define form factors
+ ff_1 = ba.FormFactorCylinder(5.0*nm, 5.0*nm)
+ ff_2 = ba.FormFactorPrism3(10.0*nm, 5.0*nm)
+
+ # Define particles
+ particle_1 = ba.Particle(material_Particle, ff_1)
+ particle_1_position = kvector_t(0.0*nm, 0.0*nm, -3.0*nm)
+ particle_1.setPosition(particle_1_position)
+ particle_2 = ba.Particle(material_Particle, ff_2)
+ particle_2_position = kvector_t(0.0*nm, 0.0*nm, -3.0*nm)
+ particle_2.setPosition(particle_2_position)
+
+ # Define interference functions
+ iff = ba.InterferenceFunctionNone()
+
+ # Define particle layouts
+ layout = ba.ParticleLayout()
+ layout.addParticle(particle_1, 0.5)
+ layout.addParticle(particle_2, 0.5)
+ layout.setInterferenceFunction(iff)
+ layout.setWeight(1)
+ layout.setTotalParticleSurfaceDensity(0.01)
+
+ # Define layers
+ layer_1 = ba.Layer(material_Vacuum)
+ layer_1.addLayout(layout)
+ layer_2 = ba.Layer(material_Substrate)
+
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer_1)
+ sample.addLayer(layer_2)
+
+ return sample
def get_simulation():
diff --git a/Examples/Python/sim03_Structures/ApproximationDA.py b/Examples/Python/sim03_Structures/ApproximationDA.py
index a6c5f9c4165b8f8845f4e3e704fa46e29e241051..c072968b804ba59f5effe416b447e620af14c038 100644
--- a/Examples/Python/sim03_Structures/ApproximationDA.py
+++ b/Examples/Python/sim03_Structures/ApproximationDA.py
@@ -2,7 +2,7 @@
Cylinders of two different sizes in Decoupling Approximation
"""
import bornagain as ba
-from bornagain import deg, angstrom, nm
+from bornagain import angstrom, deg, nm, nm2, kvector_t
def get_sample():
@@ -10,40 +10,44 @@ def get_sample():
Returns a sample with cylinders of two different sizes on a substrate.
The cylinder positions are modelled in Decoupling Approximation.
"""
- m_vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
- m_substrate = ba.HomogeneousMaterial("Substrate", 6e-6, 2e-8)
- m_particle = ba.HomogeneousMaterial("Particle", 6e-4, 2e-8)
-
- # cylindrical particle 1
- radius1 = 5*nm
- height1 = radius1
- cylinder_ff1 = ba.FormFactorCylinder(radius1, height1)
- cylinder1 = ba.Particle(m_particle, cylinder_ff1)
-
- # cylindrical particle 2
- radius2 = 8*nm
- height2 = radius2
- cylinder_ff2 = ba.FormFactorCylinder(radius2, height2)
- cylinder2 = ba.Particle(m_particle, cylinder_ff2)
-
- # interference function
- interference = ba.InterferenceFunctionRadialParaCrystal(18.0*nm, 1e3*nm)
- pdf = ba.FTDistribution1DGauss(3*nm)
- interference.setProbabilityDistribution(pdf)
-
- # assembling the sample
- particle_layout = ba.ParticleLayout()
- particle_layout.addParticle(cylinder1, 0.8)
- particle_layout.addParticle(cylinder2, 0.2)
- particle_layout.setInterferenceFunction(interference)
-
- vacuum_layer = ba.Layer(m_vacuum)
- vacuum_layer.addLayout(particle_layout)
- substrate_layer = ba.Layer(m_substrate)
- multi_layer = ba.MultiLayer()
- multi_layer.addLayer(vacuum_layer)
- multi_layer.addLayer(substrate_layer)
- return multi_layer
+
+ # Define materials
+ material_Particle = ba.HomogeneousMaterial("Particle", 0.0006, 2e-08)
+ material_Substrate = ba.HomogeneousMaterial("Substrate", 6e-06, 2e-08)
+ material_Vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
+
+ # Define form factors
+ ff_1 = ba.FormFactorCylinder(5.0*nm, 5.0*nm)
+ ff_2 = ba.FormFactorCylinder(8.0*nm, 8.0*nm)
+
+ # Define particles
+ particle_1 = ba.Particle(material_Particle, ff_1)
+ particle_2 = ba.Particle(material_Particle, ff_2)
+
+ # Define interference functions
+ iff = ba.InterferenceFunctionRadialParaCrystal(18.0*nm, 1000.0*nm)
+ iff_pdf = ba.FTDistribution1DGauss(3.0*nm)
+ iff.setProbabilityDistribution(iff_pdf)
+
+ # Define particle layouts
+ layout = ba.ParticleLayout()
+ layout.addParticle(particle_1, 0.8)
+ layout.addParticle(particle_2, 0.2)
+ layout.setInterferenceFunction(iff)
+ layout.setWeight(1)
+ layout.setTotalParticleSurfaceDensity(0.01)
+
+ # Define layers
+ layer_1 = ba.Layer(material_Vacuum)
+ layer_1.addLayout(layout)
+ layer_2 = ba.Layer(material_Substrate)
+
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer_1)
+ sample.addLayer(layer_2)
+
+ return sample
def get_simulation():
diff --git a/Examples/Python/sim03_Structures/ApproximationLMA.py b/Examples/Python/sim03_Structures/ApproximationLMA.py
index 7417b465df286e2c58b310ad6e9070d550c933d9..73c8eebb9b3947f23b9da1de47a4608d8b17e872 100644
--- a/Examples/Python/sim03_Structures/ApproximationLMA.py
+++ b/Examples/Python/sim03_Structures/ApproximationLMA.py
@@ -2,7 +2,7 @@
Cylinders of two different sizes in Local Monodisperse Approximation
"""
import bornagain as ba
-from bornagain import deg, angstrom, nm
+from bornagain import angstrom, deg, nm, nm2, kvector_t
def get_sample():
@@ -10,48 +10,52 @@ def get_sample():
Returns a sample with cylinders of two different sizes on a substrate.
The cylinder positions are modelled in Local Monodisperse Approximation.
"""
- m_vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
- m_substrate = ba.HomogeneousMaterial("Substrate", 6e-6, 2e-8)
- m_particle = ba.HomogeneousMaterial("Particle", 6e-4, 2e-8)
-
- # cylindrical particle 1
- radius1 = 5*nm
- height1 = radius1
- cylinder_ff1 = ba.FormFactorCylinder(radius1, height1)
- cylinder1 = ba.Particle(m_particle, cylinder_ff1)
-
- # cylindrical particle 2
- radius2 = 8*nm
- height2 = radius2
- cylinder_ff2 = ba.FormFactorCylinder(radius2, height2)
- cylinder2 = ba.Particle(m_particle, cylinder_ff2)
-
- # interference function1
- interference1 = ba.InterferenceFunctionRadialParaCrystal(16.8*nm, 1e3*nm)
- pdf = ba.FTDistribution1DGauss(3*nm)
- interference1.setProbabilityDistribution(pdf)
-
- # interference function2
- interference2 = ba.InterferenceFunctionRadialParaCrystal(22.8*nm, 1e3*nm)
- interference2.setProbabilityDistribution(pdf)
-
- # assembling the sample
- particle_layout1 = ba.ParticleLayout()
- particle_layout1.addParticle(cylinder1, 0.8)
- particle_layout1.setInterferenceFunction(interference1)
-
- particle_layout2 = ba.ParticleLayout()
- particle_layout2.addParticle(cylinder2, 0.2)
- particle_layout2.setInterferenceFunction(interference2)
-
- vacuum_layer = ba.Layer(m_vacuum)
- vacuum_layer.addLayout(particle_layout1)
- vacuum_layer.addLayout(particle_layout2)
- substrate_layer = ba.Layer(m_substrate)
- multi_layer = ba.MultiLayer()
- multi_layer.addLayer(vacuum_layer)
- multi_layer.addLayer(substrate_layer)
- return multi_layer
+
+ # Define materials
+ material_Particle = ba.HomogeneousMaterial("Particle", 0.0006, 2e-08)
+ material_Substrate = ba.HomogeneousMaterial("Substrate", 6e-06, 2e-08)
+ material_Vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
+
+ # Define form factors
+ ff_1 = ba.FormFactorCylinder(5.0*nm, 5.0*nm)
+ ff_2 = ba.FormFactorCylinder(8.0*nm, 8.0*nm)
+
+ # Define particles
+ particle_1 = ba.Particle(material_Particle, ff_1)
+ particle_2 = ba.Particle(material_Particle, ff_2)
+
+ # Define interference functions
+ iff_1 = ba.InterferenceFunctionRadialParaCrystal(16.8*nm, 1000.0*nm)
+ iff_1_pdf = ba.FTDistribution1DGauss(3.0*nm)
+ iff_1.setProbabilityDistribution(iff_1_pdf)
+ iff_2 = ba.InterferenceFunctionRadialParaCrystal(22.8*nm, 1000.0*nm)
+ iff_2_pdf = ba.FTDistribution1DGauss(3.0*nm)
+ iff_2.setProbabilityDistribution(iff_2_pdf)
+
+ # Define particle layouts
+ layout_1 = ba.ParticleLayout()
+ layout_1.addParticle(particle_1, 0.8)
+ layout_1.setInterferenceFunction(iff_1)
+ layout_1.setWeight(1)
+ layout_1.setTotalParticleSurfaceDensity(0.01)
+ layout_2 = ba.ParticleLayout()
+ layout_2.addParticle(particle_2, 0.2)
+ layout_2.setInterferenceFunction(iff_2)
+ layout_2.setWeight(1)
+ layout_2.setTotalParticleSurfaceDensity(0.01)
+
+ # Define layers
+ layer_1 = ba.Layer(material_Vacuum)
+ layer_1.addLayout(layout_1)
+ layer_1.addLayout(layout_2)
+ layer_2 = ba.Layer(material_Substrate)
+
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer_1)
+ sample.addLayer(layer_2)
+
+ return sample
def get_simulation():
diff --git a/Examples/Python/sim03_Structures/ApproximationSSCA.py b/Examples/Python/sim03_Structures/ApproximationSSCA.py
index 8948054bd8722e822fb8e82b6a0f0c6278b3ab2d..abc4a8cb83ba6dd9afedf41a54b0aa6ac54e53cf 100644
--- a/Examples/Python/sim03_Structures/ApproximationSSCA.py
+++ b/Examples/Python/sim03_Structures/ApproximationSSCA.py
@@ -2,7 +2,7 @@
Cylinders of two different sizes in Size-Spacing Coupling Approximation
"""
import bornagain as ba
-from bornagain import deg, angstrom, nm
+from bornagain import angstrom, deg, nm, nm2, kvector_t
def get_sample():
@@ -10,41 +10,45 @@ def get_sample():
Returns a sample with cylinders of two different sizes on a substrate.
The cylinder positions are modelled in Size-Spacing Coupling Approximation.
"""
- m_vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
- m_substrate = ba.HomogeneousMaterial("Substrate", 6e-6, 2e-8)
- m_particle = ba.HomogeneousMaterial("Particle", 6e-4, 2e-8)
-
- # cylindrical particle 1
- radius1 = 5*nm
- height1 = radius1
- cylinder_ff1 = ba.FormFactorCylinder(radius1, height1)
- cylinder1 = ba.Particle(m_particle, cylinder_ff1)
-
- # cylindrical particle 2
- radius2 = 8*nm
- height2 = radius2
- cylinder_ff2 = ba.FormFactorCylinder(radius2, height2)
- cylinder2 = ba.Particle(m_particle, cylinder_ff2)
-
- # interference function
- interference = ba.InterferenceFunctionRadialParaCrystal(18.0*nm, 1e3*nm)
- pdf = ba.FTDistribution1DGauss(3*nm)
- interference.setProbabilityDistribution(pdf)
- interference.setKappa(1.0)
-
- # assembling the sample
- particle_layout = ba.ParticleLayout()
- particle_layout.addParticle(cylinder1, 0.8)
- particle_layout.addParticle(cylinder2, 0.2)
- particle_layout.setInterferenceFunction(interference)
-
- vacuum_layer = ba.Layer(m_vacuum)
- vacuum_layer.addLayout(particle_layout)
- substrate_layer = ba.Layer(m_substrate)
- multi_layer = ba.MultiLayer()
- multi_layer.addLayer(vacuum_layer)
- multi_layer.addLayer(substrate_layer)
- return multi_layer
+
+ # Define materials
+ material_Particle = ba.HomogeneousMaterial("Particle", 0.0006, 2e-08)
+ material_Substrate = ba.HomogeneousMaterial("Substrate", 6e-06, 2e-08)
+ material_Vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
+
+ # Define form factors
+ ff_1 = ba.FormFactorCylinder(5.0*nm, 5.0*nm)
+ ff_2 = ba.FormFactorCylinder(8.0*nm, 8.0*nm)
+
+ # Define particles
+ particle_1 = ba.Particle(material_Particle, ff_1)
+ particle_2 = ba.Particle(material_Particle, ff_2)
+
+ # Define interference functions
+ iff = ba.InterferenceFunctionRadialParaCrystal(18.0*nm, 1000.0*nm)
+ iff.setKappa(1.0)
+ iff_pdf = ba.FTDistribution1DGauss(3.0*nm)
+ iff.setProbabilityDistribution(iff_pdf)
+
+ # Define particle layouts
+ layout = ba.ParticleLayout()
+ layout.addParticle(particle_1, 0.8)
+ layout.addParticle(particle_2, 0.2)
+ layout.setInterferenceFunction(iff)
+ layout.setWeight(1)
+ layout.setTotalParticleSurfaceDensity(0.01)
+
+ # Define layers
+ layer_1 = ba.Layer(material_Vacuum)
+ layer_1.addLayout(layout)
+ layer_2 = ba.Layer(material_Substrate)
+
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer_1)
+ sample.addLayer(layer_2)
+
+ return sample
def get_simulation():
diff --git a/Examples/Python/sim03_Structures/CosineRipplesAtRectLattice.py b/Examples/Python/sim03_Structures/CosineRipplesAtRectLattice.py
index 77a20d7b5340b2a4d8fcf94eb14006a52733f051..b4ae226890224b83dc62439efe2f945e434d0c20 100644
--- a/Examples/Python/sim03_Structures/CosineRipplesAtRectLattice.py
+++ b/Examples/Python/sim03_Structures/CosineRipplesAtRectLattice.py
@@ -3,7 +3,7 @@ Cosine ripple on a 2D lattice
"""
import numpy
import bornagain as ba
-from bornagain import deg, angstrom, nm
+from bornagain import angstrom, deg, nm, nm2, kvector_t
def get_sample():
@@ -11,33 +11,44 @@ def get_sample():
Returns a sample with cosine ripples on a substrate.
The structure is modelled as a 2D Lattice.
"""
- # defining materials
- m_vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
- m_substrate = ba.HomogeneousMaterial("Substrate", 6e-6, 2e-8)
- m_particle = ba.HomogeneousMaterial("Particle", 6e-4, 2e-8)
-
- # collection of particles
- ff = ba.FormFactorCosineRippleBox(100*nm, 20*nm, 4*nm)
- particle = ba.Particle(m_particle, ff)
-
- particle_layout = ba.ParticleLayout()
- particle_layout.addParticle(particle, 1.0)
-
- interference = ba.InterferenceFunction2DLattice(
- ba.BasicLattice2D(200.0*nm, 50.0*nm, 90.0*deg, 0.0*deg))
- pdf = ba.FTDecayFunction2DCauchy(160*nm, 16*nm, 0)
- interference.setDecayFunction(pdf)
- particle_layout.setInterferenceFunction(interference)
-
- # assemble the sample
- vacuum_layer = ba.Layer(m_vacuum)
- vacuum_layer.addLayout(particle_layout)
- substrate_layer = ba.Layer(m_substrate)
- multi_layer = ba.MultiLayer()
- multi_layer.addLayer(vacuum_layer)
- multi_layer.addLayer(substrate_layer)
-
- return multi_layer
+
+ # Define materials
+ material_Particle = ba.HomogeneousMaterial("Particle", 0.0006, 2e-08)
+ material_Substrate = ba.HomogeneousMaterial("Substrate", 6e-06, 2e-08)
+ material_Vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
+
+ # Define form factors
+ ff = ba.FormFactorCosineRippleBox(100.0*nm, 20.0*nm, 4.0*nm)
+
+ # Define particles
+ particle = ba.Particle(material_Particle, ff)
+
+ # Define 2D lattices
+ lattice = ba.BasicLattice2D(200.0*nm, 50.0*nm, 90.0*deg, 0.0*deg)
+
+ # Define interference functions
+ iff = ba.InterferenceFunction2DLattice(lattice)
+ iff_pdf = ba.FTDecayFunction2DCauchy(160.0*nm, 16.0*nm, 0.0*deg)
+ iff.setDecayFunction(iff_pdf)
+
+ # Define particle layouts
+ layout = ba.ParticleLayout()
+ layout.addParticle(particle, 1.0)
+ layout.setInterferenceFunction(iff)
+ layout.setWeight(1)
+ layout.setTotalParticleSurfaceDensity(0.0001)
+
+ # Define layers
+ layer_1 = ba.Layer(material_Vacuum)
+ layer_1.addLayout(layout)
+ layer_2 = ba.Layer(material_Substrate)
+
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer_1)
+ sample.addLayer(layer_2)
+
+ return sample
def get_simulation():
diff --git a/Examples/Python/sim03_Structures/Interference1DLattice.py b/Examples/Python/sim03_Structures/Interference1DLattice.py
index 5e9301b2688fa4c52d7db33ac39ea489a82180c0..a23eb272e885653ac054d3856755eda4b986c84e 100644
--- a/Examples/Python/sim03_Structures/Interference1DLattice.py
+++ b/Examples/Python/sim03_Structures/Interference1DLattice.py
@@ -4,7 +4,7 @@ Monte-carlo integration is used to get rid of
large-particle form factor oscillations.
"""
import bornagain as ba
-from bornagain import deg, angstrom, nm
+from bornagain import angstrom, deg, nm, nm2, kvector_t
def get_sample():
@@ -12,36 +12,43 @@ def get_sample():
Returns a sample with a grating on a substrate,
modelled by very long boxes forming a 1D lattice with Cauchy correlations.
"""
- # defining materials
- m_vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
- m_substrate = ba.HomogeneousMaterial("Substrate", 6e-6, 2e-8)
- m_particle = ba.HomogeneousMaterial("Particle", 6e-4, 2e-8)
-
- box_length, box_width, box_height = 10*nm, 10000*nm, 10*nm
- lattice_length = 30*nm
-
- # collection of particles
- interference = ba.InterferenceFunction1DLattice(lattice_length, 45*deg)
- pdf = ba.FTDecayFunction1DCauchy(1000.0)
- interference.setDecayFunction(pdf)
-
- box_ff = ba.FormFactorBox(box_length, box_width, box_height)
- box = ba.Particle(m_particle, box_ff)
-
- particle_layout = ba.ParticleLayout()
- particle_layout.addParticle(box, 1.0, ba.kvector_t(0.0, 0.0, 0.0),
- ba.RotationZ(45*deg))
- particle_layout.setInterferenceFunction(interference)
-
- # assembling the sample
- vacuum_layer = ba.Layer(m_vacuum)
- vacuum_layer.addLayout(particle_layout)
- substrate_layer = ba.Layer(m_substrate)
-
- multi_layer = ba.MultiLayer()
- multi_layer.addLayer(vacuum_layer)
- multi_layer.addLayer(substrate_layer)
- return multi_layer
+
+ # Define materials
+ material_Particle = ba.HomogeneousMaterial("Particle", 0.0006, 2e-08)
+ material_Substrate = ba.HomogeneousMaterial("Substrate", 6e-06, 2e-08)
+ material_Vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
+
+ # Define form factors
+ ff = ba.FormFactorBox(10.0*nm, 10000.0*nm, 10.0*nm)
+
+ # Define particles
+ particle = ba.Particle(material_Particle, ff)
+ particle_rotation = ba.RotationZ(45.0*deg)
+ particle.setRotation(particle_rotation)
+
+ # Define interference functions
+ iff = ba.InterferenceFunction1DLattice(30.0*nm, 45.0*deg)
+ iff_pdf = ba.FTDecayFunction1DCauchy(1000.0*nm)
+ iff.setDecayFunction(iff_pdf)
+
+ # Define particle layouts
+ layout = ba.ParticleLayout()
+ layout.addParticle(particle, 1.0)
+ layout.setInterferenceFunction(iff)
+ layout.setWeight(1)
+ layout.setTotalParticleSurfaceDensity(0.01)
+
+ # Define layers
+ layer_1 = ba.Layer(material_Vacuum)
+ layer_1.addLayout(layout)
+ layer_2 = ba.Layer(material_Substrate)
+
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer_1)
+ sample.addLayer(layer_2)
+
+ return sample
def get_simulation():
diff --git a/Examples/Python/sim03_Structures/Interference1DRadialParaCrystal.py b/Examples/Python/sim03_Structures/Interference1DRadialParaCrystal.py
index cadaf87bfafa237de006e92cf733744a2d712bbf..12e550a5a155149517eeca0327417769bcc95125 100644
--- a/Examples/Python/sim03_Structures/Interference1DRadialParaCrystal.py
+++ b/Examples/Python/sim03_Structures/Interference1DRadialParaCrystal.py
@@ -2,7 +2,7 @@
radial paracrystal
"""
import bornagain as ba
-from bornagain import deg, angstrom, nm
+from bornagain import angstrom, deg, nm, nm2, kvector_t
phi_min, phi_max = -2.0, 2.0
alpha_min, alpha_max = 0.0, 2.0
@@ -12,33 +12,41 @@ def get_sample():
"""
Returns a sample with cylinders on a substrate that form a radial paracrystal.
"""
- # defining materials
- m_vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
- m_substrate = ba.HomogeneousMaterial("Substrate", 6e-6, 2e-8)
- m_particle = ba.HomogeneousMaterial("Particle", 6e-4, 2e-8)
-
- # collection of particles
- cylinder_ff = ba.FormFactorCylinder(5*nm, 5*nm)
- cylinder = ba.Particle(m_particle, cylinder_ff)
-
- interference = ba.InterferenceFunctionRadialParaCrystal(20.0*nm, 1e3*nm)
- pdf = ba.FTDistribution1DGauss(7*nm)
- interference.setProbabilityDistribution(pdf)
-
- particle_layout = ba.ParticleLayout()
- particle_layout.addParticle(cylinder, 1.0)
- particle_layout.setInterferenceFunction(interference)
-
- # assembling the sample
- vacuum_layer = ba.Layer(m_vacuum)
- vacuum_layer.addLayout(particle_layout)
- substrate_layer = ba.Layer(m_substrate)
-
- multi_layer = ba.MultiLayer()
- multi_layer.addLayer(vacuum_layer)
- multi_layer.addLayer(substrate_layer)
- # print(multi_layer.treeToString())
- return multi_layer
+
+ # Define materials
+ material_Particle = ba.HomogeneousMaterial("Particle", 0.0006, 2e-08)
+ material_Substrate = ba.HomogeneousMaterial("Substrate", 6e-06, 2e-08)
+ material_Vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
+
+ # Define form factors
+ ff = ba.FormFactorCylinder(5.0*nm, 5.0*nm)
+
+ # Define particles
+ particle = ba.Particle(material_Particle, ff)
+
+ # Define interference functions
+ iff = ba.InterferenceFunctionRadialParaCrystal(20.0*nm, 1000.0*nm)
+ iff_pdf = ba.FTDistribution1DGauss(7.0*nm)
+ iff.setProbabilityDistribution(iff_pdf)
+
+ # Define particle layouts
+ layout = ba.ParticleLayout()
+ layout.addParticle(particle, 1.0)
+ layout.setInterferenceFunction(iff)
+ layout.setWeight(1)
+ layout.setTotalParticleSurfaceDensity(0.01)
+
+ # Define layers
+ layer_1 = ba.Layer(material_Vacuum)
+ layer_1.addLayout(layout)
+ layer_2 = ba.Layer(material_Substrate)
+
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer_1)
+ sample.addLayer(layer_2)
+
+ return sample
def get_simulation():
diff --git a/Examples/Python/sim03_Structures/Interference2DCenteredSquareLattice.py b/Examples/Python/sim03_Structures/Interference2DCenteredSquareLattice.py
index cba7b06e579b352b77a7be10af0f6b5a8ecd0535..03049b3fb7b280223d910daf2450e1c90088194a 100644
--- a/Examples/Python/sim03_Structures/Interference2DCenteredSquareLattice.py
+++ b/Examples/Python/sim03_Structures/Interference2DCenteredSquareLattice.py
@@ -3,7 +3,7 @@
"""
import numpy
import bornagain as ba
-from bornagain import deg, angstrom, nm
+from bornagain import angstrom, deg, nm, nm2, kvector_t
def get_sample():
@@ -11,38 +11,53 @@ def get_sample():
Returns a sample with cylinders on a substrate,
forming a 2D centered square lattice
"""
- # defining materials
- m_vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
- m_substrate = ba.HomogeneousMaterial("Substrate", 6e-6, 2e-8)
- m_particle = ba.HomogeneousMaterial("Particle", 6e-4, 2e-8)
-
- # collection of particles
- interference = ba.InterferenceFunction2DLattice(
- ba.SquareLattice2D(25.0*nm, 0))
- pdf = ba.FTDecayFunction2DCauchy(48*nm, 16*nm, 0)
- interference.setDecayFunction(pdf)
-
- particle_layout = ba.ParticleLayout()
- position1 = ba.kvector_t(0.0, 0.0, 0.0)
- position2 = ba.kvector_t(12.5*nm, 12.5*nm, 0.0)
- cylinder_ff = ba.FormFactorCylinder(3.*nm, 3.*nm)
- cylinder = ba.Particle(m_particle, cylinder_ff)
- basis = ba.ParticleComposition()
- basis.addParticles(cylinder, [position1, position2])
- particle_layout.addParticle(basis)
- particle_layout.setInterferenceFunction(interference)
-
- # assembling the sample
- vacuum_layer = ba.Layer(m_vacuum)
- vacuum_layer.addLayout(particle_layout)
- substrate_layer = ba.Layer(m_substrate)
-
- multi_layer = ba.MultiLayer()
- multi_layer.addLayer(vacuum_layer)
- multi_layer.addLayer(substrate_layer)
-
- print(multi_layer.treeToString())
- return multi_layer
+
+ # Define materials
+ material_Particle = ba.HomogeneousMaterial("Particle", 0.0006, 2e-08)
+ material_Substrate = ba.HomogeneousMaterial("Substrate", 6e-06, 2e-08)
+ material_Vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
+
+ # Define form factors
+ ff_1 = ba.FormFactorCylinder(3.0*nm, 3.0*nm)
+ ff_2 = ba.FormFactorCylinder(3.0*nm, 3.0*nm)
+
+ # Define particles
+ particle_1 = ba.Particle(material_Particle, ff_1)
+ particle_2 = ba.Particle(material_Particle, ff_2)
+ particle_2_position = kvector_t(12.5*nm, 12.5*nm, 0.0*nm)
+ particle_2.setPosition(particle_2_position)
+
+ # Define composition of particles at specific positions
+ particle_3 = ba.ParticleComposition()
+ particle_3.addParticle(particle_1)
+ particle_3.addParticle(particle_2)
+
+ # Define 2D lattices
+ lattice = ba.BasicLattice2D(25.0*nm, 25.0*nm, 90.0*deg, 0.0*deg)
+
+ # Define interference functions
+ iff = ba.InterferenceFunction2DLattice(lattice)
+ iff_pdf = ba.FTDecayFunction2DCauchy(48.0*nm, 16.0*nm, 0.0*deg)
+ iff.setDecayFunction(iff_pdf)
+
+ # Define particle layouts
+ layout = ba.ParticleLayout()
+ layout.addParticle(particle_3, 1.0)
+ layout.setInterferenceFunction(iff)
+ layout.setWeight(1)
+ layout.setTotalParticleSurfaceDensity(0.0016)
+
+ # Define layers
+ layer_1 = ba.Layer(material_Vacuum)
+ layer_1.addLayout(layout)
+ layer_2 = ba.Layer(material_Substrate)
+
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer_1)
+ sample.addLayer(layer_2)
+
+ return sample
def get_simulation():
diff --git a/Examples/Python/sim03_Structures/Interference2DParaCrystal.py b/Examples/Python/sim03_Structures/Interference2DParaCrystal.py
index 7a36395a45cc4f7ba2e4c46aebab24f9ee1c1ea8..dfae0d5c546f0ba8c4e51376d145db064bdf33d3 100644
--- a/Examples/Python/sim03_Structures/Interference2DParaCrystal.py
+++ b/Examples/Python/sim03_Structures/Interference2DParaCrystal.py
@@ -2,43 +2,54 @@
2D paracrystal
"""
import bornagain as ba
-from bornagain import deg, angstrom, nm, micrometer
+from bornagain import angstrom, deg, nm, nm2, kvector_t
def get_sample():
"""
Returns a sample with cylinders on a substrate, forming a 2D paracrystal
"""
- m_vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
- m_substrate = ba.HomogeneousMaterial("Substrate", 6e-6, 2e-8)
- m_particle = ba.HomogeneousMaterial("Particle", 6e-4, 2e-8)
-
- # collection of particles
- cylinder_ff = ba.FormFactorCylinder(4*nm, 5*nm)
- cylinder = ba.Particle(m_particle, cylinder_ff)
-
- interference = ba.InterferenceFunction2DParaCrystal(
- ba.SquareLattice2D(10.0*nm), 0.0, 20.0*micrometer, 20.0*micrometer)
- interference.setIntegrationOverXi(True)
- pdf = ba.FTDistribution2DCauchy(1.0*nm, 1.0*nm, 0)
- interference.setProbabilityDistributions(pdf, pdf)
-
- particle_layout = ba.ParticleLayout()
- particle_layout.addParticle(cylinder, 1.0)
- particle_layout.setInterferenceFunction(interference)
-
- # assembling the sample
- vacuum_layer = ba.Layer(m_vacuum)
- vacuum_layer.addLayout(particle_layout)
-
- substrate_layer = ba.Layer(m_substrate)
-
- multi_layer = ba.MultiLayer()
- multi_layer.addLayer(vacuum_layer)
- multi_layer.addLayer(substrate_layer)
- print(multi_layer.parametersToString())
- print(multi_layer.treeToString())
- return multi_layer
+
+ # Define materials
+ material_Particle = ba.HomogeneousMaterial("Particle", 0.0006, 2e-08)
+ material_Substrate = ba.HomogeneousMaterial("Substrate", 6e-06, 2e-08)
+ material_Vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
+
+ # Define form factors
+ ff = ba.FormFactorCylinder(4.0*nm, 5.0*nm)
+
+ # Define particles
+ particle = ba.Particle(material_Particle, ff)
+
+ # Define 2D lattices
+ lattice = ba.BasicLattice2D(10.0*nm, 10.0*nm, 90.0*deg, 0.0*deg)
+
+ # Define interference functions
+ iff = ba.InterferenceFunction2DParaCrystal(lattice, 0.0*nm, 20000.0*nm,
+ 20000.0*nm)
+ iff.setIntegrationOverXi(True)
+ iff_pdf_1 = ba.FTDistribution2DCauchy(1.0*nm, 1.0*nm, 0.0*deg)
+ iff_pdf_2 = ba.FTDistribution2DCauchy(1.0*nm, 1.0*nm, 0.0*deg)
+ iff.setProbabilityDistributions(iff_pdf_1, iff_pdf_2)
+
+ # Define particle layouts
+ layout = ba.ParticleLayout()
+ layout.addParticle(particle, 1.0)
+ layout.setInterferenceFunction(iff)
+ layout.setWeight(1)
+ layout.setTotalParticleSurfaceDensity(0.01)
+
+ # Define layers
+ layer_1 = ba.Layer(material_Vacuum)
+ layer_1.addLayout(layout)
+ layer_2 = ba.Layer(material_Substrate)
+
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer_1)
+ sample.addLayer(layer_2)
+
+ return sample
def get_simulation():
diff --git a/Examples/Python/sim03_Structures/Interference2DRotatedSquareLattice.py b/Examples/Python/sim03_Structures/Interference2DRotatedSquareLattice.py
index 77e4bac2e6cbeef19984ad02b5c4a6051c33f242..33ec3f79809db7ed9ef3666054966b6723863b70 100644
--- a/Examples/Python/sim03_Structures/Interference2DRotatedSquareLattice.py
+++ b/Examples/Python/sim03_Structures/Interference2DRotatedSquareLattice.py
@@ -3,41 +3,52 @@ Cylinders on a rotated 2D lattice
"""
import numpy
import bornagain as ba
-from bornagain import deg, angstrom, nm
+from bornagain import angstrom, deg, nm, nm2, kvector_t
def get_sample():
"""
Returns a sample with cylinders on a substrate, forming a rotated 2D lattice
"""
- # defining materials
- m_vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
- m_substrate = ba.HomogeneousMaterial("Substrate", 6e-6, 2e-8)
- m_particle = ba.HomogeneousMaterial("Particle", 6e-4, 2e-8)
-
- # collection of particles
- interference = ba.InterferenceFunction2DLattice(
- ba.SquareLattice2D(25.0*nm, 30.0*deg))
- pdf = ba.FTDecayFunction2DCauchy(300.0*nm/2.0/numpy.pi,
- 100.0*nm/2.0/numpy.pi, 0)
- pdf.setParameterValue("Gamma", 30.0*deg)
- interference.setDecayFunction(pdf)
-
- cylinder_ff = ba.FormFactorCylinder(3.*nm, 3.*nm)
- cylinder = ba.Particle(m_particle, cylinder_ff)
- particle_layout = ba.ParticleLayout()
- particle_layout.addParticle(cylinder)
- particle_layout.setInterferenceFunction(interference)
-
- # assembling the sample
- vacuum_layer = ba.Layer(m_vacuum)
- vacuum_layer.addLayout(particle_layout)
-
- substrate_layer = ba.Layer(m_substrate)
- multi_layer = ba.MultiLayer()
- multi_layer.addLayer(vacuum_layer)
- multi_layer.addLayer(substrate_layer)
- return multi_layer
+
+ # Define materials
+ material_Particle = ba.HomogeneousMaterial("Particle", 0.0006, 2e-08)
+ material_Substrate = ba.HomogeneousMaterial("Substrate", 6e-06, 2e-08)
+ material_Vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
+
+ # Define form factors
+ ff = ba.FormFactorCylinder(3.0*nm, 3.0*nm)
+
+ # Define particles
+ particle = ba.Particle(material_Particle, ff)
+
+ # Define 2D lattices
+ lattice = ba.BasicLattice2D(25.0*nm, 25.0*nm, 90.0*deg, 30.0*deg)
+
+ # Define interference functions
+ iff = ba.InterferenceFunction2DLattice(lattice)
+ iff_pdf = ba.FTDecayFunction2DCauchy(47.7464829276*nm, 15.9154943092*nm,
+ 30.0*deg)
+ iff.setDecayFunction(iff_pdf)
+
+ # Define particle layouts
+ layout = ba.ParticleLayout()
+ layout.addParticle(particle, 1.0)
+ layout.setInterferenceFunction(iff)
+ layout.setWeight(1)
+ layout.setTotalParticleSurfaceDensity(0.0016)
+
+ # Define layers
+ layer_1 = ba.Layer(material_Vacuum)
+ layer_1.addLayout(layout)
+ layer_2 = ba.Layer(material_Substrate)
+
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer_1)
+ sample.addLayer(layer_2)
+
+ return sample
def get_simulation():
diff --git a/Examples/Python/sim03_Structures/Interference2DSquareFiniteLattice.py b/Examples/Python/sim03_Structures/Interference2DSquareFiniteLattice.py
index 937cfc666c07e551253fb739d74c42fd86328950..e018b063c317bb9935640c8dce3f10b1fd4ffff5 100644
--- a/Examples/Python/sim03_Structures/Interference2DSquareFiniteLattice.py
+++ b/Examples/Python/sim03_Structures/Interference2DSquareFiniteLattice.py
@@ -3,40 +3,50 @@ Cylinders on a 2D square lattice
"""
import numpy
import bornagain as ba
-from bornagain import deg, angstrom, nm
+from bornagain import angstrom, deg, nm, nm2, kvector_t
def get_sample():
"""
Returns a sample with cylinders on a substrate, forming a 2D square lattice.
"""
- # defining materials
- m_vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
- m_substrate = ba.HomogeneousMaterial("Substrate", 6e-6, 2e-8)
- m_particle = ba.HomogeneousMaterial("Particle", 6e-4, 2e-8)
-
- # collection of particles
- interference = ba.InterferenceFunctionFinite2DLattice(
- ba.SquareLattice2D(25.0*nm, 0.0), 40, 40)
- interference.setPositionVariance(1.0)
-
- cylinder_ff = ba.FormFactorCylinder(3.*nm, 3.*nm)
- cylinder = ba.Particle(m_particle, cylinder_ff)
- particle_layout = ba.ParticleLayout()
- particle_layout.addParticle(cylinder)
- particle_layout.setInterferenceFunction(interference)
-
- # assembling the sample
- vacuum_layer = ba.Layer(m_vacuum)
- vacuum_layer.addLayout(particle_layout)
- substrate_layer = ba.Layer(m_substrate)
-
- multi_layer = ba.MultiLayer()
- multi_layer.addLayer(vacuum_layer)
- multi_layer.addLayer(substrate_layer)
- print(multi_layer.parametersToString())
- print(multi_layer.treeToString())
- return multi_layer
+
+ # Define materials
+ material_Particle = ba.HomogeneousMaterial("Particle", 0.0006, 2e-08)
+ material_Substrate = ba.HomogeneousMaterial("Substrate", 6e-06, 2e-08)
+ material_Vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
+
+ # Define form factors
+ ff = ba.FormFactorCylinder(3.0*nm, 3.0*nm)
+
+ # Define particles
+ particle = ba.Particle(material_Particle, ff)
+
+ # Define 2D lattices
+ lattice = ba.BasicLattice2D(25.0*nm, 25.0*nm, 90.0*deg, 0.0*deg)
+
+ # Define interference functions
+ iff = ba.InterferenceFunctionFinite2DLattice(lattice, 40, 40)
+ iff.setPositionVariance(1.0*nm2)
+
+ # Define particle layouts
+ layout = ba.ParticleLayout()
+ layout.addParticle(particle, 1.0)
+ layout.setInterferenceFunction(iff)
+ layout.setWeight(1)
+ layout.setTotalParticleSurfaceDensity(0.0016)
+
+ # Define layers
+ layer_1 = ba.Layer(material_Vacuum)
+ layer_1.addLayout(layout)
+ layer_2 = ba.Layer(material_Substrate)
+
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer_1)
+ sample.addLayer(layer_2)
+
+ return sample
def get_simulation():
diff --git a/Examples/Python/sim03_Structures/Interference2DSquareLattice.py b/Examples/Python/sim03_Structures/Interference2DSquareLattice.py
index 37d7668d657953e2367eb552d6ee82f9ecd24e99..affd7715ec59e66c594c4acf5a22332d0da8d990 100644
--- a/Examples/Python/sim03_Structures/Interference2DSquareLattice.py
+++ b/Examples/Python/sim03_Structures/Interference2DSquareLattice.py
@@ -3,41 +3,51 @@ Cylinders on a 2D square lattice
"""
import numpy
import bornagain as ba
-from bornagain import deg, angstrom, nm
+from bornagain import angstrom, deg, nm, nm2, kvector_t
def get_sample():
"""
Returns a sample with cylinders on a substrate, forming a 2D square lattice.
"""
- # defining materials
- m_vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
- m_substrate = ba.HomogeneousMaterial("Substrate", 6e-6, 2e-8)
- m_particle = ba.HomogeneousMaterial("Particle", 6e-4, 2e-8)
-
- # collection of particles
- interference = ba.InterferenceFunction2DLattice(
- ba.SquareLattice2D(25.0*nm, 0*deg))
- pdf = ba.FTDecayFunction2DCauchy(48*nm, 16*nm, 0)
- interference.setDecayFunction(pdf)
-
- cylinder_ff = ba.FormFactorCylinder(3.*nm, 3.*nm)
- cylinder = ba.Particle(m_particle, cylinder_ff)
- particle_layout = ba.ParticleLayout()
- particle_layout.addParticle(cylinder)
- particle_layout.setInterferenceFunction(interference)
-
- # assembling the sample
- vacuum_layer = ba.Layer(m_vacuum)
- vacuum_layer.addLayout(particle_layout)
- substrate_layer = ba.Layer(m_substrate)
-
- multi_layer = ba.MultiLayer()
- multi_layer.addLayer(vacuum_layer)
- multi_layer.addLayer(substrate_layer)
- print(multi_layer.parametersToString())
- print(multi_layer.treeToString())
- return multi_layer
+
+ # Define materials
+ material_Particle = ba.HomogeneousMaterial("Particle", 0.0006, 2e-08)
+ material_Substrate = ba.HomogeneousMaterial("Substrate", 6e-06, 2e-08)
+ material_Vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
+
+ # Define form factors
+ ff = ba.FormFactorCylinder(3.0*nm, 3.0*nm)
+
+ # Define particles
+ particle = ba.Particle(material_Particle, ff)
+
+ # Define 2D lattices
+ lattice = ba.BasicLattice2D(25.0*nm, 25.0*nm, 90.0*deg, 0.0*deg)
+
+ # Define interference functions
+ iff = ba.InterferenceFunction2DLattice(lattice)
+ iff_pdf = ba.FTDecayFunction2DCauchy(48.0*nm, 16.0*nm, 0.0*deg)
+ iff.setDecayFunction(iff_pdf)
+
+ # Define particle layouts
+ layout = ba.ParticleLayout()
+ layout.addParticle(particle, 1.0)
+ layout.setInterferenceFunction(iff)
+ layout.setWeight(1)
+ layout.setTotalParticleSurfaceDensity(0.0016)
+
+ # Define layers
+ layer_1 = ba.Layer(material_Vacuum)
+ layer_1.addLayout(layout)
+ layer_2 = ba.Layer(material_Substrate)
+
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer_1)
+ sample.addLayer(layer_2)
+
+ return sample
def get_simulation():
diff --git a/Examples/Python/sim03_Structures/SpheresAtHexLattice.py b/Examples/Python/sim03_Structures/SpheresAtHexLattice.py
index 61dfea83cdebb8a497f50c002c03fe99b933dcc1..7f3c49d042b4f7e05fa5b942c17900b3816cb40e 100644
--- a/Examples/Python/sim03_Structures/SpheresAtHexLattice.py
+++ b/Examples/Python/sim03_Structures/SpheresAtHexLattice.py
@@ -2,7 +2,7 @@
Spheres on a hexagonal lattice
"""
import bornagain as ba
-from bornagain import deg, angstrom, nm
+from bornagain import angstrom, deg, nm, nm2, kvector_t
def get_sample():
@@ -10,29 +10,44 @@ def get_sample():
Returns a sample with spherical particles on a substrate,
forming a hexagonal 2D lattice.
"""
- m_vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
- m_substrate = ba.HomogeneousMaterial("Substrate", 6e-6, 2e-8)
- m_particle = ba.HomogeneousMaterial("Particle", 6e-4, 2e-8)
-
- sphere_ff = ba.FormFactorFullSphere(10.0*nm)
- sphere = ba.Particle(m_particle, sphere_ff)
- particle_layout = ba.ParticleLayout()
- particle_layout.addParticle(sphere)
-
- interference = ba.InterferenceFunction2DLattice(
- ba.HexagonalLattice2D(20.0*nm, 0*deg))
- pdf = ba.FTDecayFunction2DCauchy(10*nm, 10*nm, 0)
- interference.setDecayFunction(pdf)
-
- particle_layout.setInterferenceFunction(interference)
-
- vacuum_layer = ba.Layer(m_vacuum)
- vacuum_layer.addLayout(particle_layout)
- substrate_layer = ba.Layer(m_substrate, 0)
- multi_layer = ba.MultiLayer()
- multi_layer.addLayer(vacuum_layer)
- multi_layer.addLayer(substrate_layer)
- return multi_layer
+
+ # Define materials
+ material_Particle = ba.HomogeneousMaterial("Particle", 0.0006, 2e-08)
+ material_Substrate = ba.HomogeneousMaterial("Substrate", 6e-06, 2e-08)
+ material_Vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
+
+ # Define form factors
+ ff = ba.FormFactorFullSphere(10.0*nm)
+
+ # Define particles
+ particle = ba.Particle(material_Particle, ff)
+
+ # Define 2D lattices
+ lattice = ba.BasicLattice2D(20.0*nm, 20.0*nm, 120.0*deg, 0.0*deg)
+
+ # Define interference functions
+ iff = ba.InterferenceFunction2DLattice(lattice)
+ iff_pdf = ba.FTDecayFunction2DCauchy(10.0*nm, 10.0*nm, 0.0*deg)
+ iff.setDecayFunction(iff_pdf)
+
+ # Define particle layouts
+ layout = ba.ParticleLayout()
+ layout.addParticle(particle, 1.0)
+ layout.setInterferenceFunction(iff)
+ layout.setWeight(1)
+ layout.setTotalParticleSurfaceDensity(0.00288675134595)
+
+ # Define layers
+ layer_1 = ba.Layer(material_Vacuum)
+ layer_1.addLayout(layout)
+ layer_2 = ba.Layer(material_Substrate)
+
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer_1)
+ sample.addLayer(layer_2)
+
+ return sample
def get_simulation():
diff --git a/Examples/Python/sim03_Structures/TriangularRipple.py b/Examples/Python/sim03_Structures/TriangularRipple.py
index 88f23668861e2d03f21c2b989338b40c6740c8aa..2815287042a10a22f4e7cb04849b1132650fd8b9 100644
--- a/Examples/Python/sim03_Structures/TriangularRipple.py
+++ b/Examples/Python/sim03_Structures/TriangularRipple.py
@@ -3,7 +3,7 @@
"""
import numpy
import bornagain as ba
-from bornagain import deg, angstrom, nm
+from bornagain import angstrom, deg, nm, nm2, kvector_t
def get_sample():
@@ -11,34 +11,45 @@ def get_sample():
Returns a sample with a grating on a substrate, modelled by triangular ripples
forming a 1D Paracrystal.
"""
- # defining materials
- m_vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
- m_substrate = ba.HomogeneousMaterial("Substrate", 6e-6, 2e-8)
- m_particle = ba.HomogeneousMaterial("Particle", 6e-4, 2e-8)
-
- # collection of particles
- ripple2_ff = ba.FormFactorSawtoothRippleBox(100*nm, 20*nm, 4*nm, -3.0*nm)
- ripple = ba.Particle(m_particle, ripple2_ff)
-
- particle_layout = ba.ParticleLayout()
- particle_layout.addParticle(ripple, 1.0)
-
- interference = ba.InterferenceFunction2DLattice(
- ba.BasicLattice2D(200.0*nm, 50.0*nm, 90.0*deg, 0.0*deg))
- pdf = ba.FTDecayFunction2DGauss(1000.*nm/2./numpy.pi, 100.*nm/2./numpy.pi,
- 0)
- interference.setDecayFunction(pdf)
- particle_layout.setInterferenceFunction(interference)
-
- # vacuum layer with particles and substrate form multi layer
- vacuum_layer = ba.Layer(m_vacuum)
- vacuum_layer.addLayout(particle_layout)
- substrate_layer = ba.Layer(m_substrate, 0)
- multi_layer = ba.MultiLayer()
- multi_layer.addLayer(vacuum_layer)
- multi_layer.addLayer(substrate_layer)
-
- return multi_layer
+
+ # Define materials
+ material_Particle = ba.HomogeneousMaterial("Particle", 0.0006, 2e-08)
+ material_Substrate = ba.HomogeneousMaterial("Substrate", 6e-06, 2e-08)
+ material_Vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
+
+ # Define form factors
+ ff = ba.FormFactorSawtoothRippleBox(100.0*nm, 20.0*nm, 4.0*nm, -3.0*nm)
+
+ # Define particles
+ particle = ba.Particle(material_Particle, ff)
+
+ # Define 2D lattices
+ lattice = ba.BasicLattice2D(200.0*nm, 50.0*nm, 90.0*deg, 0.0*deg)
+
+ # Define interference functions
+ iff = ba.InterferenceFunction2DLattice(lattice)
+ iff_pdf = ba.FTDecayFunction2DGauss(159.154943092*nm, 15.9154943092*nm,
+ 0.0*deg)
+ iff.setDecayFunction(iff_pdf)
+
+ # Define particle layouts
+ layout = ba.ParticleLayout()
+ layout.addParticle(particle, 1.0)
+ layout.setInterferenceFunction(iff)
+ layout.setWeight(1)
+ layout.setTotalParticleSurfaceDensity(0.0001)
+
+ # Define layers
+ layer_1 = ba.Layer(material_Vacuum)
+ layer_1.addLayout(layout)
+ layer_2 = ba.Layer(material_Substrate)
+
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer_1)
+ sample.addLayer(layer_2)
+
+ return sample
def get_simulation():
diff --git a/Examples/Python/sim04_Multilayers/BuriedParticles.py b/Examples/Python/sim04_Multilayers/BuriedParticles.py
index 64a9e3c72b5d035020ee3f7985a579930cc16073..d543571da20b84d28b480c3c3dcc2cec353c3d3c 100644
--- a/Examples/Python/sim04_Multilayers/BuriedParticles.py
+++ b/Examples/Python/sim04_Multilayers/BuriedParticles.py
@@ -2,7 +2,7 @@
Spherical particles embedded in the middle of the layer on top of substrate.
"""
import bornagain as ba
-from bornagain import deg, angstrom, nm
+from bornagain import angstrom, deg, nm, nm2, kvector_t
def get_sample():
@@ -10,30 +10,41 @@ def get_sample():
Returns a sample with spherical particles in a layer
between vacuum and substrate.
"""
- # defining materials
- m_vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
- m_interm_layer = ba.HomogeneousMaterial("IntermLayer", 3.45e-6, 5.24e-9)
- m_substrate = ba.HomogeneousMaterial("Substrate", 7.43e-6, 1.72e-7)
- m_particle = ba.HomogeneousMaterial("Particle", 0.0, 0.0)
- # collection of particles
+ # Define materials
+ material_IntermLayer = ba.HomogeneousMaterial("IntermLayer", 3.45e-06,
+ 5.24e-09)
+ material_Particle = ba.HomogeneousMaterial("Particle", 0.0, 0.0)
+ material_Substrate = ba.HomogeneousMaterial("Substrate", 7.43e-06, 1.72e-07)
+ material_Vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
+
+ # Define form factors
ff = ba.FormFactorFullSphere(10.2*nm)
- sphere = ba.Particle(m_particle, ff)
- sphere.setPosition(0.0, 0.0, -25.2)
- particle_layout = ba.ParticleLayout()
- particle_layout.addParticle(sphere, 1.0)
-
- # assembling the sample
- vacuum_layer = ba.Layer(m_vacuum)
- intermediate_layer = ba.Layer(m_interm_layer, 30.*nm)
- intermediate_layer.addLayout(particle_layout)
- substrate_layer = ba.Layer(m_substrate, 0)
-
- multi_layer = ba.MultiLayer()
- multi_layer.addLayer(vacuum_layer)
- multi_layer.addLayer(intermediate_layer)
- multi_layer.addLayer(substrate_layer)
- return multi_layer
+
+ # Define particles
+ particle = ba.Particle(material_Particle, ff)
+ particle_position = kvector_t(0.0*nm, 0.0*nm, -25.2*nm)
+ particle.setPosition(particle_position)
+
+ # Define particle layouts
+ layout = ba.ParticleLayout()
+ layout.addParticle(particle, 1.0)
+ layout.setWeight(1)
+ layout.setTotalParticleSurfaceDensity(0.01)
+
+ # Define layers
+ layer_1 = ba.Layer(material_Vacuum)
+ layer_2 = ba.Layer(material_IntermLayer, 30.0*nm)
+ layer_2.addLayout(layout)
+ layer_3 = ba.Layer(material_Substrate)
+
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer_1)
+ sample.addLayer(layer_2)
+ sample.addLayer(layer_3)
+
+ return sample
def get_simulation():
diff --git a/Examples/Python/sim04_Multilayers/HalfSpheresInAverageTopLayer.py b/Examples/Python/sim04_Multilayers/HalfSpheresInAverageTopLayer.py
index ff0079a2ef201aced9c1c09b05021c1fd984a61b..a691fbab52d14a78ea9c97b1a90129aade4781e1 100644
--- a/Examples/Python/sim04_Multilayers/HalfSpheresInAverageTopLayer.py
+++ b/Examples/Python/sim04_Multilayers/HalfSpheresInAverageTopLayer.py
@@ -3,7 +3,7 @@ Square lattice of half spheres on substrate with usage of average material
and slicing
"""
import bornagain as ba
-from bornagain import deg, angstrom, nm
+from bornagain import angstrom, deg, nm, nm2, kvector_t
sphere_radius = 5*nm
n_slices = 10
@@ -13,35 +13,45 @@ def get_sample():
"""
Returns a sample with cylinders on a substrate.
"""
- # defining materials
- m_vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
- m_layer = ba.HomogeneousMaterial("Layer", 3e-6, 2e-8)
- m_substrate = ba.HomogeneousMaterial("Substrate", 6e-6, 2e-8)
- m_particle = ba.HomogeneousMaterial("Particle", 3e-5, 2e-8)
-
- # cylindrical particle
- half_sphere_ff = ba.FormFactorTruncatedSphere(sphere_radius, sphere_radius,
- 0)
- half_sphere = ba.Particle(m_particle, half_sphere_ff)
- particle_layout = ba.ParticleLayout()
- particle_layout.addParticle(half_sphere)
-
- # interference function
- interference = ba.InterferenceFunction2DLattice(
- ba.SquareLattice2D(10*nm, 0*deg))
- pdf = ba.FTDecayFunction2DCauchy(100*nm, 100*nm, 0)
- interference.setDecayFunction(pdf)
- particle_layout.setInterferenceFunction(interference)
-
- vacuum_layer = ba.Layer(m_vacuum)
- vacuum_layer.addLayout(particle_layout)
- vacuum_layer.setNumberOfSlices(n_slices)
- substrate_layer = ba.Layer(m_substrate)
-
- multi_layer = ba.MultiLayer()
- multi_layer.addLayer(vacuum_layer)
- multi_layer.addLayer(substrate_layer)
- return multi_layer
+
+ # Define materials
+ material_Particle = ba.HomogeneousMaterial("Particle", 3e-05, 2e-08)
+ material_Substrate = ba.HomogeneousMaterial("Substrate", 6e-06, 2e-08)
+ material_Vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
+
+ # Define form factors
+ ff = ba.FormFactorTruncatedSphere(5.0*nm, 5.0*nm, 0.0*nm)
+
+ # Define particles
+ particle = ba.Particle(material_Particle, ff)
+
+ # Define 2D lattices
+ lattice = ba.BasicLattice2D(10.0*nm, 10.0*nm, 90.0*deg, 0.0*deg)
+
+ # Define interference functions
+ iff = ba.InterferenceFunction2DLattice(lattice)
+ iff_pdf = ba.FTDecayFunction2DCauchy(100.0*nm, 100.0*nm, 0.0*deg)
+ iff.setDecayFunction(iff_pdf)
+
+ # Define particle layouts
+ layout = ba.ParticleLayout()
+ layout.addParticle(particle, 1.0)
+ layout.setInterferenceFunction(iff)
+ layout.setWeight(1)
+ layout.setTotalParticleSurfaceDensity(0.01)
+
+ # Define layers
+ layer_1 = ba.Layer(material_Vacuum)
+ layer_1.setNumberOfSlices(10)
+ layer_1.addLayout(layout)
+ layer_2 = ba.Layer(material_Substrate)
+
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer_1)
+ sample.addLayer(layer_2)
+
+ return sample
def get_simulation():
diff --git a/Examples/Python/sim05_Magnetism/MagneticSpheres.py b/Examples/Python/sim05_Magnetism/MagneticSpheres.py
index e6173541e6032358a1b351c0795b51eb5a3e8f19..1b6277ea0a801dfb4a98908c95612aabf6fd4a44 100644
--- a/Examples/Python/sim05_Magnetism/MagneticSpheres.py
+++ b/Examples/Python/sim05_Magnetism/MagneticSpheres.py
@@ -2,7 +2,7 @@
Simulation demo: magnetic spheres in substrate
"""
import bornagain as ba
-from bornagain import deg, angstrom, nm
+from bornagain import angstrom, deg, nm, nm2, kvector_t
# Magnetization of the particle's material (A/m)
magnetization_particle = ba.kvector_t(0.0, 0.0, 1e7)
@@ -12,30 +12,39 @@ def get_sample():
"""
Returns a sample with magnetic spheres in the substrate.
"""
- # defining materials
- particle_material = ba.HomogeneousMaterial("Particle", 2e-5, 4e-7,
- magnetization_particle)
- vacuum_material = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
- substrate_material = ba.HomogeneousMaterial("Substrate", 7e-6, 1.8e-7)
-
- # spherical magnetic particle
- sphere_ff = ba.FormFactorFullSphere(5*nm)
- sphere = ba.Particle(particle_material, sphere_ff)
- position = ba.kvector_t(0.0, 0.0, -10.0*nm)
- particle_layout = ba.ParticleLayout()
- particle_layout.addParticle(sphere, 1.0, position)
-
- # defining layers
- vacuum_layer = ba.Layer(vacuum_material)
- substrate_layer = ba.Layer(substrate_material)
- substrate_layer.addLayout(particle_layout)
-
- # defining the multilayer
- multi_layer = ba.MultiLayer()
- multi_layer.addLayer(vacuum_layer)
- multi_layer.addLayer(substrate_layer)
-
- return multi_layer
+
+ # Define materials
+ magnetic_field = kvector_t(0, 0, 10000000)
+ material_Particle = ba.HomogeneousMaterial("Particle", 2e-05, 4e-07,
+ magnetic_field)
+ material_Substrate = ba.HomogeneousMaterial("Substrate", 7e-06, 1.8e-07)
+ material_Vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
+
+ # Define form factors
+ ff = ba.FormFactorFullSphere(5.0*nm)
+
+ # Define particles
+ particle = ba.Particle(material_Particle, ff)
+ particle_position = kvector_t(0.0*nm, 0.0*nm, -10.0*nm)
+ particle.setPosition(particle_position)
+
+ # Define particle layouts
+ layout = ba.ParticleLayout()
+ layout.addParticle(particle, 1.0)
+ layout.setWeight(1)
+ layout.setTotalParticleSurfaceDensity(0.01)
+
+ # Define layers
+ layer_1 = ba.Layer(material_Vacuum)
+ layer_2 = ba.Layer(material_Substrate)
+ layer_2.addLayout(layout)
+
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer_1)
+ sample.addLayer(layer_2)
+
+ return sample
def get_simulation():
diff --git a/Examples/Python/sim11_Device/AxesInDifferentUnits.py b/Examples/Python/sim11_Device/AxesInDifferentUnits.py
index 684c7f176215f04d59aff02ab950e64cf5b81702..95eb46e47e4608d0c88c4b765ea8a801be03adf9 100644
--- a/Examples/Python/sim11_Device/AxesInDifferentUnits.py
+++ b/Examples/Python/sim11_Device/AxesInDifferentUnits.py
@@ -3,7 +3,7 @@ In this example we demonstrate how to plot simulation results with
axes in different units (nbins, mm, degs and QyQz).
"""
import bornagain as ba
-from bornagain import deg, angstrom, nm
+from bornagain import angstrom, deg, nm, nm2, kvector_t
from matplotlib import pyplot as plt
from matplotlib import rcParams
@@ -12,25 +12,35 @@ def get_sample():
"""
Returns a sample with uncorrelated cylinders on a substrate.
"""
- # defining materials
- m_vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
- m_substrate = ba.HomogeneousMaterial("Substrate", 6e-6, 2e-8)
- m_particle = ba.HomogeneousMaterial("Particle", 6e-4, 2e-8)
-
- # collection of particles
- cylinder_ff = ba.FormFactorCylinder(5*nm, 5*nm)
- cylinder = ba.Particle(m_particle, cylinder_ff)
- particle_layout = ba.ParticleLayout()
- particle_layout.addParticle(cylinder, 1.0)
-
- vacuum_layer = ba.Layer(m_vacuum)
- vacuum_layer.addLayout(particle_layout)
- substrate_layer = ba.Layer(m_substrate)
-
- multi_layer = ba.MultiLayer()
- multi_layer.addLayer(vacuum_layer)
- multi_layer.addLayer(substrate_layer)
- return multi_layer
+
+ # Define materials
+ material_Particle = ba.HomogeneousMaterial("Particle", 0.0006, 2e-08)
+ material_Substrate = ba.HomogeneousMaterial("Substrate", 6e-06, 2e-08)
+ material_Vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
+
+ # Define form factors
+ ff = ba.FormFactorCylinder(5.0*nm, 5.0*nm)
+
+ # Define particles
+ particle = ba.Particle(material_Particle, ff)
+
+ # Define particle layouts
+ layout = ba.ParticleLayout()
+ layout.addParticle(particle, 1.0)
+ layout.setWeight(1)
+ layout.setTotalParticleSurfaceDensity(0.01)
+
+ # Define layers
+ layer_1 = ba.Layer(material_Vacuum)
+ layer_1.addLayout(layout)
+ layer_2 = ba.Layer(material_Substrate)
+
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer_1)
+ sample.addLayer(layer_2)
+
+ return sample
def get_rectangular_detector():
diff --git a/Examples/Python/sim11_Device/BeamDivergence.py b/Examples/Python/sim11_Device/BeamDivergence.py
index d3bd4386de5798ffe65fcbdf202d102f2a2dfe07..9fdc00c06f875932142dd7e433719193725f30e0 100644
--- a/Examples/Python/sim11_Device/BeamDivergence.py
+++ b/Examples/Python/sim11_Device/BeamDivergence.py
@@ -2,33 +2,42 @@
Cylinder form factor in DWBA with beam divergence
"""
import bornagain as ba
-from bornagain import deg, angstrom, nm
+from bornagain import angstrom, deg, nm, nm2, kvector_t
def get_sample():
"""
Returns a sample with uncorrelated cylinders on a substrate.
"""
- # defining materials
- m_vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
- m_substrate = ba.HomogeneousMaterial("Substrate", 6e-6, 2e-8)
- m_particle = ba.HomogeneousMaterial("Particle", 6e-4, 2e-8)
- # collection of particles
- cylinder_ff = ba.FormFactorCylinder(5*nm, 5*nm)
- cylinder = ba.Particle(m_particle, cylinder_ff)
- particle_layout = ba.ParticleLayout()
- particle_layout.addParticle(cylinder, 1.0)
+ # Define materials
+ material_Particle = ba.HomogeneousMaterial("Particle", 0.0006, 2e-08)
+ material_Substrate = ba.HomogeneousMaterial("Substrate", 6e-06, 2e-08)
+ material_Vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
- # assembling the sample
- vacuum_layer = ba.Layer(m_vacuum)
- vacuum_layer.addLayout(particle_layout)
- substrate_layer = ba.Layer(m_substrate)
+ # Define form factors
+ ff = ba.FormFactorCylinder(5.0*nm, 5.0*nm)
- multi_layer = ba.MultiLayer()
- multi_layer.addLayer(vacuum_layer)
- multi_layer.addLayer(substrate_layer)
- return multi_layer
+ # Define particles
+ particle = ba.Particle(material_Particle, ff)
+
+ # Define particle layouts
+ layout = ba.ParticleLayout()
+ layout.addParticle(particle, 1.0)
+ layout.setWeight(1)
+ layout.setTotalParticleSurfaceDensity(0.01)
+
+ # Define layers
+ layer_1 = ba.Layer(material_Vacuum)
+ layer_1.addLayout(layout)
+ layer_2 = ba.Layer(material_Substrate)
+
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer_1)
+ sample.addLayer(layer_2)
+
+ return sample
def get_simulation():
diff --git a/Examples/Python/sim11_Device/ConstantBackground.py b/Examples/Python/sim11_Device/ConstantBackground.py
index 72d4236ca93d96cdd9e7138ecb612f205e33e4f5..75935b4b9c9b04a5bcd5f0a08c49f783e4219943 100644
--- a/Examples/Python/sim11_Device/ConstantBackground.py
+++ b/Examples/Python/sim11_Device/ConstantBackground.py
@@ -2,33 +2,42 @@
Cylinder form factor in DWBA with constant background
"""
import bornagain as ba
-from bornagain import deg, angstrom, nm
+from bornagain import angstrom, deg, nm, nm2, kvector_t
def get_sample():
"""
Returns a sample with uncorrelated cylinders on a substrate.
"""
- # defining materials
- m_vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
- m_substrate = ba.HomogeneousMaterial("Substrate", 6e-6, 2e-8)
- m_particle = ba.HomogeneousMaterial("Particle", 6e-4, 2e-8)
-
- # collection of particles
- cylinder_ff = ba.FormFactorCylinder(5*nm, 5*nm)
- cylinder = ba.Particle(m_particle, cylinder_ff)
- particle_layout = ba.ParticleLayout()
- particle_layout.addParticle(cylinder, 1.0)
-
- # assembling the sample
- vacuum_layer = ba.Layer(m_vacuum)
- vacuum_layer.addLayout(particle_layout)
- substrate_layer = ba.Layer(m_substrate)
-
- multi_layer = ba.MultiLayer()
- multi_layer.addLayer(vacuum_layer)
- multi_layer.addLayer(substrate_layer)
- return multi_layer
+
+ # Define materials
+ material_Particle = ba.HomogeneousMaterial("Particle", 0.0006, 2e-08)
+ material_Substrate = ba.HomogeneousMaterial("Substrate", 6e-06, 2e-08)
+ material_Vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
+
+ # Define form factors
+ ff = ba.FormFactorCylinder(5.0*nm, 5.0*nm)
+
+ # Define particles
+ particle = ba.Particle(material_Particle, ff)
+
+ # Define particle layouts
+ layout = ba.ParticleLayout()
+ layout.addParticle(particle, 1.0)
+ layout.setWeight(1)
+ layout.setTotalParticleSurfaceDensity(0.01)
+
+ # Define layers
+ layer_1 = ba.Layer(material_Vacuum)
+ layer_1.addLayout(layout)
+ layer_2 = ba.Layer(material_Substrate)
+
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer_1)
+ sample.addLayer(layer_2)
+
+ return sample
def get_simulation():
diff --git a/Examples/Python/sim11_Device/DetectorResolutionFunction.py b/Examples/Python/sim11_Device/DetectorResolutionFunction.py
index 9b06b4de856731fee54104a25ca276f8768e5e9e..5add835f7998d4aad0f4c544aee6734061adfb37 100644
--- a/Examples/Python/sim11_Device/DetectorResolutionFunction.py
+++ b/Examples/Python/sim11_Device/DetectorResolutionFunction.py
@@ -2,33 +2,42 @@
Cylinder form factor in DWBA with detector resolution function applied
"""
import bornagain as ba
-from bornagain import deg, angstrom, nm
+from bornagain import angstrom, deg, nm, nm2, kvector_t
def get_sample():
"""
Returns a sample with uncorrelated cylinders on a substrate.
"""
- # defining materials
- m_vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
- m_substrate = ba.HomogeneousMaterial("Substrate", 6e-6, 2e-8)
- m_particle = ba.HomogeneousMaterial("Particle", 6e-4, 2e-8)
-
- # collection of particles
- cylinder_ff = ba.FormFactorCylinder(5*nm, 5*nm)
- cylinder = ba.Particle(m_particle, cylinder_ff)
- particle_layout = ba.ParticleLayout()
- particle_layout.addParticle(cylinder, 1.0)
-
- # assembling the sample
- vacuum_layer = ba.Layer(m_vacuum)
- vacuum_layer.addLayout(particle_layout)
- substrate_layer = ba.Layer(m_substrate)
-
- multi_layer = ba.MultiLayer()
- multi_layer.addLayer(vacuum_layer)
- multi_layer.addLayer(substrate_layer)
- return multi_layer
+
+ # Define materials
+ material_Particle = ba.HomogeneousMaterial("Particle", 0.0006, 2e-08)
+ material_Substrate = ba.HomogeneousMaterial("Substrate", 6e-06, 2e-08)
+ material_Vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
+
+ # Define form factors
+ ff = ba.FormFactorCylinder(5.0*nm, 5.0*nm)
+
+ # Define particles
+ particle = ba.Particle(material_Particle, ff)
+
+ # Define particle layouts
+ layout = ba.ParticleLayout()
+ layout.addParticle(particle, 1.0)
+ layout.setWeight(1)
+ layout.setTotalParticleSurfaceDensity(0.01)
+
+ # Define layers
+ layer_1 = ba.Layer(material_Vacuum)
+ layer_1.addLayout(layout)
+ layer_2 = ba.Layer(material_Substrate)
+
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer_1)
+ sample.addLayer(layer_2)
+
+ return sample
def get_simulation():
diff --git a/Examples/Python/sim11_Device/OffSpecularSimulation.py b/Examples/Python/sim11_Device/OffSpecularSimulation.py
index 825ece1c4ca575ffbf2de05ffd8cabea2af1d538..fae853ac2cf200101b68f245ea4b10bfbfcc603c 100644
--- a/Examples/Python/sim11_Device/OffSpecularSimulation.py
+++ b/Examples/Python/sim11_Device/OffSpecularSimulation.py
@@ -2,7 +2,7 @@
Long boxes at 1D lattice, ba.OffSpecular simulation
"""
import bornagain as ba
-from bornagain import deg, angstrom, nm
+from bornagain import angstrom, deg, nm, nm2, kvector_t
phi_f_min, phi_f_max = -1.0, 1.0
alpha_f_min, alpha_f_max = 0.0, 10.0
@@ -14,36 +14,43 @@ def get_sample():
Returns a sample with a grating on a substrate,
modelled by infinitely long boxes forming a 1D lattice.
"""
- # defining materials
- m_vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
- m_substrate = ba.HomogeneousMaterial("Substrate", 6e-6, 2e-8)
- m_particle = ba.HomogeneousMaterial("Particle", 6e-4, 2e-8)
-
- # collection of particles
- lattice_length = 100.0*nm
- lattice_rotation_angle = 0.0*deg
- interference = ba.InterferenceFunction1DLattice(lattice_length,
- lattice_rotation_angle)
- pdf = ba.FTDecayFunction1DCauchy(1e+6)
- interference.setDecayFunction(pdf)
-
- box_ff = ba.FormFactorBox(1000*nm, 20*nm, 10.0*nm)
- box = ba.Particle(m_particle, box_ff)
- transform = ba.RotationZ(90.0*deg)
- particle_layout = ba.ParticleLayout()
- particle_layout.addParticle(box, 1.0, ba.kvector_t(0.0, 0.0, 0.0),
- transform)
- particle_layout.setInterferenceFunction(interference)
-
- # assembling the sample
- vacuum_layer = ba.Layer(m_vacuum)
- vacuum_layer.addLayout(particle_layout)
- substrate_layer = ba.Layer(m_substrate)
-
- multi_layer = ba.MultiLayer()
- multi_layer.addLayer(vacuum_layer)
- multi_layer.addLayer(substrate_layer)
- return multi_layer
+
+ # Define materials
+ material_Particle = ba.HomogeneousMaterial("Particle", 0.0006, 2e-08)
+ material_Substrate = ba.HomogeneousMaterial("Substrate", 6e-06, 2e-08)
+ material_Vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
+
+ # Define form factors
+ ff = ba.FormFactorBox(1000.0*nm, 20.0*nm, 10.0*nm)
+
+ # Define particles
+ particle = ba.Particle(material_Particle, ff)
+ particle_rotation = ba.RotationZ(90.0*deg)
+ particle.setRotation(particle_rotation)
+
+ # Define interference functions
+ iff = ba.InterferenceFunction1DLattice(100.0*nm, 0.0*deg)
+ iff_pdf = ba.FTDecayFunction1DCauchy(1000000.0*nm)
+ iff.setDecayFunction(iff_pdf)
+
+ # Define particle layouts
+ layout = ba.ParticleLayout()
+ layout.addParticle(particle, 1.0)
+ layout.setInterferenceFunction(iff)
+ layout.setWeight(1)
+ layout.setTotalParticleSurfaceDensity(0.01)
+
+ # Define layers
+ layer_1 = ba.Layer(material_Vacuum)
+ layer_1.addLayout(layout)
+ layer_2 = ba.Layer(material_Substrate)
+
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer_1)
+ sample.addLayer(layer_2)
+
+ return sample
def get_simulation():
diff --git a/Examples/Python/sim11_Device/RectangularDetector.py b/Examples/Python/sim11_Device/RectangularDetector.py
index 97891089d3c6b71d1a8d45ea06188c56ed6cdfd4..cabd270157fa3e12a2b23ae0b7e9920dbe2968df 100644
--- a/Examples/Python/sim11_Device/RectangularDetector.py
+++ b/Examples/Python/sim11_Device/RectangularDetector.py
@@ -4,7 +4,7 @@ Results will be compared against simulation with spherical detector.
"""
import numpy
import bornagain as ba
-from bornagain import deg, angstrom, nm
+from bornagain import angstrom, deg, nm, nm2, kvector_t
import matplotlib
from matplotlib import pyplot as plt
@@ -17,26 +17,35 @@ def get_sample():
"""
Returns a sample with cylindrical particles on a substrate.
"""
- # defining materials
- m_vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
- m_substrate = ba.HomogeneousMaterial("Substrate", 6e-6, 2e-8)
- m_particle = ba.HomogeneousMaterial("Particle", 6e-4, 2e-8)
-
- # collection of particles
- edge = 40*nm
- ff = ba.FormFactorBox(edge, edge, edge)
- cylinder = ba.Particle(m_particle, ff)
- particle_layout = ba.ParticleLayout()
- particle_layout.addParticle(cylinder, 1.0)
-
- vacuum_layer = ba.Layer(m_vacuum)
- vacuum_layer.addLayout(particle_layout)
- substrate_layer = ba.Layer(m_substrate)
-
- multi_layer = ba.MultiLayer()
- multi_layer.addLayer(vacuum_layer)
- multi_layer.addLayer(substrate_layer)
- return multi_layer
+
+ # Define materials
+ material_Particle = ba.HomogeneousMaterial("Particle", 0.0006, 2e-08)
+ material_Substrate = ba.HomogeneousMaterial("Substrate", 6e-06, 2e-08)
+ material_Vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
+
+ # Define form factors
+ ff = ba.FormFactorBox(40.0*nm, 40.0*nm, 40.0*nm)
+
+ # Define particles
+ particle = ba.Particle(material_Particle, ff)
+
+ # Define particle layouts
+ layout = ba.ParticleLayout()
+ layout.addParticle(particle, 1.0)
+ layout.setWeight(1)
+ layout.setTotalParticleSurfaceDensity(0.01)
+
+ # Define layers
+ layer_1 = ba.Layer(material_Vacuum)
+ layer_1.addLayout(layout)
+ layer_2 = ba.Layer(material_Substrate)
+
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer_1)
+ sample.addLayer(layer_2)
+
+ return sample
def get_spherical_detector():
diff --git a/Examples/Python/sim21_Reflectometry/BasicPolarizedReflectometry.py b/Examples/Python/sim21_Reflectometry/BasicPolarizedReflectometry.py
index 646c07d16f9ca6b611d2a36012146b15fa86504a..a61b8f9da614369882137b1e201b4eba94cbb7da 100644
--- a/Examples/Python/sim21_Reflectometry/BasicPolarizedReflectometry.py
+++ b/Examples/Python/sim21_Reflectometry/BasicPolarizedReflectometry.py
@@ -4,7 +4,7 @@ magnetized sample.
"""
import bornagain as ba
-from bornagain import deg, angstrom
+from bornagain import angstrom, deg, nm, nm2, kvector_t
import matplotlib.pyplot as plt
@@ -14,24 +14,24 @@ def get_sample():
Defines sample and returns it
"""
- # creating materials
- m_ambient = ba.MaterialBySLD("Ambient", 0.0, 0.0)
- m_layer_mat = ba.MaterialBySLD("Layer", 1e-4, 1e-8,
- ba.kvector_t(0.0, 1e8, 0.0))
- m_substrate = ba.MaterialBySLD("Substrate", 7e-5, 2e-6)
+ # Define materials
+ material_Ambient = ba.MaterialBySLD("Ambient", 0.0, 0.0)
+ magnetic_field = kvector_t(0, 100000000, 0)
+ material_Layer = ba.MaterialBySLD("Layer", 0.0001, 1e-08, magnetic_field)
+ material_Substrate = ba.MaterialBySLD("Substrate", 7e-05, 2e-06)
- # creating layers
- ambient_layer = ba.Layer(m_ambient)
- layer = ba.Layer(m_layer_mat, 10)
- substrate_layer = ba.Layer(m_substrate)
+ # Define layers
+ layer_1 = ba.Layer(material_Ambient)
+ layer_2 = ba.Layer(material_Layer, 10.0*nm)
+ layer_3 = ba.Layer(material_Substrate)
- # creating multilayer
- multi_layer = ba.MultiLayer()
- multi_layer.addLayer(ambient_layer)
- multi_layer.addLayer(layer)
- multi_layer.addLayer(substrate_layer)
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer_1)
+ sample.addLayer(layer_2)
+ sample.addLayer(layer_3)
- return multi_layer
+ return sample
def get_simulation(scan_size=500):
diff --git a/Examples/Python/sim21_Reflectometry/PolarizedNoAnalyzer.py b/Examples/Python/sim21_Reflectometry/PolarizedNoAnalyzer.py
index 789e705a95ae5cecfaff682bbde8583ecb86b9fe..3be7294ae9257c39797f09005a7d29121cd65613 100644
--- a/Examples/Python/sim21_Reflectometry/PolarizedNoAnalyzer.py
+++ b/Examples/Python/sim21_Reflectometry/PolarizedNoAnalyzer.py
@@ -6,7 +6,7 @@ import numpy
import matplotlib.pyplot as plt
import bornagain as ba
-from bornagain import deg, angstrom
+from bornagain import angstrom, deg, nm, nm2, kvector_t
def get_sample():
@@ -14,29 +14,24 @@ def get_sample():
Defines sample and returns it
"""
- # creating materials
- magnetizationMagnitude = 1e8
- angle = 30*deg
- magnetizationVector = ba.kvector_t(magnetizationMagnitude*numpy.sin(angle),
- magnetizationMagnitude*numpy.cos(angle),
- 0)
-
- m_ambient = ba.MaterialBySLD("Ambient", 0.0, 0.0)
- m_layer_mat = ba.MaterialBySLD("Layer", 1e-4, 1e-8, magnetizationVector)
- m_substrate = ba.MaterialBySLD("Substrate", 7e-5, 2e-6)
-
- # creating layers
- ambient_layer = ba.Layer(m_ambient)
- layer = ba.Layer(m_layer_mat, 10)
- substrate_layer = ba.Layer(m_substrate)
-
- # creating multilayer
- multi_layer = ba.MultiLayer()
- multi_layer.addLayer(ambient_layer)
- multi_layer.addLayer(layer)
- multi_layer.addLayer(substrate_layer)
-
- return multi_layer
+ # Define materials
+ material_Ambient = ba.MaterialBySLD("Ambient", 0.0, 0.0)
+ magnetic_field = kvector_t(50000000, 86602540.3784, 0)
+ material_Layer = ba.MaterialBySLD("Layer", 0.0001, 1e-08, magnetic_field)
+ material_Substrate = ba.MaterialBySLD("Substrate", 7e-05, 2e-06)
+
+ # Define layers
+ layer_1 = ba.Layer(material_Ambient)
+ layer_2 = ba.Layer(material_Layer, 10.0*nm)
+ layer_3 = ba.Layer(material_Substrate)
+
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer_1)
+ sample.addLayer(layer_2)
+ sample.addLayer(layer_3)
+
+ return sample
def get_simulation(scan_size=500):
diff --git a/Examples/Python/sim21_Reflectometry/PolarizedSpinFlip.py b/Examples/Python/sim21_Reflectometry/PolarizedSpinFlip.py
index ae84979c89406edf2d92a4c0037fc6742170281a..901087ebb03e3a7cdba54d300c2bdb7811e9a2e2 100644
--- a/Examples/Python/sim21_Reflectometry/PolarizedSpinFlip.py
+++ b/Examples/Python/sim21_Reflectometry/PolarizedSpinFlip.py
@@ -5,7 +5,7 @@ a magnetized sample.
import numpy
import bornagain as ba
-from bornagain import deg, angstrom
+from bornagain import angstrom, deg, nm, nm2, kvector_t
import matplotlib.pyplot as plt
@@ -15,30 +15,24 @@ def get_sample():
Defines sample and returns it
"""
- # parametrize the magnetization
- magnetizationMagnitude = 1e8
- angle = 30*deg
- magnetizationVector = ba.kvector_t(magnetizationMagnitude*numpy.sin(angle),
- magnetizationMagnitude*numpy.cos(angle),
- 0)
-
- # creating materials
- m_ambient = ba.MaterialBySLD("Ambient", 0.0, 0.0)
- m_layer_mat = ba.MaterialBySLD("Layer", 1e-4, 1e-8, magnetizationVector)
- m_substrate = ba.MaterialBySLD("Substrate", 7e-5, 2e-6)
-
- # creating layers
- ambient_layer = ba.Layer(m_ambient)
- layer = ba.Layer(m_layer_mat, 10)
- substrate_layer = ba.Layer(m_substrate)
-
- # creating multilayer
- multi_layer = ba.MultiLayer()
- multi_layer.addLayer(ambient_layer)
- multi_layer.addLayer(layer)
- multi_layer.addLayer(substrate_layer)
-
- return multi_layer
+ # Define materials
+ material_Ambient = ba.MaterialBySLD("Ambient", 0.0, 0.0)
+ magnetic_field = kvector_t(50000000, 86602540.3784, 0)
+ material_Layer = ba.MaterialBySLD("Layer", 0.0001, 1e-08, magnetic_field)
+ material_Substrate = ba.MaterialBySLD("Substrate", 7e-05, 2e-06)
+
+ # Define layers
+ layer_1 = ba.Layer(material_Ambient)
+ layer_2 = ba.Layer(material_Layer, 10.0*nm)
+ layer_3 = ba.Layer(material_Substrate)
+
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer_1)
+ sample.addLayer(layer_2)
+ sample.addLayer(layer_3)
+
+ return sample
def get_simulation(scan_size=500):
diff --git a/Examples/Python/sim22_OffSpecular/BoxesWithSpecularPeak.py b/Examples/Python/sim22_OffSpecular/BoxesWithSpecularPeak.py
index 7a770a1f547bc6f665419b75a89293aa7f74d377..e5d58ff02734e888aacc75f51f14792d2d29278d 100644
--- a/Examples/Python/sim22_OffSpecular/BoxesWithSpecularPeak.py
+++ b/Examples/Python/sim22_OffSpecular/BoxesWithSpecularPeak.py
@@ -2,39 +2,51 @@
Square lattice of boxes on substrate (including the specular peak)
"""
import bornagain as ba
-from bornagain import deg, angstrom, nm
+from bornagain import angstrom, deg, nm, nm2, kvector_t
def get_sample():
"""
Returns a sample with boxes on a substrate.
"""
- # defining materials
- m_vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
- m_substrate = ba.HomogeneousMaterial("Substrate", 6e-6, 2e-8)
- m_particle = ba.HomogeneousMaterial("Particle", 3e-5, 2e-8)
-
- # cylindrical particle
- box_ff = ba.FormFactorBox(5*nm, 5*nm, 10*nm)
- box = ba.Particle(m_particle, box_ff)
- particle_layout = ba.ParticleLayout()
- particle_layout.addParticle(box)
-
- # interference function
- interference = ba.InterferenceFunction2DLattice(
- ba.SquareLattice2D(8*nm, 0*deg))
- pdf = ba.FTDecayFunction2DCauchy(100*nm, 100*nm, 0)
- interference.setDecayFunction(pdf)
- particle_layout.setInterferenceFunction(interference)
-
- vacuum_layer = ba.Layer(m_vacuum)
- vacuum_layer.addLayout(particle_layout)
- substrate_layer = ba.Layer(m_substrate)
-
- multi_layer = ba.MultiLayer()
- multi_layer.addLayer(vacuum_layer)
- multi_layer.addLayer(substrate_layer)
- return multi_layer
+
+ # Define materials
+ material_Particle = ba.HomogeneousMaterial("Particle", 3e-05, 2e-08)
+ material_Substrate = ba.HomogeneousMaterial("Substrate", 6e-06, 2e-08)
+ material_Vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
+
+ # Define form factors
+ ff = ba.FormFactorBox(5.0*nm, 5.0*nm, 10.0*nm)
+
+ # Define particles
+ particle = ba.Particle(material_Particle, ff)
+
+ # Define 2D lattices
+ lattice = ba.BasicLattice2D(8.0*nm, 8.0*nm, 90.0*deg, 0.0*deg)
+
+ # Define interference functions
+ iff = ba.InterferenceFunction2DLattice(lattice)
+ iff_pdf = ba.FTDecayFunction2DCauchy(100.0*nm, 100.0*nm, 0.0*deg)
+ iff.setDecayFunction(iff_pdf)
+
+ # Define particle layouts
+ layout = ba.ParticleLayout()
+ layout.addParticle(particle, 1.0)
+ layout.setInterferenceFunction(iff)
+ layout.setWeight(1)
+ layout.setTotalParticleSurfaceDensity(0.015625)
+
+ # Define layers
+ layer_1 = ba.Layer(material_Vacuum)
+ layer_1.addLayout(layout)
+ layer_2 = ba.Layer(material_Substrate)
+
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer_1)
+ sample.addLayer(layer_2)
+
+ return sample
def get_simulation():
diff --git a/Examples/Python/sim23_SAS/PolarizedSANS.py b/Examples/Python/sim23_SAS/PolarizedSANS.py
index 3728621a111b0faae8ea7e01f1dce245b905daeb..eff2352864e67fb3f72c99bb0de61adebcfe60f9 100644
--- a/Examples/Python/sim23_SAS/PolarizedSANS.py
+++ b/Examples/Python/sim23_SAS/PolarizedSANS.py
@@ -4,7 +4,7 @@ simulated with BornAgain.
"""
import bornagain as ba
-from bornagain import deg, nm, kvector_t
+from bornagain import angstrom, deg, nm, nm2, kvector_t
# Magnetization of the particle's core material (A/m)
magnetization_core = kvector_t(0.0, 0.0, 1e7)
@@ -14,31 +14,41 @@ def get_sample():
"""
Returns a sample with a magnetic core-shell particle in a solvent.
"""
- # Defining Materials
- mat_solvent = ba.HomogeneousMaterial("Solvent", 5e-6, 0.0)
- mat_core = ba.HomogeneousMaterial("Core", 6e-6, 2e-8, magnetization_core)
- mat_shell = ba.HomogeneousMaterial("Shell", 1e-7, 2e-8)
- # Defining Layer
- solvent_layer = ba.Layer(mat_solvent)
+ # Define materials
+ magnetic_field = kvector_t(0, 0, 10000000)
+ material_Core = ba.HomogeneousMaterial("Core", 6e-06, 2e-08, magnetic_field)
+ material_Shell = ba.HomogeneousMaterial("Shell", 1e-07, 2e-08)
+ material_Solvent = ba.HomogeneousMaterial("Solvent", 5e-06, 0.0)
- # Defining particle layout with a core-shell particle
+ # Define form factors
+ ff_1 = ba.FormFactorFullSphere(10.0*nm)
+ ff_2 = ba.FormFactorFullSphere(12.0*nm)
+
+ # Define particles
+ particle_1 = ba.Particle(material_Core, ff_1)
+ particle_1_position = kvector_t(0.0*nm, 0.0*nm, 2.0*nm)
+ particle_1.setPosition(particle_1_position)
+ particle_2 = ba.Particle(material_Shell, ff_2)
+
+ # Define core shell particles
+ particle_3 = ba.ParticleCoreShell(particle_2, particle_1)
+
+ # Define particle layouts
layout = ba.ParticleLayout()
- core_sphere_ff = ba.FormFactorFullSphere(10*nm)
- shell_sphere_ff = ba.FormFactorFullSphere(12*nm)
- core = ba.Particle(mat_core, core_sphere_ff)
- shell = ba.Particle(mat_shell, shell_sphere_ff)
- position = kvector_t(0.0, 0.0, 2.0)
- particleCoreShell = ba.ParticleCoreShell(shell, core, position)
- layout.addParticle(particleCoreShell)
-
- # Adding layout to layer
- solvent_layer.addLayout(layout)
-
- # Defining Multilayer with single layer
- multiLayer = ba.MultiLayer()
- multiLayer.addLayer(solvent_layer)
- return multiLayer
+ layout.addParticle(particle_3, 1.0)
+ layout.setWeight(1)
+ layout.setTotalParticleSurfaceDensity(0.01)
+
+ # Define layers
+ layer = ba.Layer(material_Solvent)
+ layer.addLayout(layout)
+
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer)
+
+ return sample
def get_simulation():
diff --git a/Examples/Python/sim29_DepthProbe/DepthProbe.py b/Examples/Python/sim29_DepthProbe/DepthProbe.py
index a27ef8e6b33217183d5f58949172ec9483553ba8..22de45ab8dae4c9be84d3f70d54bbcf8cb0b22af 100644
--- a/Examples/Python/sim29_DepthProbe/DepthProbe.py
+++ b/Examples/Python/sim29_DepthProbe/DepthProbe.py
@@ -24,7 +24,7 @@ to the sample, z = 0 corresponding to the sample
surface
"""
import bornagain as ba
-from bornagain import deg, angstrom, nm
+from bornagain import angstrom, deg, nm, nm2, kvector_t
# layer thicknesses in angstroms
t_Ti = 130.0*angstrom
@@ -55,31 +55,29 @@ def get_sample():
Constructs a sample with one resonating Ti/Pt layer
"""
- # define materials
- m_Si = ba.HomogeneousMaterial("Si", 3.3009e-05, 0.0)
- m_Ti = ba.HomogeneousMaterial("Ti", -3.0637e-05, 1.5278e-08)
- m_TiO2 = ba.HomogeneousMaterial("TiO2", 4.1921e-05, 8.1293e-09)
- m_Pt = ba.HomogeneousMaterial("Pt", 1.0117e-04, 3.01822e-08)
- m_D2O = ba.HomogeneousMaterial("D2O", 1.0116e-04, 1.8090e-12)
-
- # create layers
- l_Si = ba.Layer(m_Si)
- l_Ti = ba.Layer(m_Ti, 130.0*angstrom)
- l_Pt = ba.Layer(m_Pt, 320.0*angstrom)
- l_Ti_top = ba.Layer(m_Ti, 100.0*angstrom)
- l_TiO2 = ba.Layer(m_TiO2, 30.0*angstrom)
- l_D2O = ba.Layer(m_D2O)
-
- # construct sample
+ # Define materials
+ material_D2O = ba.HomogeneousMaterial("D2O", 0.00010116, 1.809e-12)
+ material_Pt = ba.HomogeneousMaterial("Pt", 0.00010117, 3.01822e-08)
+ material_Si = ba.HomogeneousMaterial("Si", 3.3009e-05, 0.0)
+ material_Ti = ba.HomogeneousMaterial("Ti", -3.0637e-05, 1.5278e-08)
+ material_TiO2 = ba.HomogeneousMaterial("TiO2", 4.1921e-05, 8.1293e-09)
+
+ # Define layers
+ layer_1 = ba.Layer(material_Si)
+ layer_2 = ba.Layer(material_Ti, 13.0*nm)
+ layer_3 = ba.Layer(material_Pt, 32.0*nm)
+ layer_4 = ba.Layer(material_Ti, 10.0*nm)
+ layer_5 = ba.Layer(material_TiO2, 3.0*nm)
+ layer_6 = ba.Layer(material_D2O)
+
+ # Define sample
sample = ba.MultiLayer()
- sample.addLayer(l_Si)
-
- sample.addLayer(l_Ti)
- sample.addLayer(l_Pt)
-
- sample.addLayer(l_Ti_top)
- sample.addLayer(l_TiO2)
- sample.addLayer(l_D2O)
+ sample.addLayer(layer_1)
+ sample.addLayer(layer_2)
+ sample.addLayer(layer_3)
+ sample.addLayer(layer_4)
+ sample.addLayer(layer_5)
+ sample.addLayer(layer_6)
return sample
diff --git a/Examples/Python/sim31_Parameterization/AccessingSimulationResults.py b/Examples/Python/sim31_Parameterization/AccessingSimulationResults.py
index a443761cff9f7168a8f16093be9c2cf5ff957eb5..5dd69b16e2eebbd937a9ad7d7c6944ce0eb5346e 100644
--- a/Examples/Python/sim31_Parameterization/AccessingSimulationResults.py
+++ b/Examples/Python/sim31_Parameterization/AccessingSimulationResults.py
@@ -5,7 +5,7 @@ The standard "Cylinders in DWBA" sample is used to setup the simulation.
import math
import random
import bornagain as ba
-from bornagain import deg, angstrom, nm
+from bornagain import angstrom, deg, nm, nm2, kvector_t
from matplotlib import pyplot as plt
from matplotlib import rcParams
@@ -14,25 +14,35 @@ def get_sample():
"""
Returns a sample with uncorrelated cylinders on a substrate.
"""
- # defining materials
- m_vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
- m_substrate = ba.HomogeneousMaterial("Substrate", 6e-6, 2e-8)
- m_particle = ba.HomogeneousMaterial("Particle", 6e-4, 2e-8)
-
- # collection of particles
- cylinder_ff = ba.FormFactorCylinder(5*nm, 5*nm)
- cylinder = ba.Particle(m_particle, cylinder_ff)
- particle_layout = ba.ParticleLayout()
- particle_layout.addParticle(cylinder, 1.0)
-
- vacuum_layer = ba.Layer(m_vacuum)
- vacuum_layer.addLayout(particle_layout)
- substrate_layer = ba.Layer(m_substrate)
-
- multi_layer = ba.MultiLayer()
- multi_layer.addLayer(vacuum_layer)
- multi_layer.addLayer(substrate_layer)
- return multi_layer
+
+ # Define materials
+ material_Particle = ba.HomogeneousMaterial("Particle", 0.0006, 2e-08)
+ material_Substrate = ba.HomogeneousMaterial("Substrate", 6e-06, 2e-08)
+ material_Vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
+
+ # Define form factors
+ ff = ba.FormFactorCylinder(5.0*nm, 5.0*nm)
+
+ # Define particles
+ particle = ba.Particle(material_Particle, ff)
+
+ # Define particle layouts
+ layout = ba.ParticleLayout()
+ layout.addParticle(particle, 1.0)
+ layout.setWeight(1)
+ layout.setTotalParticleSurfaceDensity(0.01)
+
+ # Define layers
+ layer_1 = ba.Layer(material_Vacuum)
+ layer_1.addLayout(layout)
+ layer_2 = ba.Layer(material_Substrate)
+
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer_1)
+ sample.addLayer(layer_2)
+
+ return sample
def get_simulation():
diff --git a/Examples/Python/sim31_Parameterization/SimulationParameters.py b/Examples/Python/sim31_Parameterization/SimulationParameters.py
index 39481ad8338e6fdf91bc41dee35cebbc5c9088cf..195210a3e01c8bc8fc9c285a56b841766fbd7b7e 100644
--- a/Examples/Python/sim31_Parameterization/SimulationParameters.py
+++ b/Examples/Python/sim31_Parameterization/SimulationParameters.py
@@ -7,7 +7,7 @@ these parameters on the fly during runtime.
from __future__ import print_function
import bornagain as ba
-from bornagain import deg, angstrom, nm
+from bornagain import angstrom, deg, nm, nm2, kvector_t
def get_sample():
@@ -15,31 +15,42 @@ def get_sample():
Returns a sample with uncorrelated cylinders and prisms on a substrate.
Parameter set is fixed.
"""
- # defining materials
- m_vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
- m_substrate = ba.HomogeneousMaterial("Substrate", 6e-6, 2e-8)
- m_particle = ba.HomogeneousMaterial("Particle", 6e-4, 2e-8)
- # collection of particles
- cylinder_ff = ba.FormFactorCylinder(5*nm, 5*nm)
- cylinder = ba.Particle(m_particle, cylinder_ff)
- prism_ff = ba.FormFactorPrism3(5*nm, 5*nm)
- prism = ba.Particle(m_particle, prism_ff)
+ # Define materials
+ material_Particle = ba.HomogeneousMaterial("Particle", 0.0006, 2e-08)
+ material_Substrate = ba.HomogeneousMaterial("Substrate", 6e-06, 2e-08)
+ material_Vacuum = ba.HomogeneousMaterial("Vacuum", 0.0, 0.0)
+ # Define form factors
+ ff_1 = ba.FormFactorCylinder(5.0*nm, 5.0*nm)
+ ff_2 = ba.FormFactorPrism3(5.0*nm, 5.0*nm)
+
+ # Define particles
+ particle_1 = ba.Particle(material_Particle, ff_1)
+ particle_2 = ba.Particle(material_Particle, ff_2)
+
+ # Define interference functions
+ iff = ba.InterferenceFunctionNone()
+
+ # Define particle layouts
layout = ba.ParticleLayout()
- layout.addParticle(cylinder, 0.5)
- layout.addParticle(prism, 0.5)
- interference = ba.InterferenceFunctionNone()
- layout.setInterferenceFunction(interference)
-
- # vacuum layer with particles and substrate form multi layer
- vacuum_layer = ba.Layer(m_vacuum)
- vacuum_layer.addLayout(layout)
- substrate_layer = ba.Layer(m_substrate, 0)
- multi_layer = ba.MultiLayer()
- multi_layer.addLayer(vacuum_layer)
- multi_layer.addLayer(substrate_layer)
- return multi_layer
+ layout.addParticle(particle_1, 0.5)
+ layout.addParticle(particle_2, 0.5)
+ layout.setInterferenceFunction(iff)
+ layout.setWeight(1)
+ layout.setTotalParticleSurfaceDensity(0.01)
+
+ # Define layers
+ layer_1 = ba.Layer(material_Vacuum)
+ layer_1.addLayout(layout)
+ layer_2 = ba.Layer(material_Substrate)
+
+ # Define sample
+ sample = ba.MultiLayer()
+ sample.addLayer(layer_1)
+ sample.addLayer(layer_2)
+
+ return sample
def get_simulation():
diff --git a/Examples/cpp/modules/FindBornAgain.cmake b/Examples/cpp/modules/FindBornAgain.cmake
index 81253fd49a1b1fe7eccc6ab77601b2e4030323ea..e3cf00544e45560cfec863b047485353ec6e5e5c 100644
--- a/Examples/cpp/modules/FindBornAgain.cmake
+++ b/Examples/cpp/modules/FindBornAgain.cmake
@@ -17,22 +17,20 @@
set(BORNAGAINSYS $ENV{BORNAGAINSYS})
+set(CoreComponents "Base;Param;Sample;Device;Core")
if(BORNAGAINSYS)
- set(BORNAGAIN_LIBRARY_DIR ${BORNAGAINSYS}/lib/BornAgain-1.6)
- set(BORNAGAIN_INCLUDE_DIR ${BORNAGAINSYS}/include/BornAgain-1.6)
+ set(BORNAGAIN_LIBRARY_DIR ${BORNAGAINSYS}/lib/BornAgain-1.18)
+ set(BORNAGAIN_INCLUDE_DIR ${BORNAGAINSYS}/include/BornAgain-1.18)
endif()
-find_library (BORNAGAIN_CORE _libBornAgainCore.so
+foreach(lib ${CoreComponents})
+ message(STATUS ${lib})
+ find_library (BORNAGAIN_${lib} _libBornAgain${lib}.so
PATHS ${BORNAGAIN_LIBRARY_DIR}
- HINTS ${BORNAGAIN_LIBRARY_DIR}
-)
-
-find_library (BORNAGAIN_FIT _libBornAgainFit.so
- PATHS ${BORNAGAIN_LIBRARY_DIR}
- HINTS ${BORNAGAIN_LIBRARY_DIR}
-)
-set(BORNAGAIN_LIBRARIES ${BORNAGAIN_CORE} ${BORNAGAIN_FIT})
+ HINTS ${BORNAGAIN_LIBRARY_DIR})
+ list(APPEND BORNAGAIN_LIBRARIES ${BORNAGAIN_${lib}})
+endforeach()
find_path(BORNAGAIN_INCLUDE_DIR BAVersion.h
PATHS /usr/include /usr/local/include /opt/local/include ${BORNAGAIN_INCLUDE_DIR}
diff --git a/Sample/Fresnel/MatrixFresnelMap.cpp b/Sample/Fresnel/MatrixFresnelMap.cpp
index 4004bbbf5c0d23a1d65ec92c23702c3169b8fa10..f1ba1dc45033bff2f220a1152466606e35fb6f3d 100644
--- a/Sample/Fresnel/MatrixFresnelMap.cpp
+++ b/Sample/Fresnel/MatrixFresnelMap.cpp
@@ -13,9 +13,9 @@
// ************************************************************************************************
#include "Sample/Fresnel/MatrixFresnelMap.h"
+#include "Sample/RT/ILayerRTCoefficients.h"
#include "Base/Pixel/SimulationElement.h"
#include "Sample/Slice/Slice.h"
-#include "Sample/Specular/SpecularMagneticOldStrategy.h"
#include <functional>
MatrixFresnelMap::MatrixFresnelMap(std::unique_ptr<ISpecularStrategy> strategy)
diff --git a/Sample/Fresnel/MatrixFresnelMap.h b/Sample/Fresnel/MatrixFresnelMap.h
index 2616af9f5b8fb1915c51c999054a99a0816659f5..1a669e4879fb43ead39ef9dbe38c68a4e95ea622 100644
--- a/Sample/Fresnel/MatrixFresnelMap.h
+++ b/Sample/Fresnel/MatrixFresnelMap.h
@@ -16,8 +16,6 @@
#define BORNAGAIN_SAMPLE_FRESNEL_MATRIXFRESNELMAP_H
#include "Sample/Fresnel/IFresnelMap.h"
-#include "Sample/RT/MatrixRTCoefficients.h"
-#include "Sample/Specular/SpecularMagneticStrategy.h"
#include <cstddef>
#include <memory>
#include <unordered_map>
diff --git a/Sample/RT/MatrixRTCoefficients.cpp b/Sample/RT/MatrixRTCoefficients.cpp
index fe1d3181cc6d862a4de858503f9d320c740845db..aeef6194df18ad8d005b7117797d2abc812f112e 100644
--- a/Sample/RT/MatrixRTCoefficients.cpp
+++ b/Sample/RT/MatrixRTCoefficients.cpp
@@ -2,302 +2,163 @@
//
// BornAgain: simulate and fit scattering at grazing incidence
//
-//! @file Sample/RT/MatrixRTCoefficients.cpp
-//! @brief Implements class MatrixRTCoefficients.
+//! @file Sample/RT/MatrixRTCoefficients.cpp
+//! @brief Implements class MatrixRTCoefficients.
//!
//! @homepage http://www.bornagainproject.org
//! @license GNU General Public License v3 or higher (see COPYING)
-//! @copyright Forschungszentrum Jülich GmbH 2018
+//! @copyright Forschungszentrum Jülich GmbH 2020
//! @authors Scientific Computing Group at MLZ (see CITATION, AUTHORS)
//
// ************************************************************************************************
#include "Sample/RT/MatrixRTCoefficients.h"
+#include "Base/Utils/Assert.h"
+
+namespace {
+complex_t GetImExponential(complex_t exponent);
+const auto eps = std::numeric_limits<double>::epsilon() * 10.;
+} // namespace
+
+MatrixRTCoefficients::MatrixRTCoefficients(double kz_sign, Eigen::Vector2cd eigenvalues,
+ kvector_t b, double magnetic_SLD)
+ : m_kz_sign(kz_sign)
+ , m_lambda(std::move(eigenvalues))
+ , m_b(std::move(b))
+ , m_magnetic_SLD(magnetic_SLD) {
+ ASSERT(std::abs(m_b.mag() - 1) < eps || (m_b.mag() < eps && magnetic_SLD < eps));
+
+ m_T << 1, 0, 0, 1;
+ m_R << -1, 0, 0, -1;
+}
+
+MatrixRTCoefficients::MatrixRTCoefficients(const MatrixRTCoefficients& other) = default;
+
+MatrixRTCoefficients::~MatrixRTCoefficients() = default;
MatrixRTCoefficients* MatrixRTCoefficients::clone() const {
return new MatrixRTCoefficients(*this);
}
+Eigen::Matrix2cd MatrixRTCoefficients::TransformationMatrix(Eigen::Vector2d selection) const {
+ const Eigen::Matrix2cd exp2 = Eigen::DiagonalMatrix<complex_t, 2>(selection);
+
+ if (std::abs(m_b.mag() - 1.) < eps) {
+ Eigen::Matrix2cd Q;
+ const double factor1 = 2. * (1. + m_b.z());
+ Q << (1. + m_b.z()), (I * m_b.y() - m_b.x()), (m_b.x() + I * m_b.y()), (m_b.z() + 1.);
+ return Q * exp2 * Q.adjoint() / factor1;
+
+ } else if (m_b.mag() < eps)
+ return exp2;
+
+ throw std::runtime_error("Broken magnetic field vector");
+}
+
+Eigen::Matrix2cd MatrixRTCoefficients::T1Matrix() const {
+ return TransformationMatrix({0., 1.});
+}
+
+Eigen::Matrix2cd MatrixRTCoefficients::T2Matrix() const {
+ return TransformationMatrix({1., 0.});
+}
+
Eigen::Vector2cd MatrixRTCoefficients::T1plus() const {
- Eigen::Vector2cd result;
- Eigen::Matrix<complex_t, 4, 1> m = T1m * phi_psi_plus;
- result(0) = m(2);
- result(1) = m(3);
- if (lambda(0) == 0.0 && result == Eigen::Vector2cd::Zero())
- result(0) = 0.5;
- return result;
+ return T1Matrix() * m_T.col(0);
}
Eigen::Vector2cd MatrixRTCoefficients::R1plus() const {
- Eigen::Vector2cd result;
- Eigen::Matrix<complex_t, 4, 1> m = R1m * phi_psi_plus;
- result(0) = m(2);
- result(1) = m(3);
- Eigen::Matrix<complex_t, 4, 1> mT = T1m * phi_psi_plus;
- if (lambda(0) == 0.0 && mT(2) == 0.0 && mT(3) == 0.0)
- result(0) = -0.5;
- return result;
+ return T1Matrix() * m_R.col(0);
}
Eigen::Vector2cd MatrixRTCoefficients::T2plus() const {
- Eigen::Vector2cd result;
- Eigen::Matrix<complex_t, 4, 1> m = T2m * phi_psi_plus;
- result(0) = m(2);
- result(1) = m(3);
- if (lambda(1) == 0.0 && result == Eigen::Vector2cd::Zero())
- result(0) = 0.5;
- return result;
+ return T2Matrix() * m_T.col(0);
}
Eigen::Vector2cd MatrixRTCoefficients::R2plus() const {
- Eigen::Vector2cd result;
- Eigen::Matrix<complex_t, 4, 1> m = R2m * phi_psi_plus;
- result(0) = m(2);
- result(1) = m(3);
- Eigen::Matrix<complex_t, 4, 1> mT = T2m * phi_psi_plus;
- if (lambda(1) == 0.0 && mT(2) == 0.0 && mT(3) == 0.0)
- result(0) = -0.5;
- return result;
+ return T2Matrix() * m_R.col(0);
}
Eigen::Vector2cd MatrixRTCoefficients::T1min() const {
- Eigen::Vector2cd result;
- Eigen::Matrix<complex_t, 4, 1> m = T1m * phi_psi_min;
- result(0) = m(2);
- result(1) = m(3);
- if (lambda(0) == 0.0 && result == Eigen::Vector2cd::Zero())
- result(1) = 0.5;
- return result;
+ return T1Matrix() * m_T.col(1);
}
Eigen::Vector2cd MatrixRTCoefficients::R1min() const {
- Eigen::Vector2cd result;
- Eigen::Matrix<complex_t, 4, 1> m = R1m * phi_psi_min;
- result(0) = m(2);
- result(1) = m(3);
- Eigen::Matrix<complex_t, 4, 1> mT = T1m * phi_psi_min;
- if (lambda(0) == 0.0 && mT(2) == 0.0 && mT(3) == 0.0)
- result(1) = -0.5;
- return result;
+ return T1Matrix() * m_R.col(1);
}
Eigen::Vector2cd MatrixRTCoefficients::T2min() const {
- Eigen::Vector2cd result;
- Eigen::Matrix<complex_t, 4, 1> m = T2m * phi_psi_min;
- result(0) = m(2);
- result(1) = m(3);
- if (lambda(1) == 0.0 && result == Eigen::Vector2cd::Zero())
- result(1) = 0.5;
- return result;
+ return T2Matrix() * m_T.col(1);
}
Eigen::Vector2cd MatrixRTCoefficients::R2min() const {
- Eigen::Vector2cd result;
- Eigen::Matrix<complex_t, 4, 1> m = R2m * phi_psi_min;
- result(0) = m(2);
- result(1) = m(3);
- Eigen::Matrix<complex_t, 4, 1> mT = T2m * phi_psi_min;
- if (lambda(1) == 0.0 && mT(2) == 0.0 && mT(3) == 0.0)
- result(1) = -0.5;
+ return T2Matrix() * m_R.col(1);
+}
+
+Eigen::Vector2cd MatrixRTCoefficients::getKz() const {
+ return m_kz_sign * m_lambda;
+}
+
+Eigen::Matrix2cd MatrixRTCoefficients::pMatrixHelper(double sign) const {
+ const complex_t alpha = m_lambda(1) + m_lambda(0);
+ const complex_t beta = m_lambda(1) - m_lambda(0);
+
+ Eigen::Matrix2cd result;
+
+ kvector_t b = m_b;
+
+ result << alpha + sign * beta * b.z(), sign * beta * (b.x() - I * b.y()),
+ sign * beta * (b.x() + I * b.y()), alpha - sign * beta * b.z();
+
+ return m_kz_sign * result;
+}
+
+Eigen::Matrix2cd MatrixRTCoefficients::computeP() const {
+ Eigen::Matrix2cd result = pMatrixHelper(1.);
+ result *= 0.5;
+
return result;
}
-void MatrixRTCoefficients::calculateTRMatrices() {
- if (m_b_mag == 0.0) {
- calculateTRWithoutMagnetization();
- return;
- }
-
- if (lambda(0) == 0.0) {
- complex_t ikt = mul_I(m_kt);
- // Lambda1 component contained only in T1m (R1m=0)
- // row 0:
- T1m(0, 0) = (1.0 - m_bz / m_b_mag) / 2.0;
- T1m(0, 1) = -m_scatt_matrix(0, 1) / 2.0 / m_b_mag;
- // row 1:
- T1m(1, 0) = -m_scatt_matrix(1, 0) / 2.0 / m_b_mag;
- T1m(1, 1) = (1.0 + m_bz / m_b_mag) / 2.0;
- // row 2:
- T1m(2, 0) = ikt * (1.0 - m_bz / m_b_mag) / 2.0;
- T1m(2, 1) = -ikt * m_scatt_matrix(0, 1) / 2.0 / m_b_mag;
- T1m(2, 2) = T1m(0, 0);
- T1m(2, 3) = T1m(0, 1);
- // row 3:
- T1m(3, 0) = -ikt * m_scatt_matrix(1, 0) / 2.0 / m_b_mag;
- T1m(3, 1) = ikt * (1.0 + m_bz / m_b_mag) / 2.0;
- T1m(3, 2) = T1m(1, 0);
- T1m(3, 3) = T1m(1, 1);
- } else {
- // 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);
- }
-
- if (lambda(1) == 0.0) {
- complex_t ikt = mul_I(m_kt);
- // Lambda2 component contained only in T2m (R2m=0)
- // row 0:
- T2m(0, 0) = (1.0 + m_bz / m_b_mag) / 2.0;
- T2m(0, 1) = m_scatt_matrix(0, 1) / 2.0 / m_b_mag;
- // row 1:
- T2m(1, 0) = m_scatt_matrix(1, 0) / 2.0 / m_b_mag;
- T2m(1, 1) = (1.0 - m_bz / m_b_mag) / 2.0;
- // row 2:
- T2m(2, 0) = ikt * (1.0 + m_bz / m_b_mag) / 2.0;
- T2m(2, 1) = ikt * m_scatt_matrix(0, 1) / 2.0 / m_b_mag;
- T2m(2, 2) = T2m(0, 0);
- T2m(2, 3) = T2m(0, 1);
- // row 3:
- T2m(3, 0) = ikt * m_scatt_matrix(1, 0) / 2.0 / m_b_mag;
- T2m(3, 1) = ikt * (1.0 - m_bz / m_b_mag) / 2.0;
- T2m(3, 2) = T2m(1, 0);
- T2m(3, 3) = T2m(1, 1);
- } else {
- // 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);
- }
+Eigen::Matrix2cd MatrixRTCoefficients::computeInverseP() const {
+ const complex_t alpha = m_lambda(1) + m_lambda(0);
+ const complex_t beta = m_lambda(1) - m_lambda(0);
+
+ if (std::abs(alpha * alpha - beta * beta) == 0.)
+ return Eigen::Matrix2cd::Zero();
+
+ Eigen::Matrix2cd result = pMatrixHelper(-1.);
+ result *= 2. / (alpha * alpha - beta * beta);
+
+ return result;
}
-void MatrixRTCoefficients::calculateTRWithoutMagnetization() {
- T1m.setZero();
- R1m.setZero();
- T2m.setZero();
- R2m.setZero();
-
- if (m_a == 0.0) {
- // Spin down component contained only in T1 (R1=0)
- T1m(1, 1) = 1.0;
- T1m(3, 1) = mul_I(m_kt);
- T1m(3, 3) = 1.0;
-
- // Spin up component contained only in T2 (R2=0)
- T2m(0, 0) = 1.0;
- T2m(2, 0) = mul_I(m_kt);
- T2m(2, 2) = 1.0;
- return;
- }
-
- // T1m:
- 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(1, 1) = 0.5;
- R1m(1, 3) = std::sqrt(m_a) / 2.0;
- R1m(3, 1) = 1.0 / (2.0 * std::sqrt(m_a));
- R1m(3, 3) = 0.5;
-
- // T2m:
- T2m(0, 0) = 0.5;
- T2m(0, 2) = -std::sqrt(m_a) / 2.0;
- T2m(2, 0) = -1.0 / (2.0 * std::sqrt(m_a));
- T2m(2, 2) = 0.5;
-
- // R2m:
- 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;
+Eigen::Matrix2cd MatrixRTCoefficients::computeDeltaMatrix(double thickness) {
+ Eigen::Matrix2cd result;
+ const complex_t alpha = 0.5 * thickness * (m_lambda(1) + m_lambda(0));
+
+ const Eigen::Matrix2cd exp2 = Eigen::DiagonalMatrix<complex_t, 2>(
+ {GetImExponential(m_kz_sign * thickness * m_lambda(1)),
+ GetImExponential(m_kz_sign * thickness * m_lambda(0))});
+
+ // Compute resulting phase matrix according to exp(i p_m d_m) = exp1 * Q * exp2 * Q.adjoint();
+ if (std::abs(m_b.mag() - 1.) < eps) {
+ Eigen::Matrix2cd Q;
+ const double factor1 = 2. * (1. + m_b.z());
+ Q << (1. + m_b.z()), (I * m_b.y() - m_b.x()), (m_b.x() + I * m_b.y()), (m_b.z() + 1.);
+
+ return Q * exp2 * Q.adjoint() / factor1;
+
+ } else if (m_b.mag() < eps)
+ return Eigen::Matrix2cd::Identity() * GetImExponential(m_kz_sign * alpha);
+
+ throw std::runtime_error("Broken magnetic field vector");
}
-void MatrixRTCoefficients::initializeBottomLayerPhiPsi() {
- if (m_b_mag == 0.0) {
- phi_psi_min << 0.0, -std::sqrt(m_a), 0.0, 1.0;
- phi_psi_plus << -std::sqrt(m_a), 0.0, 1.0, 0.0;
- return;
- }
- // First basis vector that has no upward going wave amplitude
- phi_psi_min(0) = m_scatt_matrix(0, 1) * (lambda(0) - lambda(1)) / 2.0 / m_b_mag;
- phi_psi_min(1) = (m_bz * (lambda(1) - lambda(0)) / m_b_mag - lambda(1) - lambda(0)) / 2.0;
- phi_psi_min(2) = 0.0;
- phi_psi_min(3) = 1.0;
-
- // Second basis vector that has no upward going wave amplitude
- phi_psi_plus(0) = -(m_scatt_matrix(0, 0) + lambda(0) * lambda(1)) / (lambda(0) + lambda(1));
- phi_psi_plus(1) = m_scatt_matrix(1, 0) * (lambda(0) - lambda(1)) / 2.0 / m_b_mag;
- phi_psi_plus(2) = 1.0;
- phi_psi_plus(3) = 0.0;
+namespace {
+complex_t GetImExponential(complex_t exponent) {
+ if (exponent.imag() > -std::log(std::numeric_limits<double>::min()))
+ return 0.0;
+ return std::exp(I * exponent);
}
+} // namespace
diff --git a/Sample/RT/MatrixRTCoefficients.h b/Sample/RT/MatrixRTCoefficients.h
index 964a47d8a351c6be53befa459680d852dfe28262..0f467eda100970224f81c3a690803f87f488ee78 100644
--- a/Sample/RT/MatrixRTCoefficients.h
+++ b/Sample/RT/MatrixRTCoefficients.h
@@ -2,12 +2,12 @@
//
// BornAgain: simulate and fit scattering at grazing incidence
//
-//! @file Sample/RT/MatrixRTCoefficients.h
+//! @file Sample/RT/MatrixRTCoefficients.h
//! @brief Defines class MatrixRTCoefficients.
//!
//! @homepage http://www.bornagainproject.org
//! @license GNU General Public License v3 or higher (see COPYING)
-//! @copyright Forschungszentrum Jülich GmbH 2018
+//! @copyright Forschungszentrum Jülich GmbH 2020
//! @authors Scientific Computing Group at MLZ (see CITATION, AUTHORS)
//
// ************************************************************************************************
@@ -15,55 +15,63 @@
#ifndef BORNAGAIN_SAMPLE_RT_MATRIXRTCOEFFICIENTS_H
#define BORNAGAIN_SAMPLE_RT_MATRIXRTCOEFFICIENTS_H
+#include "Base/Vector/Vectors3D.h"
#include "Sample/RT/ILayerRTCoefficients.h"
+#include <vector>
//! Specular reflection and transmission coefficients in a layer in case
-//! of 2x2 matrix interactions between the layers and the scattered particle.
+//! of magnetic interactions between the scattered particle and the layer.
//! @ingroup algorithms_internal
class MatrixRTCoefficients : public ILayerRTCoefficients {
public:
- MatrixRTCoefficients() : m_kt(0.0) {}
- virtual ~MatrixRTCoefficients() {}
+ friend class SpecularMagneticStrategy;
+ friend class SpecularMagneticNCStrategy;
+ friend class SpecularMagneticTanhStrategy;
- virtual MatrixRTCoefficients* clone() const;
+ MatrixRTCoefficients(double kz_sign, Eigen::Vector2cd eigenvalues, kvector_t b,
+ double magnetic_SLD);
+ MatrixRTCoefficients(const MatrixRTCoefficients& other);
+ ~MatrixRTCoefficients() override;
+
+ MatrixRTCoefficients* clone() const override;
//! The following functions return the transmitted and reflected amplitudes
//! for different incoming beam polarizations and eigenmodes
- virtual Eigen::Vector2cd T1plus() const;
- virtual Eigen::Vector2cd R1plus() const;
- virtual Eigen::Vector2cd T2plus() const;
- virtual Eigen::Vector2cd R2plus() const;
- virtual Eigen::Vector2cd T1min() const;
- virtual Eigen::Vector2cd R1min() const;
- virtual Eigen::Vector2cd T2min() const;
- virtual Eigen::Vector2cd R2min() const;
+ Eigen::Vector2cd T1plus() const override;
+ Eigen::Vector2cd R1plus() const override;
+ Eigen::Vector2cd T2plus() const override;
+ Eigen::Vector2cd R2plus() const override;
+ Eigen::Vector2cd T1min() const override;
+ Eigen::Vector2cd R1min() const override;
+ Eigen::Vector2cd T2min() const override;
+ Eigen::Vector2cd R2min() const override;
//! Returns z-part of the two wavevector eigenmodes
- virtual Eigen::Vector2cd getKz() const { return kz; }
+ Eigen::Vector2cd getKz() const override;
+ double magneticSLD() const { return m_magnetic_SLD; }
+
+ Eigen::Matrix2cd computeP() const;
+ Eigen::Matrix2cd computeInverseP() const;
+
+ Eigen::Matrix2cd computeDeltaMatrix(double thickness);
+
+ Eigen::Matrix2cd getReflectionMatrix() const override { return m_R; };
+
+private:
+ double m_kz_sign; //! wave propagation direction (-1 for direct one, 1 for time reverse)
+ Eigen::Vector2cd m_lambda; //!< eigenvalues for wave propagation
+ kvector_t m_b; //!< unit magnetic field vector
+ double m_magnetic_SLD;
+
+ Eigen::Matrix2cd m_T;
+ Eigen::Matrix2cd m_R;
- void calculateTRMatrices();
- void calculateTRWithoutMagnetization();
- void initializeBottomLayerPhiPsi();
+ // helper functions to compute DWBA compatible amplitudes used in the T1plus() etc. functions
+ Eigen::Matrix2cd TransformationMatrix(Eigen::Vector2d selection) const;
+ Eigen::Matrix2cd T1Matrix() const;
+ Eigen::Matrix2cd T2Matrix() const;
- // NOTE: exceptionally, this class has member variables without prefix m_
- Eigen::Vector2cd kz; //!< z-part of the two wavevector eigenmodes
- Eigen::Vector2cd lambda; //!< positive eigenvalues of transfer matrix
- Eigen::Vector4cd phi_psi_plus; //!< boundary values for up-polarization
- Eigen::Vector4cd phi_psi_min; //!< boundary values for down-polarization
- Eigen::Matrix4cd T1m; //!< matrix selecting the transmitted part of
- //!< the first eigenmode
- Eigen::Matrix4cd R1m; //!< matrix selecting the reflected part of
- //!< the first eigenmode
- Eigen::Matrix4cd T2m; //!< matrix selecting the transmitted part of
- //!< the second eigenmode
- Eigen::Matrix4cd R2m; //!< matrix selecting the reflected part of
- //!< the second eigenmode
- Eigen::Matrix2cd m_scatt_matrix; //!< scattering matrix
- complex_t m_a; //!< polarization independent part of scattering matrix
- complex_t m_b_mag; //!< magnitude of magnetic interaction term
- complex_t m_bz; //!< z-part of magnetic interaction term
- double m_kt; //!< wavevector length times thickness of layer for use when
- //!< lambda=0
+ Eigen::Matrix2cd pMatrixHelper(double sign) const;
};
#endif // BORNAGAIN_SAMPLE_RT_MATRIXRTCOEFFICIENTS_H
diff --git a/Sample/RT/MatrixRTCoefficients_v3.cpp b/Sample/RT/MatrixRTCoefficients_v3.cpp
deleted file mode 100644
index 0e3ccbf5bbd659907c7acbeaf70d9a7851e76957..0000000000000000000000000000000000000000
--- a/Sample/RT/MatrixRTCoefficients_v3.cpp
+++ /dev/null
@@ -1,163 +0,0 @@
-// ************************************************************************************************
-//
-// BornAgain: simulate and fit scattering at grazing incidence
-//
-//! @file Sample/RT/MatrixRTCoefficients_v3.cpp
-//! @brief Implements class MatrixRTCoefficients_v3.
-//!
-//! @homepage http://www.bornagainproject.org
-//! @license GNU General Public License v3 or higher (see COPYING)
-//! @copyright Forschungszentrum Jülich GmbH 2020
-//! @authors Scientific Computing Group at MLZ (see CITATION, AUTHORS)
-//
-// ************************************************************************************************
-
-#include "Sample/RT/MatrixRTCoefficients_v3.h"
-#include "Base/Utils/Assert.h"
-
-namespace {
-complex_t GetImExponential(complex_t exponent);
-const auto eps = std::numeric_limits<double>::epsilon() * 10.;
-} // namespace
-
-MatrixRTCoefficients_v3::MatrixRTCoefficients_v3(double kz_sign, Eigen::Vector2cd eigenvalues,
- kvector_t b, double magnetic_SLD)
- : m_kz_sign(kz_sign)
- , m_lambda(std::move(eigenvalues))
- , m_b(std::move(b))
- , m_magnetic_SLD(magnetic_SLD) {
- ASSERT(std::abs(m_b.mag() - 1) < eps || (m_b.mag() < eps && magnetic_SLD < eps));
-
- m_T << 1, 0, 0, 1;
- m_R << -1, 0, 0, -1;
-}
-
-MatrixRTCoefficients_v3::MatrixRTCoefficients_v3(const MatrixRTCoefficients_v3& other) = default;
-
-MatrixRTCoefficients_v3::~MatrixRTCoefficients_v3() = default;
-
-MatrixRTCoefficients_v3* MatrixRTCoefficients_v3::clone() const {
- return new MatrixRTCoefficients_v3(*this);
-}
-
-Eigen::Matrix2cd MatrixRTCoefficients_v3::TransformationMatrix(Eigen::Vector2d selection) const {
- const Eigen::Matrix2cd exp2 = Eigen::DiagonalMatrix<complex_t, 2>(selection);
-
- if (std::abs(m_b.mag() - 1.) < eps) {
- Eigen::Matrix2cd Q;
- const double factor1 = 2. * (1. + m_b.z());
- Q << (1. + m_b.z()), (I * m_b.y() - m_b.x()), (m_b.x() + I * m_b.y()), (m_b.z() + 1.);
- return Q * exp2 * Q.adjoint() / factor1;
-
- } else if (m_b.mag() < eps)
- return exp2;
-
- throw std::runtime_error("Broken magnetic field vector");
-}
-
-Eigen::Matrix2cd MatrixRTCoefficients_v3::T1Matrix() const {
- return TransformationMatrix({0., 1.});
-}
-
-Eigen::Matrix2cd MatrixRTCoefficients_v3::T2Matrix() const {
- return TransformationMatrix({1., 0.});
-}
-
-Eigen::Vector2cd MatrixRTCoefficients_v3::T1plus() const {
- return T1Matrix() * m_T.col(0);
-}
-
-Eigen::Vector2cd MatrixRTCoefficients_v3::R1plus() const {
- return T1Matrix() * m_R.col(0);
-}
-
-Eigen::Vector2cd MatrixRTCoefficients_v3::T2plus() const {
- return T2Matrix() * m_T.col(0);
-}
-
-Eigen::Vector2cd MatrixRTCoefficients_v3::R2plus() const {
- return T2Matrix() * m_R.col(0);
-}
-
-Eigen::Vector2cd MatrixRTCoefficients_v3::T1min() const {
- return T1Matrix() * m_T.col(1);
-}
-
-Eigen::Vector2cd MatrixRTCoefficients_v3::R1min() const {
- return T1Matrix() * m_R.col(1);
-}
-
-Eigen::Vector2cd MatrixRTCoefficients_v3::T2min() const {
- return T2Matrix() * m_T.col(1);
-}
-
-Eigen::Vector2cd MatrixRTCoefficients_v3::R2min() const {
- return T2Matrix() * m_R.col(1);
-}
-
-Eigen::Vector2cd MatrixRTCoefficients_v3::getKz() const {
- return m_kz_sign * m_lambda;
-}
-
-Eigen::Matrix2cd MatrixRTCoefficients_v3::pMatrixHelper(double sign) const {
- const complex_t alpha = m_lambda(1) + m_lambda(0);
- const complex_t beta = m_lambda(1) - m_lambda(0);
-
- Eigen::Matrix2cd result;
-
- kvector_t b = m_b;
-
- result << alpha + sign * beta * b.z(), sign * beta * (b.x() - I * b.y()),
- sign * beta * (b.x() + I * b.y()), alpha - sign * beta * b.z();
-
- return result;
-}
-
-Eigen::Matrix2cd MatrixRTCoefficients_v3::computeP() const {
- Eigen::Matrix2cd result = pMatrixHelper(1.);
- result *= 0.5;
-
- return result;
-}
-
-Eigen::Matrix2cd MatrixRTCoefficients_v3::computeInverseP() const {
- const complex_t alpha = m_lambda(1) + m_lambda(0);
- const complex_t beta = m_lambda(1) - m_lambda(0);
-
- if (std::abs(alpha * alpha - beta * beta) == 0.)
- return Eigen::Matrix2cd::Zero();
-
- Eigen::Matrix2cd result = pMatrixHelper(-1.);
- result *= 2. / (alpha * alpha - beta * beta);
-
- return result;
-}
-
-Eigen::Matrix2cd MatrixRTCoefficients_v3::computeDeltaMatrix(double thickness) {
- Eigen::Matrix2cd result;
- const complex_t alpha = 0.5 * thickness * (m_lambda(1) + m_lambda(0));
-
- const Eigen::Matrix2cd exp2 = Eigen::DiagonalMatrix<complex_t, 2>(
- {GetImExponential(thickness * m_lambda(1)), GetImExponential(thickness * m_lambda(0))});
-
- // Compute resulting phase matrix according to exp(i p_m d_m) = exp1 * Q * exp2 * Q.adjoint();
- if (std::abs(m_b.mag() - 1.) < eps) {
- Eigen::Matrix2cd Q;
- const double factor1 = 2. * (1. + m_b.z());
- Q << (1. + m_b.z()), (I * m_b.y() - m_b.x()), (m_b.x() + I * m_b.y()), (m_b.z() + 1.);
-
- return Q * exp2 * Q.adjoint() / factor1;
-
- } else if (m_b.mag() < eps)
- return Eigen::Matrix2cd::Identity() * GetImExponential(alpha);
-
- throw std::runtime_error("Broken magnetic field vector");
-}
-
-namespace {
-complex_t GetImExponential(complex_t exponent) {
- if (exponent.imag() > -std::log(std::numeric_limits<double>::min()))
- return 0.0;
- return std::exp(I * exponent);
-}
-} // namespace
diff --git a/Sample/RT/MatrixRTCoefficients_v3.h b/Sample/RT/MatrixRTCoefficients_v3.h
deleted file mode 100644
index 7a239c1f011c2c589561cef9b6c51062f755ff39..0000000000000000000000000000000000000000
--- a/Sample/RT/MatrixRTCoefficients_v3.h
+++ /dev/null
@@ -1,82 +0,0 @@
-// ************************************************************************************************
-//
-// BornAgain: simulate and fit scattering at grazing incidence
-//
-//! @file Sample/RT/MatrixRTCoefficients_v3.h
-//! @brief Defines class MatrixRTCoefficients_v3.
-//!
-//! @homepage http://www.bornagainproject.org
-//! @license GNU General Public License v3 or higher (see COPYING)
-//! @copyright Forschungszentrum Jülich GmbH 2020
-//! @authors Scientific Computing Group at MLZ (see CITATION, AUTHORS)
-//
-// ************************************************************************************************
-
-#ifndef BORNAGAIN_SAMPLE_RT_MATRIXRTCOEFFICIENTS_V3_H
-#define BORNAGAIN_SAMPLE_RT_MATRIXRTCOEFFICIENTS_V3_H
-
-#include "Base/Vector/Vectors3D.h"
-#include "Sample/RT/ILayerRTCoefficients.h"
-#include <vector>
-
-//! Specular reflection and transmission coefficients in a layer in case
-//! of magnetic interactions between the scattered particle and the layer.
-//! @ingroup algorithms_internal
-
-class MatrixRTCoefficients_v3 : public ILayerRTCoefficients {
-public:
- friend class SpecularMagneticNewStrategy;
- friend class SpecularMagneticNewNCStrategy;
- friend class SpecularMagneticNewNCTestingStrategy;
- friend class SpecularMagneticNewTanhStrategy;
- friend class SpecularMagnetic_v3ConsistencyTest;
- friend class SpecularMagnetic_v3ConsistencyTest_ScalarMagneticAmplitudes_Test;
- friend class SpecularMagnetic_v3ConsistencyTest_AmplitudesBackwardsBackwards_Test;
- template <class sampleClass> friend class TestSimulation;
-
- MatrixRTCoefficients_v3(double kz_sign, Eigen::Vector2cd eigenvalues, kvector_t b,
- double magnetic_SLD);
- MatrixRTCoefficients_v3(const MatrixRTCoefficients_v3& other);
- ~MatrixRTCoefficients_v3() override;
-
- MatrixRTCoefficients_v3* clone() const override;
-
- //! The following functions return the transmitted and reflected amplitudes
- //! for different incoming beam polarizations and eigenmodes
- Eigen::Vector2cd T1plus() const override;
- Eigen::Vector2cd R1plus() const override;
- Eigen::Vector2cd T2plus() const override;
- Eigen::Vector2cd R2plus() const override;
- Eigen::Vector2cd T1min() const override;
- Eigen::Vector2cd R1min() const override;
- Eigen::Vector2cd T2min() const override;
- Eigen::Vector2cd R2min() const override;
- //! Returns z-part of the two wavevector eigenmodes
- Eigen::Vector2cd getKz() const override;
- double magneticSLD() const { return m_magnetic_SLD; }
-
- Eigen::Matrix2cd computeP() const;
- Eigen::Matrix2cd computeInverseP() const;
-
- Eigen::Matrix2cd computeDeltaMatrix(double thickness);
-
- Eigen::Matrix2cd getReflectionMatrix() const override { return m_R; };
-
-private:
- double m_kz_sign; //! wave propagation direction (-1 for direct one, 1 for time reverse)
- Eigen::Vector2cd m_lambda; //!< eigenvalues for wave propagation
- kvector_t m_b; //!< unit magnetic field vector
- double m_magnetic_SLD;
-
- Eigen::Matrix2cd m_T;
- Eigen::Matrix2cd m_R;
-
- // helper functions to compute DWBA compatible amplitudes used in the T1plus() etc. functions
- Eigen::Matrix2cd TransformationMatrix(Eigen::Vector2d selection) const;
- Eigen::Matrix2cd T1Matrix() const;
- Eigen::Matrix2cd T2Matrix() const;
-
- Eigen::Matrix2cd pMatrixHelper(double sign) const;
-};
-
-#endif // BORNAGAIN_SAMPLE_RT_MATRIXRTCOEFFICIENTS_V3_H
diff --git a/Sample/Specular/SpecularMagneticNewNCStrategy.cpp b/Sample/Specular/SpecularMagneticNCStrategy.cpp
similarity index 89%
rename from Sample/Specular/SpecularMagneticNewNCStrategy.cpp
rename to Sample/Specular/SpecularMagneticNCStrategy.cpp
index 5fd3633ac4ff66eab9e94a8769f82da7530561c2..85746b1327d6c4af8600f0ecabf4d01f06e68cbc 100644
--- a/Sample/Specular/SpecularMagneticNewNCStrategy.cpp
+++ b/Sample/Specular/SpecularMagneticNCStrategy.cpp
@@ -2,8 +2,8 @@
//
// BornAgain: simulate and fit scattering at grazing incidence
//
-//! @file Sample/Specular/SpecularMagneticNewNCStrategy.cpp
-//! @brief Implements class SpecularMagneticNewNCStrategy.
+//! @file Sample/Specular/SpecularMagneticNCStrategy.cpp
+//! @brief Implements class SpecularMagneticNCStrategy.
//!
//! @homepage http://www.bornagainproject.org
//! @license GNU General Public License v3 or higher (see COPYING)
@@ -12,15 +12,15 @@
//
// ************************************************************************************************
-#include "Sample/Specular/SpecularMagneticNewNCStrategy.h"
+#include "Sample/Specular/SpecularMagneticNCStrategy.h"
namespace {
complex_t checkForUnderflow(complex_t val);
}
std::pair<Eigen::Matrix2cd, Eigen::Matrix2cd>
-SpecularMagneticNewNCStrategy::computeRoughnessMatrices(const MatrixRTCoefficients_v3& coeff_i,
- const MatrixRTCoefficients_v3& coeff_i1,
+SpecularMagneticNCStrategy::computeRoughnessMatrices(const MatrixRTCoefficients& coeff_i,
+ const MatrixRTCoefficients& coeff_i1,
double sigma) const {
complex_t beta_i = coeff_i.m_lambda(1) - coeff_i.m_lambda(0);
complex_t beta_i1 = coeff_i1.m_lambda(1) - coeff_i1.m_lambda(0);
@@ -65,8 +65,8 @@ SpecularMagneticNewNCStrategy::computeRoughnessMatrices(const MatrixRTCoefficien
}
std::pair<Eigen::Matrix2cd, Eigen::Matrix2cd>
-SpecularMagneticNewNCStrategy::computeBackwardsSubmatrices(const MatrixRTCoefficients_v3& coeff_i,
- const MatrixRTCoefficients_v3& coeff_i1,
+SpecularMagneticNCStrategy::computeBackwardsSubmatrices(const MatrixRTCoefficients& coeff_i,
+ const MatrixRTCoefficients& coeff_i1,
double sigma) const {
Eigen::Matrix2cd roughness_sum{Eigen::Matrix2cd::Identity()};
Eigen::Matrix2cd roughness_diff{Eigen::Matrix2cd::Identity()};
diff --git a/Sample/Specular/SpecularMagneticNewNCStrategy.h b/Sample/Specular/SpecularMagneticNCStrategy.h
similarity index 59%
rename from Sample/Specular/SpecularMagneticNewNCStrategy.h
rename to Sample/Specular/SpecularMagneticNCStrategy.h
index 66d402287521912e9bd4f42ce3f3a346829ed5e4..573c5b9e7a322dd09705f8694beb526aa29b49af 100644
--- a/Sample/Specular/SpecularMagneticNewNCStrategy.h
+++ b/Sample/Specular/SpecularMagneticNCStrategy.h
@@ -2,8 +2,8 @@
//
// BornAgain: simulate and fit scattering at grazing incidence
//
-//! @file Sample/Specular/SpecularMagneticNewNCStrategy.h
-//! @brief Defines class SpecularMagneticNewNCStrategy.
+//! @file Sample/Specular/SpecularMagneticNCStrategy.h
+//! @brief Defines class SpecularMagneticNCStrategy.
//!
//! @homepage http://www.bornagainproject.org
//! @license GNU General Public License v3 or higher (see COPYING)
@@ -12,10 +12,10 @@
//
// ************************************************************************************************
-#ifndef BORNAGAIN_SAMPLE_SPECULAR_SPECULARMAGNETICNEWNCSTRATEGY_H
-#define BORNAGAIN_SAMPLE_SPECULAR_SPECULARMAGNETICNEWNCSTRATEGY_H
+#ifndef BORNAGAIN_SAMPLE_SPECULAR_SPECULARMAGNETICNCSTRATEGY_H
+#define BORNAGAIN_SAMPLE_SPECULAR_SPECULARMAGNETICNCSTRATEGY_H
-#include "Sample/Specular/SpecularMagneticNewStrategy.h"
+#include "Sample/Specular/SpecularMagneticStrategy.h"
#include <memory>
#include <vector>
@@ -27,15 +27,15 @@
//! document "Polarized Implementation of the Transfer Matrix Method"
//!
//! @ingroup algorithms_internal
-class SpecularMagneticNewNCStrategy : public SpecularMagneticNewStrategy {
+class SpecularMagneticNCStrategy : public SpecularMagneticStrategy {
private:
std::pair<Eigen::Matrix2cd, Eigen::Matrix2cd>
- computeRoughnessMatrices(const MatrixRTCoefficients_v3& coeff_i,
- const MatrixRTCoefficients_v3& coeff_i1, double sigma) const;
+ computeRoughnessMatrices(const MatrixRTCoefficients& coeff_i,
+ const MatrixRTCoefficients& coeff_i1, double sigma) const;
virtual std::pair<Eigen::Matrix2cd, Eigen::Matrix2cd>
- computeBackwardsSubmatrices(const MatrixRTCoefficients_v3& coeff_i,
- const MatrixRTCoefficients_v3& coeff_i1, double sigma) const;
+ computeBackwardsSubmatrices(const MatrixRTCoefficients& coeff_i,
+ const MatrixRTCoefficients& coeff_i1, double sigma) const;
};
-#endif // BORNAGAIN_SAMPLE_SPECULAR_SPECULARMAGNETICNEWNCSTRATEGY_H
+#endif // BORNAGAIN_SAMPLE_SPECULAR_SPECULARMAGNETICNCSTRATEGY_H
diff --git a/Sample/Specular/SpecularMagneticNewStrategy.cpp b/Sample/Specular/SpecularMagneticNewStrategy.cpp
deleted file mode 100644
index cb71ab7c766222053a34063244e123d5c65eadbd..0000000000000000000000000000000000000000
--- a/Sample/Specular/SpecularMagneticNewStrategy.cpp
+++ /dev/null
@@ -1,183 +0,0 @@
-// ************************************************************************************************
-//
-// BornAgain: simulate and fit scattering at grazing incidence
-//
-//! @file Sample/Specular/SpecularMagneticNewStrategy.cpp
-//! @brief Implements class SpecularMagneticNewStrategy.
-//!
-//! @homepage http://www.bornagainproject.org
-//! @license GNU General Public License v3 or higher (see COPYING)
-//! @copyright Forschungszentrum Jülich GmbH 2020
-//! @authors Scientific Computing Group at MLZ (see CITATION, AUTHORS)
-//
-// ************************************************************************************************
-
-#include "Sample/Specular/SpecularMagneticNewStrategy.h"
-#include "Base/Const/PhysicalConstants.h"
-#include "Sample/Slice/KzComputation.h"
-#include "Sample/Slice/LayerRoughness.h"
-#include "Sample/Slice/Slice.h"
-
-namespace {
-double magneticSLD(kvector_t B_field);
-Eigen::Vector2cd eigenvalues(complex_t kz, double b_mag);
-Eigen::Vector2cd checkForUnderflow(const Eigen::Vector2cd& eigenvs);
-
-// The factor 1e-18 is here to have unit: 1/T*nm^-2
-constexpr double magnetic_prefactor = PhysConsts::m_n * PhysConsts::g_factor_n * PhysConsts::mu_N
- / PhysConsts::h_bar / PhysConsts::h_bar / 4. / M_PI * 1e-18;
-const auto eps = std::numeric_limits<double>::epsilon() * 10.;
-const LayerRoughness* GetBottomRoughness(const std::vector<Slice>& slices,
- const size_t slice_index);
-} // namespace
-
-ISpecularStrategy::coeffs_t SpecularMagneticNewStrategy::Execute(const std::vector<Slice>& slices,
- const kvector_t& k) const {
- return Execute(slices, KzComputation::computeReducedKz(slices, k));
-}
-
-ISpecularStrategy::coeffs_t
-SpecularMagneticNewStrategy::Execute(const std::vector<Slice>& slices,
- const std::vector<complex_t>& kz) const {
- if (slices.size() != kz.size())
- throw std::runtime_error("Number of slices does not match the size of the kz-vector");
-
- ISpecularStrategy::coeffs_t result;
- for (auto& coeff : computeTR(slices, kz))
- result.push_back(std::make_unique<MatrixRTCoefficients_v3>(coeff));
-
- return result;
-}
-
-std::vector<MatrixRTCoefficients_v3>
-SpecularMagneticNewStrategy::computeTR(const std::vector<Slice>& slices,
- const std::vector<complex_t>& kzs) const {
- const size_t N = slices.size();
-
- if (slices.size() != kzs.size())
- throw std::runtime_error(
- "Error in SpecularMagnetic_::execute: kz vector and slices size shall coinside.");
- if (slices.empty())
- return {};
-
- std::vector<MatrixRTCoefficients_v3> result;
- result.reserve(N);
-
- const double kz_sign = kzs.front().real() >= 0.0 ? 1.0 : -1.0; // save sign to restore it later
-
- auto B_0 = slices.front().bField();
- result.emplace_back(kz_sign, eigenvalues(kzs.front(), 0.0), kvector_t{0.0, 0.0, 0.0}, 0.0);
- for (size_t i = 1, size = slices.size(); i < size; ++i) {
- auto B = slices[i].bField() - B_0;
- auto magnetic_SLD = magneticSLD(B);
- result.emplace_back(kz_sign, checkForUnderflow(eigenvalues(kzs[i], magnetic_SLD)),
- B.mag() > eps ? B / B.mag() : kvector_t{0.0, 0.0, 0.0}, magnetic_SLD);
- }
-
- if (N == 1) {
- result[0].m_T = Eigen::Matrix2cd::Identity();
- result[0].m_R = Eigen::Matrix2cd::Zero();
- return result;
-
- } else if (kzs[0] == 0.) {
- result[0].m_T = Eigen::Matrix2cd::Identity();
- result[0].m_R = -Eigen::Matrix2cd::Identity();
- for (size_t i = 1; i < N; ++i) {
- result[i].m_T.setZero();
- result[i].m_R.setZero();
- }
- return result;
- }
-
- calculateUpwards(result, slices);
-
- return result;
-}
-
-void SpecularMagneticNewStrategy::calculateUpwards(std::vector<MatrixRTCoefficients_v3>& coeff,
- const std::vector<Slice>& slices) const {
- const auto N = slices.size();
- std::vector<Eigen::Matrix2cd> SMatrices(N - 1);
- std::vector<complex_t> Normalization(N - 1);
-
- // bottom boundary condition
- coeff.back().m_T = Eigen::Matrix2cd::Identity();
- coeff.back().m_R = Eigen::Matrix2cd::Zero();
-
- for (int i = N - 2; i >= 0; --i) {
- double sigma = 0.;
- if (const auto roughness = GetBottomRoughness(slices, i))
- sigma = roughness->getSigma();
-
- // compute the 2x2 submatrices in the write-up denoted as P+, P- and delta
- const auto mpmm = computeBackwardsSubmatrices(coeff[i], coeff[i + 1], sigma);
- const Eigen::Matrix2cd mp = mpmm.first;
- const Eigen::Matrix2cd mm = mpmm.second;
-
- const Eigen::Matrix2cd delta = coeff[i].computeDeltaMatrix(slices[i].thickness());
-
- // compute the rotation matrix
- Eigen::Matrix2cd S, Si;
- Si = mp + mm * coeff[i + 1].m_R;
- S << Si(1, 1), -Si(0, 1), -Si(1, 0), Si(0, 0);
- const complex_t norm = S(0, 0) * S(1, 1) - S(0, 1) * S(1, 0);
- S = S * delta;
-
- // store the rotation matrix and normalization constant in order to rotate
- // the coefficients for all lower slices at the end of the computation
- SMatrices[i] = S;
- Normalization[i] = norm;
-
- // compute the reflection matrix and
- // rotate the polarization such that we have pure incoming states (T = I)
- S /= norm;
-
- // T is always equal to the identity at this point, no need to store
- coeff[i].m_R = delta * (mm + mp * coeff[i + 1].m_R) * S;
- }
-
- // now correct all amplitudes in forward direction by dividing with the remaining
- // normalization constants. In addition rotate the polarization by the amount
- // that was rotated above the current interface
- // if the normalization overflows, all amplitudes below that point are set to zero
- complex_t dumpingFactor = 1;
- Eigen::Matrix2cd S = Eigen::Matrix2cd::Identity();
- for (size_t i = 1; i < N; ++i) {
- dumpingFactor = dumpingFactor * Normalization[i - 1];
- S = SMatrices[i - 1] * S;
-
- if (std::isinf(std::norm(dumpingFactor))) {
- std::for_each(coeff.begin() + i, coeff.end(), [](auto& coeff) {
- coeff.m_T = Eigen::Matrix2cd::Zero();
- coeff.m_R = Eigen::Matrix2cd::Zero();
- });
- break;
- }
-
- coeff[i].m_T = S / dumpingFactor; // T * S omitted, since T is always I
- coeff[i].m_R *= S / dumpingFactor;
- }
-}
-
-namespace {
-double magneticSLD(kvector_t B_field) {
- return magnetic_prefactor * B_field.mag();
-}
-
-Eigen::Vector2cd eigenvalues(complex_t kz, double magnetic_SLD) {
- const complex_t a = kz * kz;
- return {std::sqrt(a - 4. * M_PI * magnetic_SLD), std::sqrt(a + 4. * M_PI * magnetic_SLD)};
-}
-
-Eigen::Vector2cd checkForUnderflow(const Eigen::Vector2cd& eigenvs) {
- auto lambda = [](complex_t value) { return std::abs(value) < 1e-40 ? 1e-40 : value; };
- return {lambda(eigenvs(0)), lambda(eigenvs(1))};
-}
-
-const LayerRoughness* GetBottomRoughness(const std::vector<Slice>& slices,
- const size_t slice_index) {
- if (slice_index + 1 < slices.size())
- return slices[slice_index + 1].topRoughness();
- return nullptr;
-}
-} // namespace
diff --git a/Sample/Specular/SpecularMagneticNewStrategy.h b/Sample/Specular/SpecularMagneticNewStrategy.h
deleted file mode 100644
index fc7e85f3afcee0b1061bfb2cbe43acff2f5991f1..0000000000000000000000000000000000000000
--- a/Sample/Specular/SpecularMagneticNewStrategy.h
+++ /dev/null
@@ -1,56 +0,0 @@
-// ************************************************************************************************
-//
-// BornAgain: simulate and fit scattering at grazing incidence
-//
-//! @file Sample/Specular/SpecularMagneticNewStrategy.h
-//! @brief Defines class SpecularMagneticNewStrategy.
-//!
-//! @homepage http://www.bornagainproject.org
-//! @license GNU General Public License v3 or higher (see COPYING)
-//! @copyright Forschungszentrum Jülich GmbH 2020
-//! @authors Scientific Computing Group at MLZ (see CITATION, AUTHORS)
-//
-// ************************************************************************************************
-
-#ifndef BORNAGAIN_SAMPLE_SPECULAR_SPECULARMAGNETICNEWSTRATEGY_H
-#define BORNAGAIN_SAMPLE_SPECULAR_SPECULARMAGNETICNEWSTRATEGY_H
-
-#include "Sample/RT/MatrixRTCoefficients_v3.h"
-#include "Sample/Specular/ISpecularStrategy.h"
-#include <memory>
-#include <vector>
-
-class Slice;
-
-//! Implements the magnetic Fresnel computation with Nevot-Croce roughness
-//!
-//! Implements the transfer matrix formalism for the calculation of wave
-//! amplitudes of the coherent wave solution in a multilayer with magnetization.
-//! For a description, see internal
-//! document "Polarized Implementation of the Transfer Matrix Method"
-//!
-//! @ingroup algorithms_internal
-class SpecularMagneticNewStrategy : public ISpecularStrategy {
-public:
- //! Computes refraction angle reflection/transmission coefficients
- //! for given sliced multilayer and wavevector k
- ISpecularStrategy::coeffs_t Execute(const std::vector<Slice>& slices, const kvector_t& k) const;
-
- //! Computes refraction angle reflection/transmission coefficients
- //! for given sliced multilayer and a set of kz projections corresponding to each slice
- ISpecularStrategy::coeffs_t Execute(const std::vector<Slice>& slices,
- const std::vector<complex_t>& kz) const;
-
-private:
- std::vector<MatrixRTCoefficients_v3> computeTR(const std::vector<Slice>& slices,
- const std::vector<complex_t>& kzs) const;
-
- virtual std::pair<Eigen::Matrix2cd, Eigen::Matrix2cd>
- computeBackwardsSubmatrices(const MatrixRTCoefficients_v3& coeff_i,
- const MatrixRTCoefficients_v3& coeff_i1, double sigma) const = 0;
-
- void calculateUpwards(std::vector<MatrixRTCoefficients_v3>& coeff,
- const std::vector<Slice>& slices) const;
-};
-
-#endif // BORNAGAIN_SAMPLE_SPECULAR_SPECULARMAGNETICNEWSTRATEGY_H
diff --git a/Sample/Specular/SpecularMagneticOldStrategy.h b/Sample/Specular/SpecularMagneticOldStrategy.h
deleted file mode 100644
index 88e6b17d87a4fe36201f586c7049bb7fbd595f39..0000000000000000000000000000000000000000
--- a/Sample/Specular/SpecularMagneticOldStrategy.h
+++ /dev/null
@@ -1,39 +0,0 @@
-// ************************************************************************************************
-//
-// BornAgain: simulate and fit scattering at grazing incidence
-//
-//! @file Sample/Specular/SpecularMagneticOldStrategy.h
-//! @brief Defines class SpecularMagneticOldStrategy.
-//!
-//! @homepage http://www.bornagainproject.org
-//! @license GNU General Public License v3 or higher (see COPYING)
-//! @copyright Forschungszentrum Jülich GmbH 2018
-//! @authors Scientific Computing Group at MLZ (see CITATION, AUTHORS)
-//
-// ************************************************************************************************
-
-#ifndef BORNAGAIN_SAMPLE_SPECULAR_SPECULARMAGNETICOLDSTRATEGY_H
-#define BORNAGAIN_SAMPLE_SPECULAR_SPECULARMAGNETICOLDSTRATEGY_H
-
-#include "Sample/RT/MatrixRTCoefficients.h"
-#include "Sample/Specular/ISpecularStrategy.h"
-#include <memory>
-#include <vector>
-
-class Slice;
-
-//! Implements the matrix formalism for the calculation of wave amplitudes of
-//! the coherent wave solution in a multilayer with magnetization.
-//! @ingroup algorithms_internal
-
-class SpecularMagneticOldStrategy : public ISpecularStrategy {
-public:
- //! Computes refraction angle reflection/transmission coefficients
- //! for given sliced multilayer and wavevector k
- coeffs_t Execute(const std::vector<Slice>& slices, const kvector_t& k) const;
-
- coeffs_t Execute(const std::vector<Slice>& slices, const std::vector<complex_t>& kz) const;
-
-}; // class SpecularMagneticOldStrategy
-
-#endif // BORNAGAIN_SAMPLE_SPECULAR_SPECULARMAGNETICOLDSTRATEGY_H
diff --git a/Sample/Specular/SpecularMagneticStrategy.cpp b/Sample/Specular/SpecularMagneticStrategy.cpp
index b46fca7567fd7a7700ecc2e9945649524eb651f3..e940ff0093f780eea1efa006bd760c139848e5ce 100644
--- a/Sample/Specular/SpecularMagneticStrategy.cpp
+++ b/Sample/Specular/SpecularMagneticStrategy.cpp
@@ -7,7 +7,7 @@
//!
//! @homepage http://www.bornagainproject.org
//! @license GNU General Public License v3 or higher (see COPYING)
-//! @copyright Forschungszentrum Jülich GmbH 2018
+//! @copyright Forschungszentrum Jülich GmbH 2020
//! @authors Scientific Computing Group at MLZ (see CITATION, AUTHORS)
//
// ************************************************************************************************
@@ -15,42 +15,44 @@
#include "Sample/Specular/SpecularMagneticStrategy.h"
#include "Base/Const/PhysicalConstants.h"
#include "Sample/Slice/KzComputation.h"
+#include "Sample/Slice/LayerRoughness.h"
#include "Sample/Slice/Slice.h"
namespace {
-kvector_t magneticImpact(kvector_t B_field);
+double magneticSLD(kvector_t B_field);
Eigen::Vector2cd eigenvalues(complex_t kz, double b_mag);
Eigen::Vector2cd checkForUnderflow(const Eigen::Vector2cd& eigenvs);
-complex_t GetImExponential(complex_t exponent);
// The factor 1e-18 is here to have unit: 1/T*nm^-2
constexpr double magnetic_prefactor = PhysConsts::m_n * PhysConsts::g_factor_n * PhysConsts::mu_N
- / PhysConsts::h_bar / PhysConsts::h_bar * 1e-18;
+ / PhysConsts::h_bar / PhysConsts::h_bar / 4. / M_PI * 1e-18;
+const auto eps = std::numeric_limits<double>::epsilon() * 10.;
+const LayerRoughness* GetBottomRoughness(const std::vector<Slice>& slices,
+ const size_t slice_index);
} // namespace
ISpecularStrategy::coeffs_t SpecularMagneticStrategy::Execute(const std::vector<Slice>& slices,
- const kvector_t& k) const {
+ const kvector_t& k) const {
return Execute(slices, KzComputation::computeReducedKz(slices, k));
}
ISpecularStrategy::coeffs_t
SpecularMagneticStrategy::Execute(const std::vector<Slice>& slices,
- const std::vector<complex_t>& kz) const {
+ const std::vector<complex_t>& kz) const {
if (slices.size() != kz.size())
throw std::runtime_error("Number of slices does not match the size of the kz-vector");
ISpecularStrategy::coeffs_t result;
for (auto& coeff : computeTR(slices, kz))
- result.push_back(std::make_unique<MatrixRTCoefficients_v2>(coeff));
+ result.push_back(std::make_unique<MatrixRTCoefficients>(coeff));
return result;
}
-std::vector<MatrixRTCoefficients_v2>
+std::vector<MatrixRTCoefficients>
SpecularMagneticStrategy::computeTR(const std::vector<Slice>& slices,
- const std::vector<complex_t>& kzs) {
- if (kzs[0] == 0.)
- throw std::runtime_error("Edge case k_z = 0 not implemented");
+ const std::vector<complex_t>& kzs) const {
+ const size_t N = slices.size();
if (slices.size() != kzs.size())
throw std::runtime_error(
@@ -58,206 +60,113 @@ SpecularMagneticStrategy::computeTR(const std::vector<Slice>& slices,
if (slices.empty())
return {};
- std::vector<MatrixRTCoefficients_v2> result;
- result.reserve(slices.size());
+ std::vector<MatrixRTCoefficients> result;
+ result.reserve(N);
- const double kz_sign = kzs.front().real() > 0.0 ? 1.0 : -1.0; // save sign to restore it later
- const kvector_t b_0 = magneticImpact(slices.front().bField());
- result.emplace_back(kz_sign, eigenvalues(kzs.front(), 0.0), kvector_t{0.0, 0.0, 0.0});
+ const double kz_sign = kzs.front().real() >= 0.0 ? 1.0 : -1.0; // save sign to restore it later
+
+ auto B_0 = slices.front().bField();
+ result.emplace_back(kz_sign, eigenvalues(kzs.front(), 0.0), kvector_t{0.0, 0.0, 0.0}, 0.0);
for (size_t i = 1, size = slices.size(); i < size; ++i) {
- kvector_t b = magneticImpact(slices[i].bField()) - b_0;
- result.emplace_back(kz_sign, checkForUnderflow(eigenvalues(kzs[i], b.mag())), b);
+ auto B = slices[i].bField() - B_0;
+ auto magnetic_SLD = magneticSLD(B);
+ result.emplace_back(kz_sign, checkForUnderflow(eigenvalues(kzs[i], magnetic_SLD)),
+ B.mag() > eps ? B / B.mag() : kvector_t{0.0, 0.0, 0.0}, magnetic_SLD);
}
- if (result.front().m_lambda == Eigen::Vector2cd::Zero()) {
- std::for_each(result.begin(), result.end(), [](auto& coeff) { setNoTransmission(coeff); });
+ if (N == 1) {
+ result[0].m_T = Eigen::Matrix2cd::Identity();
+ result[0].m_R = Eigen::Matrix2cd::Zero();
return result;
- }
-
- std::for_each(result.begin(), result.end(), [](auto& coeff) { calculateTR(coeff); });
- nullifyBottomReflection(result.back());
- propagateBackwardsForwards(result, slices);
-
- return result;
-}
-void SpecularMagneticStrategy::calculateTR(MatrixRTCoefficients_v2& coeff) {
- const double b = coeff.m_b.mag();
- if (b == 0.0) {
- calculateZeroFieldTR(coeff);
- return;
+ } else if (kzs[0] == 0.) {
+ result[0].m_T = Eigen::Matrix2cd::Identity();
+ result[0].m_R = -Eigen::Matrix2cd::Identity();
+ for (size_t i = 1; i < N; ++i) {
+ result[i].m_T.setZero();
+ result[i].m_R.setZero();
+ }
+ return result;
}
- const double bpbz = b + coeff.m_b.z();
- const double bmbz = b - coeff.m_b.z();
- const complex_t bxmby = coeff.m_b.x() - I * coeff.m_b.y();
- const complex_t bxpby = coeff.m_b.x() + I * coeff.m_b.y();
- const complex_t l_1 = coeff.m_lambda(0);
- const complex_t l_2 = coeff.m_lambda(1);
-
- auto& T1 = coeff.T1;
- T1 << bmbz, -bxmby, -bmbz * l_1, bxmby * l_1, -bxpby, bpbz, bxpby * l_1, -bpbz * l_1,
- -bmbz / l_1, bxmby / l_1, bmbz, -bxmby, bxpby / l_1, -bpbz / l_1, -bxpby, bpbz;
- T1 /= 4.0 * b;
-
- auto& R1 = coeff.R1;
- R1 << T1(0, 0), T1(0, 1), -T1(0, 2), -T1(0, 3), T1(1, 0), T1(1, 1), -T1(1, 2), -T1(1, 3),
- -T1(2, 0), -T1(2, 1), T1(2, 2), T1(2, 3), -T1(3, 0), -T1(3, 1), T1(3, 2), T1(3, 3);
-
- auto& T2 = coeff.T2;
- T2 << bpbz, bxmby, -bpbz * l_2, -bxmby * l_2, bxpby, bmbz, -bxpby * l_2, -bmbz * l_2,
- -bpbz / l_2, -bxmby / l_2, bpbz, bxmby, -bxpby / l_2, -bmbz / l_2, bxpby, bmbz;
- T2 /= 4.0 * b;
-
- auto& R2 = coeff.R2;
- R2 << T2(0, 0), T2(0, 1), -T2(0, 2), -T2(0, 3), T2(1, 0), T2(1, 1), -T2(1, 2), -T2(1, 3),
- -T2(2, 0), -T2(2, 1), T2(2, 2), T2(2, 3), -T2(3, 0), -T2(3, 1), T2(3, 2), T2(3, 3);
-}
-
-void SpecularMagneticStrategy::calculateZeroFieldTR(MatrixRTCoefficients_v2& coeff) {
- coeff.T1 = Eigen::Matrix4cd::Zero();
- coeff.R1 = Eigen::Matrix4cd::Zero();
- coeff.T2 = Eigen::Matrix4cd::Zero();
- coeff.R2 = Eigen::Matrix4cd::Zero();
-
- // lambda_1 == lambda_2, no difference which one to use
- const complex_t eigen_value = coeff.m_lambda(0);
-
- Eigen::Matrix3cd Tblock;
- Tblock << 0.5, 0.0, -0.5 * eigen_value, 0.0, 0.0, 0.0, -0.5 / eigen_value, 0.0, 0.5;
-
- Eigen::Matrix3cd Rblock;
- Rblock << 0.5, 0.0, 0.5 * eigen_value, 0.0, 0.0, 0.0, 0.5 / eigen_value, 0.0, 0.5;
-
- coeff.T1.block<3, 3>(1, 1) = Tblock;
- coeff.R1.block<3, 3>(1, 1) = Rblock;
- coeff.T2.block<3, 3>(0, 0) = Tblock;
- coeff.R2.block<3, 3>(0, 0) = Rblock;
-}
+ calculateUpwards(result, slices);
-void SpecularMagneticStrategy::setNoTransmission(MatrixRTCoefficients_v2& coeff) {
- coeff.m_w_plus = Eigen::Vector4cd::Zero();
- coeff.m_w_min = Eigen::Vector4cd::Zero();
- coeff.T1 = Eigen::Matrix4cd::Identity() / 4.0;
- coeff.R1 = coeff.T1;
- coeff.T2 = coeff.T1;
- coeff.R2 = coeff.T1;
-}
-
-void SpecularMagneticStrategy::nullifyBottomReflection(MatrixRTCoefficients_v2& coeff) {
- const complex_t l_1 = coeff.m_lambda(0);
- const complex_t l_2 = coeff.m_lambda(1);
- const double b_mag = coeff.m_b.mag();
- const kvector_t& b = coeff.m_b;
-
- if (b_mag == 0.0) {
- // both eigenvalues are the same, no difference which one to take
- coeff.m_w_plus << -l_1, 0.0, 1.0, 0.0;
- coeff.m_w_min << 0.0, -l_1, 0.0, 1.0;
- return;
- }
-
- // First basis vector that has no upward going wave amplitude
- coeff.m_w_min(0) = (b.x() - I * b.y()) * (l_1 - l_2) / 2.0 / b_mag;
- coeff.m_w_min(1) = b.z() * (l_2 - l_1) / 2.0 / b_mag - (l_1 + l_2) / 2.0;
- coeff.m_w_min(2) = 0.0;
- coeff.m_w_min(3) = 1.0;
-
- // Second basis vector that has no upward going wave amplitude
- coeff.m_w_plus(0) = -(l_1 + l_2) / 2.0 - b.z() / (l_1 + l_2);
- coeff.m_w_plus(1) = (b.x() + I * b.y()) * (l_1 - l_2) / 2.0 / b_mag;
- coeff.m_w_plus(2) = 1.0;
- coeff.m_w_plus(3) = 0.0;
+ return result;
}
-void SpecularMagneticStrategy::propagateBackwardsForwards(
- std::vector<MatrixRTCoefficients_v2>& coeff, const std::vector<Slice>& slices) {
- const int size = static_cast<int>(coeff.size());
- std::vector<Eigen::Matrix2cd> SMatrices(coeff.size());
- std::vector<complex_t> Normalization(coeff.size());
-
- for (int index = size - 2; index >= 0; --index) {
- const size_t i = static_cast<size_t>(index);
- const double t = slices[i].thickness();
- const auto kz = coeff[i].getKz();
- const Eigen::Matrix4cd l = coeff[i].R1 * GetImExponential(kz(0) * t)
- + coeff[i].T1 * GetImExponential(-kz(0) * t)
- + coeff[i].R2 * GetImExponential(kz(1) * t)
- + coeff[i].T2 * GetImExponential(-kz(1) * t);
- coeff[i].m_w_plus = l * coeff[i + 1].m_w_plus;
- coeff[i].m_w_min = l * coeff[i + 1].m_w_min;
-
- // rotate and normalize polarization
- const auto Snorm = findNormalizationCoefficients(coeff[i]);
- auto S = Snorm.first;
- auto norm = Snorm.second;
-
+void SpecularMagneticStrategy::calculateUpwards(std::vector<MatrixRTCoefficients>& coeff,
+ const std::vector<Slice>& slices) const {
+ const auto N = slices.size();
+ std::vector<Eigen::Matrix2cd> SMatrices(N - 1);
+ std::vector<complex_t> Normalization(N - 1);
+
+ // bottom boundary condition
+ coeff.back().m_T = Eigen::Matrix2cd::Identity();
+ coeff.back().m_R = Eigen::Matrix2cd::Zero();
+
+ for (int i = N - 2; i >= 0; --i) {
+ double sigma = 0.;
+ if (const auto roughness = GetBottomRoughness(slices, i))
+ sigma = roughness->getSigma();
+
+ // compute the 2x2 submatrices in the write-up denoted as P+, P- and delta
+ const auto mpmm = computeBackwardsSubmatrices(coeff[i], coeff[i + 1], sigma);
+ const Eigen::Matrix2cd mp = mpmm.first;
+ const Eigen::Matrix2cd mm = mpmm.second;
+
+ const Eigen::Matrix2cd delta = coeff[i].computeDeltaMatrix(slices[i].thickness());
+
+ // compute the rotation matrix
+ Eigen::Matrix2cd S, Si;
+ Si = mp + mm * coeff[i + 1].m_R;
+ S << Si(1, 1), -Si(0, 1), -Si(1, 0), Si(0, 0);
+ const complex_t norm = S(0, 0) * S(1, 1) - S(0, 1) * S(1, 0);
+ S = S * delta;
+
+ // store the rotation matrix and normalization constant in order to rotate
+ // the coefficients for all lower slices at the end of the computation
SMatrices[i] = S;
Normalization[i] = norm;
- const complex_t a_plus = S(0, 0) / norm;
- const complex_t b_plus = S(1, 0) / norm;
- const complex_t a_min = S(0, 1) / norm;
- const complex_t b_min = S(1, 1) / norm;
-
- const Eigen::Vector4cd w_plus = a_plus * coeff[i].m_w_plus + b_plus * coeff[i].m_w_min;
- const Eigen::Vector4cd w_min = a_min * coeff[i].m_w_plus + b_min * coeff[i].m_w_min;
+ // compute the reflection matrix and
+ // rotate the polarization such that we have pure incoming states (T = I)
+ S /= norm;
- coeff[i].m_w_plus = std::move(w_plus);
- coeff[i].m_w_min = std::move(w_min);
+ // T is always equal to the identity at this point, no need to store
+ coeff[i].m_R = delta * (mm + mp * coeff[i + 1].m_R) * S;
}
+ // now correct all amplitudes in forward direction by dividing with the remaining
+ // normalization constants. In addition rotate the polarization by the amount
+ // that was rotated above the current interface
+ // if the normalization overflows, all amplitudes below that point are set to zero
complex_t dumpingFactor = 1;
Eigen::Matrix2cd S = Eigen::Matrix2cd::Identity();
- for (size_t i = 1; i < coeff.size(); ++i) {
+ for (size_t i = 1; i < N; ++i) {
dumpingFactor = dumpingFactor * Normalization[i - 1];
S = SMatrices[i - 1] * S;
if (std::isinf(std::norm(dumpingFactor))) {
- // not entirely sure, whether this is the correct edge case
- std::for_each(coeff.begin() + i, coeff.end(),
- [](auto& coeff) { setNoTransmission(coeff); });
+ std::for_each(coeff.begin() + i, coeff.end(), [](auto& coeff) {
+ coeff.m_T = Eigen::Matrix2cd::Zero();
+ coeff.m_R = Eigen::Matrix2cd::Zero();
+ });
break;
}
- const complex_t a_plus = S(0, 0) / dumpingFactor;
- const complex_t b_plus = S(1, 0) / dumpingFactor;
- const complex_t a_min = S(0, 1) / dumpingFactor;
- const complex_t b_min = S(1, 1) / dumpingFactor;
-
- Eigen::Vector4cd w_plus = a_plus * coeff[i].m_w_plus + b_plus * coeff[i].m_w_min;
- Eigen::Vector4cd w_min = a_min * coeff[i].m_w_plus + b_min * coeff[i].m_w_min;
-
- coeff[i].m_w_plus = std::move(w_plus);
- coeff[i].m_w_min = std::move(w_min);
+ coeff[i].m_T = S / dumpingFactor; // T * S omitted, since T is always I
+ coeff[i].m_R *= S / dumpingFactor;
}
}
-std::pair<Eigen::Matrix2cd, complex_t>
-SpecularMagneticStrategy::findNormalizationCoefficients(const MatrixRTCoefficients_v2& coeff) {
- const Eigen::Vector2cd Ta = coeff.T1plus() + coeff.T2plus();
- const Eigen::Vector2cd Tb = coeff.T1min() + coeff.T2min();
-
- Eigen::Matrix2cd S;
- S << Ta(0), Tb(0), Ta(1), Tb(1);
-
- Eigen::Matrix2cd SInverse;
- SInverse << S(1, 1), -S(0, 1), -S(1, 0), S(0, 0);
- const complex_t d1 = S(1, 1) - S(0, 1);
- const complex_t d2 = S(1, 0) - S(0, 0);
- const complex_t denominator = S(0, 0) * d1 - d2 * S(0, 1);
-
- return {SInverse, denominator};
-}
-
namespace {
-kvector_t magneticImpact(kvector_t B_field) {
- return -magnetic_prefactor * B_field;
+double magneticSLD(kvector_t B_field) {
+ return magnetic_prefactor * B_field.mag();
}
-Eigen::Vector2cd eigenvalues(complex_t kz, double b_mag) {
+Eigen::Vector2cd eigenvalues(complex_t kz, double magnetic_SLD) {
const complex_t a = kz * kz;
- return {I * std::sqrt(a + b_mag), I * std::sqrt(a - b_mag)};
+ return {std::sqrt(a - 4. * M_PI * magnetic_SLD), std::sqrt(a + 4. * M_PI * magnetic_SLD)};
}
Eigen::Vector2cd checkForUnderflow(const Eigen::Vector2cd& eigenvs) {
@@ -265,9 +174,10 @@ Eigen::Vector2cd checkForUnderflow(const Eigen::Vector2cd& eigenvs) {
return {lambda(eigenvs(0)), lambda(eigenvs(1))};
}
-complex_t GetImExponential(complex_t exponent) {
- if (exponent.imag() > -std::log(std::numeric_limits<double>::min()))
- return 0.0;
- return std::exp(I * exponent);
+const LayerRoughness* GetBottomRoughness(const std::vector<Slice>& slices,
+ const size_t slice_index) {
+ if (slice_index + 1 < slices.size())
+ return slices[slice_index + 1].topRoughness();
+ return nullptr;
}
} // namespace
diff --git a/Sample/Specular/SpecularMagneticStrategy.h b/Sample/Specular/SpecularMagneticStrategy.h
index 31a5c12d73652606c60f857d247a813494b713b1..4906c817a347788beb1d0d649d3f152bf81c8877 100644
--- a/Sample/Specular/SpecularMagneticStrategy.h
+++ b/Sample/Specular/SpecularMagneticStrategy.h
@@ -7,7 +7,7 @@
//!
//! @homepage http://www.bornagainproject.org
//! @license GNU General Public License v3 or higher (see COPYING)
-//! @copyright Forschungszentrum Jülich GmbH 2018
+//! @copyright Forschungszentrum Jülich GmbH 2020
//! @authors Scientific Computing Group at MLZ (see CITATION, AUTHORS)
//
// ************************************************************************************************
@@ -15,25 +15,33 @@
#ifndef BORNAGAIN_SAMPLE_SPECULAR_SPECULARMAGNETICSTRATEGY_H
#define BORNAGAIN_SAMPLE_SPECULAR_SPECULARMAGNETICSTRATEGY_H
-#include "Sample/RT/MatrixRTCoefficients_v2.h"
+#include "Sample/RT/MatrixRTCoefficients.h"
#include "Sample/Specular/ISpecularStrategy.h"
#include <memory>
#include <vector>
class Slice;
-//! Implements the magnetic Fresnel computation without roughness
+//! Implements the magnetic Fresnel computation with Nevot-Croce roughness
//!
-//! Implements the matrix formalism for the calculation of wave amplitudes of
-//! the coherent wave solution in a multilayer with magnetization.
-//! For a detailed description see internal document "Polarized Specular Reflectometry"
+//! Implements the transfer matrix formalism for the calculation of wave
+//! amplitudes of the coherent wave solution in a multilayer with magnetization.
+//! For a description, see internal
+//! document "Polarized Implementation of the Transfer Matrix Method"
//!
//! @ingroup algorithms_internal
class SpecularMagneticStrategy : public ISpecularStrategy {
public:
+ // TODO remove once external test code is not needed anmyore
+ // for the moment i need them!
+ using coefficient_type = MatrixRTCoefficients;
+ using coefficient_pointer_type = std::unique_ptr<const coefficient_type>;
+ using coeffs_t = std::vector<coefficient_pointer_type>;
+
//! Computes refraction angle reflection/transmission coefficients
//! for given sliced multilayer and wavevector k
- ISpecularStrategy::coeffs_t Execute(const std::vector<Slice>& slices, const kvector_t& k) const;
+ ISpecularStrategy::coeffs_t Execute(const std::vector<Slice>& slices,
+ const kvector_t& k) const;
//! Computes refraction angle reflection/transmission coefficients
//! for given sliced multilayer and a set of kz projections corresponding to each slice
@@ -41,30 +49,15 @@ public:
const std::vector<complex_t>& kz) const;
private:
- static std::vector<MatrixRTCoefficients_v2> computeTR(const std::vector<Slice>& slices,
- const std::vector<complex_t>& kzs);
-
- //! Computes frobenius matrices for multilayer solution
- static void calculateTR(MatrixRTCoefficients_v2& coeff);
- static void calculateZeroFieldTR(MatrixRTCoefficients_v2& coeff);
-
- static void setNoTransmission(MatrixRTCoefficients_v2& coeff);
-
- //! initializes reflectionless bottom boundary condition.
- static void nullifyBottomReflection(MatrixRTCoefficients_v2& coeff);
+ std::vector<MatrixRTCoefficients> computeTR(const std::vector<Slice>& slices,
+ const std::vector<complex_t>& kzs) const;
- //! Propagates boundary conditions from the bottom to the top of the layer stack.
- //! Used to compute boundary conditions from the bottom one (with nullified reflection)
- //! simultaneously propagates amplitudes forward again
- //! Due to the use of temporary objects this is combined into one function now
- static void propagateBackwardsForwards(std::vector<MatrixRTCoefficients_v2>& coeff,
- const std::vector<Slice>& slices);
+ virtual std::pair<Eigen::Matrix2cd, Eigen::Matrix2cd>
+ computeBackwardsSubmatrices(const MatrixRTCoefficients& coeff_i,
+ const MatrixRTCoefficients& coeff_i1, double sigma) const = 0;
- //! finds linear coefficients for normalizing transmitted wave to unity.
- //! The left column of the returned matrix corresponds to the coefficients for pure spin-up
- //! wave, while the right column - to the coefficients for the spin-down one.
- static std::pair<Eigen::Matrix2cd, complex_t>
- findNormalizationCoefficients(const MatrixRTCoefficients_v2& coeff);
+ void calculateUpwards(std::vector<MatrixRTCoefficients>& coeff,
+ const std::vector<Slice>& slices) const;
};
#endif // BORNAGAIN_SAMPLE_SPECULAR_SPECULARMAGNETICSTRATEGY_H
diff --git a/Sample/Specular/SpecularMagneticNewTanhStrategy.cpp b/Sample/Specular/SpecularMagneticTanhStrategy.cpp
similarity index 86%
rename from Sample/Specular/SpecularMagneticNewTanhStrategy.cpp
rename to Sample/Specular/SpecularMagneticTanhStrategy.cpp
index 79708afac0ff317c4570834079e18dc2dd7810bc..6ddc26a6b29f2b01335ad9947af7070f1922a750 100644
--- a/Sample/Specular/SpecularMagneticNewTanhStrategy.cpp
+++ b/Sample/Specular/SpecularMagneticTanhStrategy.cpp
@@ -2,8 +2,8 @@
//
// BornAgain: simulate and fit scattering at grazing incidence
//
-//! @file Sample/Specular/SpecularMagneticNewTanhStrategy.cpp
-//! @brief Implements class SpecularMagneticNewTanhStrategy.
+//! @file Sample/Specular/SpecularMagneticTanhStrategy.cpp
+//! @brief Implements class SpecularMagneticTanhStrategy.
//!
//! @homepage http://www.bornagainproject.org
//! @license GNU General Public License v3 or higher (see COPYING)
@@ -12,7 +12,7 @@
//
// ************************************************************************************************
-#include "Sample/Specular/SpecularMagneticNewTanhStrategy.h"
+#include "Sample/Specular/SpecularMagneticTanhStrategy.h"
#include "Base/Math/Constants.h"
#include "Base/Math/Functions.h"
@@ -21,7 +21,7 @@ const double pi2_15 = std::pow(M_PI_2, 1.5);
} // namespace
Eigen::Matrix2cd
-SpecularMagneticNewTanhStrategy::computeRoughnessMatrix(const MatrixRTCoefficients_v3& coeff,
+SpecularMagneticTanhStrategy::computeRoughnessMatrix(const MatrixRTCoefficients& coeff,
double sigma, bool inverse) const {
if (sigma < 10 * std::numeric_limits<double>::epsilon())
return Eigen::Matrix2cd{Eigen::Matrix2cd::Identity()};
@@ -60,8 +60,8 @@ SpecularMagneticNewTanhStrategy::computeRoughnessMatrix(const MatrixRTCoefficien
}
std::pair<Eigen::Matrix2cd, Eigen::Matrix2cd>
-SpecularMagneticNewTanhStrategy::computeBackwardsSubmatrices(
- const MatrixRTCoefficients_v3& coeff_i, const MatrixRTCoefficients_v3& coeff_i1,
+SpecularMagneticTanhStrategy::computeBackwardsSubmatrices(
+ const MatrixRTCoefficients& coeff_i, const MatrixRTCoefficients& coeff_i1,
double sigma) const {
Eigen::Matrix2cd R{Eigen::Matrix2cd::Identity()};
Eigen::Matrix2cd RInv{Eigen::Matrix2cd::Identity()};
diff --git a/Sample/Specular/SpecularMagneticNewTanhStrategy.h b/Sample/Specular/SpecularMagneticTanhStrategy.h
similarity index 62%
rename from Sample/Specular/SpecularMagneticNewTanhStrategy.h
rename to Sample/Specular/SpecularMagneticTanhStrategy.h
index 9bcde6f1ab34f9fe0fa8b925c9e7749947e9ab0c..cb10b10636b46930554656c8183a99b63a70b0c6 100644
--- a/Sample/Specular/SpecularMagneticNewTanhStrategy.h
+++ b/Sample/Specular/SpecularMagneticTanhStrategy.h
@@ -2,8 +2,8 @@
//
// BornAgain: simulate and fit scattering at grazing incidence
//
-//! @file Sample/Specular/SpecularMagneticNewTanhStrategy.h
-//! @brief Defines class SpecularMagneticNewTanhStrategy.
+//! @file Sample/Specular/SpecularMagneticTanhStrategy.h
+//! @brief Defines class SpecularMagneticTanhStrategy.
//!
//! @homepage http://www.bornagainproject.org
//! @license GNU General Public License v3 or higher (see COPYING)
@@ -12,10 +12,10 @@
//
// ************************************************************************************************
-#ifndef BORNAGAIN_SAMPLE_SPECULAR_SPECULARMAGNETICNEWTANHSTRATEGY_H
-#define BORNAGAIN_SAMPLE_SPECULAR_SPECULARMAGNETICNEWTANHSTRATEGY_H
+#ifndef BORNAGAIN_SAMPLE_SPECULAR_SPECULARMAGNETICTANHSTRATEGY_H
+#define BORNAGAIN_SAMPLE_SPECULAR_SPECULARMAGNETICTANHSTRATEGY_H
-#include "Sample/Specular/SpecularMagneticNewStrategy.h"
+#include "Sample/Specular/SpecularMagneticStrategy.h"
//! Implements the magnetic Fresnel computation with the analytical Tanh roughness
//!
@@ -25,14 +25,14 @@
//! document "Polarized Implementation of the Transfer Matrix Method"
//!
//! @ingroup algorithms_internal
-class SpecularMagneticNewTanhStrategy : public SpecularMagneticNewStrategy {
+class SpecularMagneticTanhStrategy : public SpecularMagneticStrategy {
private:
virtual std::pair<Eigen::Matrix2cd, Eigen::Matrix2cd>
- computeBackwardsSubmatrices(const MatrixRTCoefficients_v3& coeff_i,
- const MatrixRTCoefficients_v3& coeff_i1, double sigma) const;
+ computeBackwardsSubmatrices(const MatrixRTCoefficients& coeff_i,
+ const MatrixRTCoefficients& coeff_i1, double sigma) const;
- Eigen::Matrix2cd computeRoughnessMatrix(const MatrixRTCoefficients_v3& coeff, double sigma,
+ Eigen::Matrix2cd computeRoughnessMatrix(const MatrixRTCoefficients& coeff, double sigma,
bool inverse = false) const;
};
-#endif // BORNAGAIN_SAMPLE_SPECULAR_SPECULARMAGNETICNEWTANHSTRATEGY_H
+#endif // BORNAGAIN_SAMPLE_SPECULAR_SPECULARMAGNETICTANHSTRATEGY_H
diff --git a/Sample/Specular/SpecularScalarStrategy.h b/Sample/Specular/SpecularScalarStrategy.h
index eda61cef769dcc30b9e9ddcba0bcc2fb47afb789..15d8ae24e7b64f55e31e0fd8fe4f54f987bce2b6 100644
--- a/Sample/Specular/SpecularScalarStrategy.h
+++ b/Sample/Specular/SpecularScalarStrategy.h
@@ -32,6 +32,12 @@ class Slice;
//! @ingroup algorithms_internal
class SpecularScalarStrategy : public ISpecularStrategy {
public:
+ // TODO remove once external test code is not needed anmyore
+ // for the moment i need them!
+ using coefficient_type = ScalarRTCoefficients;
+ using coefficient_pointer_type = std::unique_ptr<const coefficient_type>;
+ using coeffs_t = std::vector<coefficient_pointer_type>;
+
//! Computes refraction angles and transmission/reflection coefficients
//! for given coherent wave propagation in a multilayer.
virtual ISpecularStrategy::coeffs_t Execute(const std::vector<Slice>& slices,
diff --git a/Sample/Specular/SpecularStrategyBuilder.cpp b/Sample/Specular/SpecularStrategyBuilder.cpp
index b79d0decf6fbe99dff11a75683ab566a87a7bca8..767d53588522e50607c84bfdaf1edf59f73a8838 100644
--- a/Sample/Specular/SpecularStrategyBuilder.cpp
+++ b/Sample/Specular/SpecularStrategyBuilder.cpp
@@ -14,8 +14,8 @@
#include "Sample/Specular/SpecularStrategyBuilder.h"
#include "Sample/Multilayer/MultiLayerUtils.h"
-#include "Sample/Specular/SpecularMagneticNewNCStrategy.h"
-#include "Sample/Specular/SpecularMagneticNewTanhStrategy.h"
+#include "Sample/Specular/SpecularMagneticNCStrategy.h"
+#include "Sample/Specular/SpecularMagneticTanhStrategy.h"
#include "Sample/Specular/SpecularScalarNCStrategy.h"
#include "Sample/Specular/SpecularScalarTanhStrategy.h"
@@ -26,11 +26,11 @@ std::unique_ptr<ISpecularStrategy> SpecularStrategyBuilder::build(const MultiLay
if (magnetic) {
if (roughnessModel == RoughnessModel::TANH || roughnessModel == RoughnessModel::DEFAULT) {
- return std::make_unique<SpecularMagneticNewTanhStrategy>();
+ return std::make_unique<SpecularMagneticTanhStrategy>();
} else if (roughnessModel == RoughnessModel::NEVOT_CROCE) {
- return std::make_unique<SpecularMagneticNewNCStrategy>();
+ return std::make_unique<SpecularMagneticNCStrategy>();
} else
throw std::logic_error("Invalid roughness model");
diff --git a/Tests/GTestWrapper/google_test.h b/Tests/GTestWrapper/google_test.h
index 1cdb8e03690b95278d089ec39b7213edee426a82..b48f1d4f6dde7968d008458103086f9ddc07393f 100644
--- a/Tests/GTestWrapper/google_test.h
+++ b/Tests/GTestWrapper/google_test.h
@@ -38,4 +38,16 @@
EXPECT_NEAR_COMPLEX(v1(0), v2(0), eps); \
EXPECT_NEAR_COMPLEX(v1(1), v2(1), eps);
+#define EXPECT_NEAR_MATRIXCOEFF(c1, c2, eps) \
+ EXPECT_NEAR_VECTOR2CD(c1->T1plus(), c2->T1plus(), eps); \
+ EXPECT_NEAR_VECTOR2CD(c1->R1plus(), c2->R1plus(), eps); \
+ EXPECT_NEAR_VECTOR2CD(c1->T2plus(), c2->T2plus(), eps); \
+ EXPECT_NEAR_VECTOR2CD(c1->R2plus(), c2->R2plus(), eps); \
+ EXPECT_NEAR_VECTOR2CD(c1->T1min(), c2->T1min(), eps); \
+ EXPECT_NEAR_VECTOR2CD(c1->R1min(), c2->R1min(), eps); \
+ EXPECT_NEAR_VECTOR2CD(c1->T2min(), c2->T2min(), eps); \
+ EXPECT_NEAR_VECTOR2CD(c1->R2min(), c2->R2min(), eps); \
+ EXPECT_NEAR_VECTOR2CD(c1->getKz(), c2->getKz(), eps)
+
+
#endif // BORNAGAIN_TESTS_GTESTWRAPPER_GOOGLE_TEST_H
diff --git a/Tests/UnitTests/Core/Fresnel/MatrixRTCoefficientsTest.cpp b/Tests/UnitTests/Core/Fresnel/MatrixRTCoefficientsTest.cpp
index bcee4702fe57ffae5e474146f4f1ce71b3806bab..ba2843d41185ff3cc07635d0bf5cefd5a399df90 100644
--- a/Tests/UnitTests/Core/Fresnel/MatrixRTCoefficientsTest.cpp
+++ b/Tests/UnitTests/Core/Fresnel/MatrixRTCoefficientsTest.cpp
@@ -3,55 +3,55 @@
class MatrixRTCoefficientsTest : public ::testing::Test {
protected:
- MatrixRTCoefficients mrtcDefault;
+ MatrixRTCoefficients mrtcDefault{-1., {0., 0.}, kvector_t{0.0, 0.0, 0.0}, 0.};
};
TEST_F(MatrixRTCoefficientsTest, T1plus) {
Eigen::Vector2cd vector = mrtcDefault.T1plus();
- EXPECT_EQ(complex_t(0.5, 0.0), vector(0));
+ EXPECT_EQ(0.0, vector(0));
EXPECT_EQ(0.0, vector(1));
}
TEST_F(MatrixRTCoefficientsTest, T1min) {
Eigen::Vector2cd vector = mrtcDefault.T1min();
EXPECT_EQ(0.0, vector(0));
- EXPECT_EQ(complex_t(0.5, 0.0), vector(1));
+ EXPECT_EQ(complex_t(1.0, 0.0), vector(1));
}
TEST_F(MatrixRTCoefficientsTest, T2plus) {
Eigen::Vector2cd vector = mrtcDefault.T2plus();
- EXPECT_EQ(complex_t(0.5, 0.0), vector(0));
+ EXPECT_EQ(complex_t(1.0, 0.0), vector(0));
EXPECT_EQ(0.0, vector(1));
}
TEST_F(MatrixRTCoefficientsTest, T2min) {
Eigen::Vector2cd vector = mrtcDefault.T2min();
EXPECT_EQ(0.0, vector(0));
- EXPECT_EQ(complex_t(0.5, 0.0), vector(1));
+ EXPECT_EQ(0.0, vector(1));
}
TEST_F(MatrixRTCoefficientsTest, R1plus) {
Eigen::Vector2cd vector = mrtcDefault.R1plus();
- EXPECT_EQ(complex_t(-0.5, 0.0), vector(0));
+ EXPECT_EQ(0.0, vector(0));
EXPECT_EQ(0.0, vector(1));
}
TEST_F(MatrixRTCoefficientsTest, R1min) {
Eigen::Vector2cd vector = mrtcDefault.R1min();
EXPECT_EQ(0.0, vector(0));
- EXPECT_EQ(complex_t(-0.5, 0.0), vector(1));
+ EXPECT_EQ(complex_t(-1.0, 0.0), vector(1));
}
TEST_F(MatrixRTCoefficientsTest, R2plus) {
Eigen::Vector2cd vector = mrtcDefault.R2plus();
- EXPECT_EQ(complex_t(-0.5, 0.0), vector(0));
+ EXPECT_EQ(complex_t(-1.0, 0.0), vector(0));
EXPECT_EQ(0.0, vector(1));
}
TEST_F(MatrixRTCoefficientsTest, R2min) {
Eigen::Vector2cd vector = mrtcDefault.R2min();
EXPECT_EQ(0.0, vector(0));
- EXPECT_EQ(complex_t(-0.5, 0.0), vector(1));
+ EXPECT_EQ(0.0, vector(1));
}
TEST_F(MatrixRTCoefficientsTest, getKz) {
diff --git a/Tests/UnitTests/Core/Fresnel/MatrixRTCoefficients_v3Test.cpp b/Tests/UnitTests/Core/Fresnel/MatrixRTCoefficients_v3Test.cpp
deleted file mode 100644
index e237eb7ba99880cb71b9575364b84710308a0501..0000000000000000000000000000000000000000
--- a/Tests/UnitTests/Core/Fresnel/MatrixRTCoefficients_v3Test.cpp
+++ /dev/null
@@ -1,61 +0,0 @@
-#include "Sample/RT/MatrixRTCoefficients_v3.h"
-#include "Tests/GTestWrapper/google_test.h"
-
-class MatrixRTCoefficients_v3Test : public ::testing::Test {
-protected:
- MatrixRTCoefficients_v3 mrtcDefault{-1., {0., 0.}, kvector_t{0.0, 0.0, 0.0}, 0.};
-};
-
-TEST_F(MatrixRTCoefficients_v3Test, T1plus) {
- Eigen::Vector2cd vector = mrtcDefault.T1plus();
- EXPECT_EQ(0.0, vector(0));
- EXPECT_EQ(0.0, vector(1));
-}
-
-TEST_F(MatrixRTCoefficients_v3Test, T1min) {
- Eigen::Vector2cd vector = mrtcDefault.T1min();
- EXPECT_EQ(0.0, vector(0));
- EXPECT_EQ(complex_t(1.0, 0.0), vector(1));
-}
-
-TEST_F(MatrixRTCoefficients_v3Test, T2plus) {
- Eigen::Vector2cd vector = mrtcDefault.T2plus();
- EXPECT_EQ(complex_t(1.0, 0.0), vector(0));
- EXPECT_EQ(0.0, vector(1));
-}
-
-TEST_F(MatrixRTCoefficients_v3Test, T2min) {
- Eigen::Vector2cd vector = mrtcDefault.T2min();
- EXPECT_EQ(0.0, vector(0));
- EXPECT_EQ(0.0, vector(1));
-}
-
-TEST_F(MatrixRTCoefficients_v3Test, R1plus) {
- Eigen::Vector2cd vector = mrtcDefault.R1plus();
- EXPECT_EQ(0.0, vector(0));
- EXPECT_EQ(0.0, vector(1));
-}
-
-TEST_F(MatrixRTCoefficients_v3Test, R1min) {
- Eigen::Vector2cd vector = mrtcDefault.R1min();
- EXPECT_EQ(0.0, vector(0));
- EXPECT_EQ(complex_t(-1.0, 0.0), vector(1));
-}
-
-TEST_F(MatrixRTCoefficients_v3Test, R2plus) {
- Eigen::Vector2cd vector = mrtcDefault.R2plus();
- EXPECT_EQ(complex_t(-1.0, 0.0), vector(0));
- EXPECT_EQ(0.0, vector(1));
-}
-
-TEST_F(MatrixRTCoefficients_v3Test, R2min) {
- Eigen::Vector2cd vector = mrtcDefault.R2min();
- EXPECT_EQ(0.0, vector(0));
- EXPECT_EQ(0.0, vector(1));
-}
-
-TEST_F(MatrixRTCoefficients_v3Test, getKz) {
- Eigen::Vector2cd vector = mrtcDefault.getKz();
- EXPECT_EQ(complex_t(0.0, 0.0), vector(0));
- EXPECT_EQ(complex_t(0.0, 0.0), vector(1));
-}
diff --git a/Tests/UnitTests/Core/Fresnel/SpecularMagneticTest.cpp b/Tests/UnitTests/Core/Fresnel/SpecularMagneticTest.cpp
index bff76f55318ac388cbeca38d3bfc359b64ee6715..9d57f448becfc1511b5fc9cc14a49a1133a8a1fc 100644
--- a/Tests/UnitTests/Core/Fresnel/SpecularMagneticTest.cpp
+++ b/Tests/UnitTests/Core/Fresnel/SpecularMagneticTest.cpp
@@ -1,30 +1,31 @@
#include "Base/Const/Units.h"
+#include "Core/Legacy/SpecularMagneticStrategy_v2.h"
#include "Sample/Material/MaterialFactoryFuncs.h"
#include "Sample/Multilayer/Layer.h"
#include "Sample/Multilayer/MultiLayer.h"
#include "Sample/Processed/ProcessedSample.h"
#include "Sample/RT/SimulationOptions.h"
#include "Sample/Slice/KzComputation.h"
-#include "Sample/Specular/SpecularMagneticNewTanhStrategy.h"
-#include "Sample/Specular/SpecularMagneticStrategy.h"
+#include "Sample/Specular/SpecularMagneticTanhStrategy.h"
#include "Sample/Specular/SpecularScalarTanhStrategy.h"
#include "Tests/GTestWrapper/google_test.h"
#include <utility>
-constexpr double eps = 1e-10;
class SpecularMagneticTest : public ::testing::Test {
protected:
- std::unique_ptr<ProcessedSample> sample_zerofield();
+ auto static constexpr eps = 1.e-10;
+
+ std::unique_ptr<ProcessedSample> sample_zerofield(bool slab);
std::unique_ptr<ProcessedSample> sample_degenerate();
template <typename Strategy> void test_degenerate();
template <typename Strategy>
void testZeroField(const kvector_t& k, const ProcessedSample& sample);
- template <typename Strategy> void testcase_zerofield(std::vector<double>&& angles);
+ template <typename Strategy> void testcase_zerofield(std::vector<double>&& angles, bool slab=false);
};
-template <> void SpecularMagneticTest::test_degenerate<SpecularMagneticNewTanhStrategy>() {
+template <> void SpecularMagneticTest::test_degenerate<SpecularMagneticTanhStrategy>() {
kvector_t v;
Eigen::Vector2cd T1p{0.0, 0.0};
@@ -37,7 +38,7 @@ template <> void SpecularMagneticTest::test_degenerate<SpecularMagneticNewTanhSt
Eigen::Vector2cd R2m{0.0, 0.0};
auto sample = sample_degenerate();
- auto result = std::make_unique<SpecularMagneticNewTanhStrategy>()->Execute(sample->slices(), v);
+ auto result = std::make_unique<SpecularMagneticTanhStrategy>()->Execute(sample->slices(), v);
for (auto& coeff : result) {
EXPECT_NEAR_VECTOR2CD(coeff->T1plus(), T1p, eps);
EXPECT_NEAR_VECTOR2CD(coeff->T2plus(), T2p, eps);
@@ -85,16 +86,23 @@ std::unique_ptr<ProcessedSample> SpecularMagneticTest::sample_degenerate() {
return std::make_unique<ProcessedSample>(mLayer, SimulationOptions());
}
-TEST_F(SpecularMagneticTest, degenerate_new) {
- test_degenerate<SpecularMagneticNewTanhStrategy>();
+TEST_F(SpecularMagneticTest, degenerate_) {
+ test_degenerate<SpecularMagneticTanhStrategy>();
}
-std::unique_ptr<ProcessedSample> SpecularMagneticTest::sample_zerofield() {
+std::unique_ptr<ProcessedSample> SpecularMagneticTest::sample_zerofield(bool slab) {
MultiLayer multi_layer_scalar;
Material substr_material_scalar = HomogeneousMaterial("Substrate", 7e-6, 2e-8);
Layer vacuum_layer(HomogeneousMaterial("Vacuum", 0.0, 0.0));
Layer substr_layer_scalar(substr_material_scalar);
multi_layer_scalar.addLayer(vacuum_layer);
+
+ if(slab) {
+ Material layer_material = HomogeneousMaterial("Layer", 3e-6, 1e-8);
+ Layer layer(layer_material, 10. * Units::nm);
+ multi_layer_scalar.addLayer(layer);
+ }
+
multi_layer_scalar.addLayer(substr_layer_scalar);
SimulationOptions options;
@@ -104,18 +112,43 @@ std::unique_ptr<ProcessedSample> SpecularMagneticTest::sample_zerofield() {
}
template <typename Strategy>
-void SpecularMagneticTest::testcase_zerofield(std::vector<double>&& angles) {
+void SpecularMagneticTest::testcase_zerofield(std::vector<double>&& angles, bool slab) {
for (auto& angle : angles) {
- auto sample = sample_zerofield();
+ auto sample = sample_zerofield(slab);
kvector_t k = vecOfLambdaAlphaPhi(1.0, angle * Units::deg, 0.0);
testZeroField<Strategy>(k, *sample);
}
}
-TEST_F(SpecularMagneticTest, zerofield) {
- testcase_zerofield<SpecularMagneticStrategy>({-0.1, -2.0, -10.0});
+TEST_F(SpecularMagneticTest, zerofield_v2_negative_k) {
+ testcase_zerofield<SpecularMagneticStrategy_v2>({-0.1, -2.0, -10.0});
+}
+
+TEST_F(SpecularMagneticTest, zerofield_v2_positive_k) {
+ testcase_zerofield<SpecularMagneticStrategy_v2>({0.1, 2.0, 10.0});
+}
+
+TEST_F(SpecularMagneticTest, zerofield2_v2_negative_k) {
+ testcase_zerofield<SpecularMagneticStrategy_v2>({-0.1, -2.0, -10.0}, true);
+}
+
+TEST_F(SpecularMagneticTest, zerofield2_v2_positive_k) {
+ testcase_zerofield<SpecularMagneticStrategy_v2>({0.1, 2.0, 10.0}, true);
+}
+
+
+TEST_F(SpecularMagneticTest, zerofield_positive_k) {
+ testcase_zerofield<SpecularMagneticTanhStrategy>({0.0, 1.e-9, 1.e-5, 0.1, 2.0, 10.0});
+}
+
+TEST_F(SpecularMagneticTest, zerofield_negative_k) {
+ testcase_zerofield<SpecularMagneticTanhStrategy>({-0.0, -1.e-9, -1.e-5, -0.1, -2.0, -10.0});
+}
+
+TEST_F(SpecularMagneticTest, zerofield2_positive_k) {
+ testcase_zerofield<SpecularMagneticTanhStrategy>({0.0, 1.e-9, 1.e-5, 0.1, 2.0, 10.0}, true);
}
-TEST_F(SpecularMagneticTest, zerofield_new) {
- testcase_zerofield<SpecularMagneticNewTanhStrategy>({-0.0, -1.e-9, -1.e-5, -0.1, -2.0, -10.0});
+TEST_F(SpecularMagneticTest, zerofield2_negative_k) {
+ testcase_zerofield<SpecularMagneticTanhStrategy>({-0.0, -1.e-9, -1.e-5, -0.1, -2.0, -10.0}, true);
}
diff --git a/Tests/UnitTests/Core/Legacy/MatrixRTCoefficients_v1Test.cpp b/Tests/UnitTests/Core/Legacy/MatrixRTCoefficients_v1Test.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1b8d11b56eb535adb469acfc3dffb5258da8023d
--- /dev/null
+++ b/Tests/UnitTests/Core/Legacy/MatrixRTCoefficients_v1Test.cpp
@@ -0,0 +1,61 @@
+#include "Core/Legacy/MatrixRTCoefficients_v1.h"
+#include "Tests/GTestWrapper/google_test.h"
+
+class MatrixRTCoefficients_v1Test : public ::testing::Test {
+protected:
+ MatrixRTCoefficients_v1 mrtcDefault;
+};
+
+TEST_F(MatrixRTCoefficients_v1Test, T1plus) {
+ Eigen::Vector2cd vector = mrtcDefault.T1plus();
+ EXPECT_EQ(complex_t(0.5, 0.0), vector(0));
+ EXPECT_EQ(0.0, vector(1));
+}
+
+TEST_F(MatrixRTCoefficients_v1Test, T1min) {
+ Eigen::Vector2cd vector = mrtcDefault.T1min();
+ EXPECT_EQ(0.0, vector(0));
+ EXPECT_EQ(complex_t(0.5, 0.0), vector(1));
+}
+
+TEST_F(MatrixRTCoefficients_v1Test, T2plus) {
+ Eigen::Vector2cd vector = mrtcDefault.T2plus();
+ EXPECT_EQ(complex_t(0.5, 0.0), vector(0));
+ EXPECT_EQ(0.0, vector(1));
+}
+
+TEST_F(MatrixRTCoefficients_v1Test, T2min) {
+ Eigen::Vector2cd vector = mrtcDefault.T2min();
+ EXPECT_EQ(0.0, vector(0));
+ EXPECT_EQ(complex_t(0.5, 0.0), vector(1));
+}
+
+TEST_F(MatrixRTCoefficients_v1Test, R1plus) {
+ Eigen::Vector2cd vector = mrtcDefault.R1plus();
+ EXPECT_EQ(complex_t(-0.5, 0.0), vector(0));
+ EXPECT_EQ(0.0, vector(1));
+}
+
+TEST_F(MatrixRTCoefficients_v1Test, R1min) {
+ Eigen::Vector2cd vector = mrtcDefault.R1min();
+ EXPECT_EQ(0.0, vector(0));
+ EXPECT_EQ(complex_t(-0.5, 0.0), vector(1));
+}
+
+TEST_F(MatrixRTCoefficients_v1Test, R2plus) {
+ Eigen::Vector2cd vector = mrtcDefault.R2plus();
+ EXPECT_EQ(complex_t(-0.5, 0.0), vector(0));
+ EXPECT_EQ(0.0, vector(1));
+}
+
+TEST_F(MatrixRTCoefficients_v1Test, R2min) {
+ Eigen::Vector2cd vector = mrtcDefault.R2min();
+ EXPECT_EQ(0.0, vector(0));
+ EXPECT_EQ(complex_t(-0.5, 0.0), vector(1));
+}
+
+TEST_F(MatrixRTCoefficients_v1Test, getKz) {
+ Eigen::Vector2cd vector = mrtcDefault.getKz();
+ EXPECT_EQ(complex_t(0.0, 0.0), vector(0));
+ EXPECT_EQ(complex_t(0.0, 0.0), vector(1));
+}
diff --git a/Tests/UnitTests/Core/Fresnel/MatrixRTCoefficients_v2Test.cpp b/Tests/UnitTests/Core/Legacy/MatrixRTCoefficients_v2Test.cpp
similarity index 97%
rename from Tests/UnitTests/Core/Fresnel/MatrixRTCoefficients_v2Test.cpp
rename to Tests/UnitTests/Core/Legacy/MatrixRTCoefficients_v2Test.cpp
index 0fa5397e86bfe38e71e0354e79db25bd31f8ba97..6934736917de848b4e8d7f8eb33e05671035a704 100644
--- a/Tests/UnitTests/Core/Fresnel/MatrixRTCoefficients_v2Test.cpp
+++ b/Tests/UnitTests/Core/Legacy/MatrixRTCoefficients_v2Test.cpp
@@ -1,4 +1,4 @@
-#include "Sample/RT/MatrixRTCoefficients_v2.h"
+#include "Core/Legacy/MatrixRTCoefficients_v2.h"
#include "Tests/GTestWrapper/google_test.h"
class MatrixRTCoefficients_v2Test : public ::testing::Test {
diff --git a/Tests/UnitTests/Core/Legacy/SpecularMagneticConsistencyTest.cpp b/Tests/UnitTests/Core/Legacy/SpecularMagneticConsistencyTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4b974632db4aada7df0498462a9502d991eb19e4
--- /dev/null
+++ b/Tests/UnitTests/Core/Legacy/SpecularMagneticConsistencyTest.cpp
@@ -0,0 +1,72 @@
+#include "Base/Const/Units.h"
+#include "Core/Legacy/SpecularMagneticStrategy_v2.h"
+#include "Sample/Material/MaterialFactoryFuncs.h"
+#include "Sample/Multilayer/Layer.h"
+#include "Sample/Multilayer/MultiLayer.h"
+#include "Sample/Processed/ProcessedSample.h"
+#include "Sample/RT/SimulationOptions.h"
+#include "Sample/Slice/KzComputation.h"
+#include "Sample/Specular/SpecularMagneticTanhStrategy.h"
+#include "Tests/GTestWrapper/google_test.h"
+
+class SpecularMagneticConsistencyTest : public ::testing::Test
+{
+protected:
+ auto static constexpr eps = 1.e-10;
+
+ std::unique_ptr<ProcessedSample> sample_1();
+
+ template <typename Strategy1, typename Strategy2>
+ void testcase(const std::vector<Slice>& slices, double k);
+};
+
+std::unique_ptr<ProcessedSample> SpecularMagneticConsistencyTest::sample_1()
+{
+ MultiLayer multi_layer_scalar;
+ Material substr_material_scalar = HomogeneousMaterial("Substrate", 7e-6, 2e-8);
+ Material layer_material = HomogeneousMaterial("Layer", 3e-6, 1e-8, kvector_t{1.e7, 0, 1.e7});
+
+ Layer vacuum_layer(HomogeneousMaterial("Vacuum", 0.0, 0.0));
+ Layer substr_layer_scalar(substr_material_scalar);
+ Layer layer(layer_material, 10. * Units::nm);
+
+ multi_layer_scalar.addLayer(vacuum_layer);
+ multi_layer_scalar.addLayer(layer);
+ multi_layer_scalar.addLayer(substr_layer_scalar);
+
+ SimulationOptions options;
+ auto sample_scalar = std::make_unique<ProcessedSample>(multi_layer_scalar, options);
+
+ return sample_scalar;
+}
+
+template <typename Strategy1, typename Strategy2>
+void SpecularMagneticConsistencyTest::testcase(const std::vector<Slice>& slices, double k)
+{
+
+ const auto kz = kvector_t{0., 0., k};
+ const auto coeffs1 = std::make_unique<Strategy1>()->Execute(
+ slices, KzComputation::computeKzFromRefIndices(slices, kz));
+ const auto coeffs2 = std::make_unique<Strategy2>()->Execute(
+ slices, KzComputation::computeKzFromRefIndices(slices, kz));
+
+ for (size_t i = 0; i < coeffs1.size(); ++i) {
+ EXPECT_NEAR_MATRIXCOEFF(coeffs1[i], coeffs2[i], eps);
+ }
+}
+
+TEST_F(SpecularMagneticConsistencyTest, NewOld)
+{
+ using Strategy1 = SpecularMagneticTanhStrategy;
+ using Strategy2 = SpecularMagneticStrategy_v2;
+
+ const auto sample = sample_1();
+ const auto slices = sample->slices();
+
+ const std::vector<double> kzs = {1.e-9, 1.e-5, 0.1, 2.0, 10.0};
+
+ for (const auto& k : kzs) {
+ testcase<Strategy1, Strategy2>(slices, k);
+ testcase<Strategy1, Strategy2>(slices, -k);
+ }
+}
diff --git a/Tests/UnitTests/Core/Fresnel/SpecularMagneticOldTest.cpp b/Tests/UnitTests/Core/Legacy/SpecularMagnetic_v1Test.cpp
similarity index 87%
rename from Tests/UnitTests/Core/Fresnel/SpecularMagneticOldTest.cpp
rename to Tests/UnitTests/Core/Legacy/SpecularMagnetic_v1Test.cpp
index bbcd0cd712bf2f42555d1a3f4c17cbe750a2e0e0..adc1ff579bcb32fb7a5268b9ee918e1e4bba7c3d 100644
--- a/Tests/UnitTests/Core/Fresnel/SpecularMagneticOldTest.cpp
+++ b/Tests/UnitTests/Core/Legacy/SpecularMagnetic_v1Test.cpp
@@ -4,13 +4,13 @@
#include "Sample/Multilayer/MultiLayer.h"
#include "Sample/Processed/ProcessedSample.h"
#include "Sample/RT/SimulationOptions.h"
-#include "Sample/Specular/SpecularMagneticOldStrategy.h"
+#include "Core/Legacy/SpecularMagneticStrategy_v1.h"
#include "Sample/Specular/SpecularScalarTanhStrategy.h"
#include "Tests/GTestWrapper/google_test.h"
-class SpecularMagneticOldTest : public ::testing::Test {};
+class SpecularMagnetic_v1Test : public ::testing::Test {};
-TEST_F(SpecularMagneticOldTest, initial) {
+TEST_F(SpecularMagnetic_v1Test, initial) {
MultiLayer mLayer;
kvector_t v;
@@ -22,10 +22,10 @@ TEST_F(SpecularMagneticOldTest, initial) {
mLayer.addLayer(layer0);
SimulationOptions options;
ProcessedSample sample(mLayer, options);
- std::make_unique<SpecularMagneticOldStrategy>()->Execute(sample.slices(), v);
+ std::make_unique<SpecularMagneticStrategy_v1>()->Execute(sample.slices(), v);
}
-TEST_F(SpecularMagneticOldTest, zerofield) {
+TEST_F(SpecularMagnetic_v1Test, zerofield) {
double eps = 1e-10;
kvector_t substr_field(0.0, 0.0, 0.0);
@@ -61,9 +61,9 @@ TEST_F(SpecularMagneticOldTest, zerofield) {
Eigen::Vector2cd RMS = RTScalar.R1min() + RTScalar.R2min();
auto coeffs_zerofield =
- std::make_unique<SpecularMagneticOldStrategy>()->Execute(sample_zerofield.slices(), k1);
- MatrixRTCoefficients RTMatrix =
- *dynamic_cast<const MatrixRTCoefficients*>(coeffs_zerofield[1].get());
+ std::make_unique<SpecularMagneticStrategy_v1>()->Execute(sample_zerofield.slices(), k1);
+ MatrixRTCoefficients_v1 RTMatrix =
+ *dynamic_cast<const MatrixRTCoefficients_v1*>(coeffs_zerofield[1].get());
Eigen::Vector2cd TPM = RTMatrix.T1plus() + RTMatrix.T2plus();
Eigen::Vector2cd RPM = RTMatrix.R1plus() + RTMatrix.R2plus();
Eigen::Vector2cd TMM = RTMatrix.T1min() + RTMatrix.T2min();
@@ -88,8 +88,8 @@ TEST_F(SpecularMagneticOldTest, zerofield) {
RMS = RTScalar.R1min() + RTScalar.R2min();
coeffs_zerofield =
- std::make_unique<SpecularMagneticOldStrategy>()->Execute(sample_zerofield.slices(), k2);
- RTMatrix = *dynamic_cast<const MatrixRTCoefficients*>(coeffs_zerofield[1].get());
+ std::make_unique<SpecularMagneticStrategy_v1>()->Execute(sample_zerofield.slices(), k2);
+ RTMatrix = *dynamic_cast<const MatrixRTCoefficients_v1*>(coeffs_zerofield[1].get());
TPM = RTMatrix.T1plus() + RTMatrix.T2plus();
RPM = RTMatrix.R1plus() + RTMatrix.R2plus();
TMM = RTMatrix.T1min() + RTMatrix.T2min();
@@ -114,8 +114,8 @@ TEST_F(SpecularMagneticOldTest, zerofield) {
RMS = RTScalar.R1min() + RTScalar.R2min();
coeffs_zerofield =
- std::make_unique<SpecularMagneticOldStrategy>()->Execute(sample_zerofield.slices(), k3);
- RTMatrix = *dynamic_cast<const MatrixRTCoefficients*>(coeffs_zerofield[1].get());
+ std::make_unique<SpecularMagneticStrategy_v1>()->Execute(sample_zerofield.slices(), k3);
+ RTMatrix = *dynamic_cast<const MatrixRTCoefficients_v1*>(coeffs_zerofield[1].get());
TPM = RTMatrix.T1plus() + RTMatrix.T2plus();
RPM = RTMatrix.R1plus() + RTMatrix.R2plus();
TMM = RTMatrix.T1min() + RTMatrix.T2min();