From a3a4bb59bc50c2fe415fd6aab4c88dd2b055899b Mon Sep 17 00:00:00 2001
From: Randolf Beerwerth <r.beerwerth@fz-juelich.de>
Date: Tue, 22 Sep 2020 12:12:17 +0200
Subject: [PATCH] Correctly implement k = 0 corner case
---
.../SpecularMagneticNewStrategy.cpp | 22 ++++++++++++-------
Core/RT/MatrixRTCoefficients_v3.cpp | 7 +++---
2 files changed, 18 insertions(+), 11 deletions(-)
diff --git a/Core/Multilayer/SpecularMagneticNewStrategy.cpp b/Core/Multilayer/SpecularMagneticNewStrategy.cpp
index e6ace084c63..6dac8e0bf6e 100644
--- a/Core/Multilayer/SpecularMagneticNewStrategy.cpp
+++ b/Core/Multilayer/SpecularMagneticNewStrategy.cpp
@@ -56,6 +56,8 @@ 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.");
@@ -63,9 +65,9 @@ SpecularMagneticNewStrategy::computeTR(const std::vector<Slice>& slices,
return {};
std::vector<MatrixRTCoefficients_v3> result;
- result.reserve(slices.size());
+ result.reserve(N);
- const double kz_sign = kzs.front().real() > 0.0 ? 1.0 : -1.0; // save sign to restore it later
+ 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);
@@ -76,12 +78,16 @@ SpecularMagneticNewStrategy::computeTR(const std::vector<Slice>& slices,
B.mag() > eps ? B / B.mag() : kvector_t{0.0, 0.0, 0.0}, magnetic_SLD);
}
- if (std::abs(kzs[0]) < std::sqrt(eps)) {
- for (size_t i = 0, size = slices.size(); i < size; ++i) {
- if (std::abs(kzs[i]) >= std::sqrt(eps)) {
- result[i].m_T << 0, 0, 0, 0;
- result[i].m_R << 0, 0, 0, 0;
- }
+ if(N == 1){
+ result[0].m_T = Eigen::Matrix2cd::Identity();
+ result[0].m_R = Eigen::Matrix2cd::Zero();
+ return result;
+ }else if(kzs[0] == 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;
}
diff --git a/Core/RT/MatrixRTCoefficients_v3.cpp b/Core/RT/MatrixRTCoefficients_v3.cpp
index 18b4ec09c8c..b72cfead875 100644
--- a/Core/RT/MatrixRTCoefficients_v3.cpp
+++ b/Core/RT/MatrixRTCoefficients_v3.cpp
@@ -54,9 +54,9 @@ Eigen::Matrix2cd MatrixRTCoefficients_v3::TransformationMatrix(complex_t eigenva
} else if (m_b.mag() < eps && eigenvalue != 0.)
return exp2;
- // result = Eigen::Matrix2cd(Eigen::DiagonalMatrix<complex_t, 2>({0.5, 0.5}));
+
else if (m_b.mag() < eps && eigenvalue == 0.)
- return Eigen::Matrix2cd(Eigen::DiagonalMatrix<complex_t, 2>({0.5, 0.5}));
+ return exp2;
throw std::runtime_error("Broken magnetic field vector");
}
@@ -145,7 +145,8 @@ Eigen::Matrix2cd MatrixRTCoefficients_v3::computeInverseP() const
const complex_t beta = m_lambda(1) - m_lambda(0);
if (std::abs(alpha * alpha - beta * beta) < eps)
- throw std::runtime_error("Singular p_m");
+// throw std::runtime_error("Singular p_m");
+ return Eigen::Matrix2cd::Zero();
Eigen::Matrix2cd result = pMatrixHelper(-1.);
result *= 2. / (alpha * alpha - beta * beta);
--
GitLab