From 7d5cb8fb3332e1a7f4a820324db988389b759c47 Mon Sep 17 00:00:00 2001
From: Randolf Beerwerth <r.beerwerth@fz-juelich.de>
Date: Fri, 4 Sep 2020 16:19:28 +0200
Subject: [PATCH] Provide stabilized computation of magnetic Fresnel
computation
---
Core/Multilayer/SpecularMagneticStrategy.cpp | 78 ++++++++++++++------
Core/Multilayer/SpecularMagneticStrategy.h | 18 ++---
2 files changed, 63 insertions(+), 33 deletions(-)
diff --git a/Core/Multilayer/SpecularMagneticStrategy.cpp b/Core/Multilayer/SpecularMagneticStrategy.cpp
index f58c2bb910d..97f032ec451 100644
--- a/Core/Multilayer/SpecularMagneticStrategy.cpp
+++ b/Core/Multilayer/SpecularMagneticStrategy.cpp
@@ -77,8 +77,7 @@ SpecularMagneticStrategy::computeTR(const std::vector<Slice>& slices,
std::for_each(result.begin(), result.end(), [](auto& coeff) { calculateTR(coeff); });
nullifyBottomReflection(result.back());
- propagateBackwards(result, slices);
- propagateForwards(result, findNormalizationCoefficients(result.front()));
+ propagateBackwardsForwards(result, slices);
return result;
}
@@ -176,10 +175,13 @@ void SpecularMagneticStrategy::nullifyBottomReflection(MatrixRTCoefficients_v2&
coeff.m_w_plus(3) = 0.0;
}
-void SpecularMagneticStrategy::propagateBackwards(std::vector<MatrixRTCoefficients_v2>& coeff,
- const std::vector<Slice>& slices)
+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();
@@ -190,10 +192,54 @@ void SpecularMagneticStrategy::propagateBackwards(std::vector<MatrixRTCoefficien
+ 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
+ auto r = findNormalizationCoefficients(coeff[i]);
+ auto S = std::get<0>(r);
+ auto norm = std::get<1>(r);
+
+ 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;
+
+ 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);
+ }
+
+ auto dumpingFactor = complex_t(1, 0);
+ 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);
}
}
-Eigen::Matrix2cd
+std::pair<Eigen::Matrix2cd, complex_t>
SpecularMagneticStrategy::findNormalizationCoefficients(const MatrixRTCoefficients_v2& coeff)
{
const Eigen::Vector2cd Ta = coeff.T1plus() + coeff.T2plus();
@@ -204,25 +250,11 @@ SpecularMagneticStrategy::findNormalizationCoefficients(const MatrixRTCoefficien
Eigen::Matrix2cd result;
result << S(1, 1), -S(0, 1), -S(1, 0), S(0, 0);
- result /= S(0, 0) * S(1, 1) - S(1, 0) * S(0, 1);
-
- return result;
-}
+ auto d1 = S(1, 1) - S(0, 1);
+ auto d2 = S(1, 0) - S(0, 0);
+ auto denom = S(0, 0) * d1 - d2 * S(0, 1);
-void SpecularMagneticStrategy::propagateForwards(std::vector<MatrixRTCoefficients_v2>& coeff,
- const Eigen::Matrix2cd& weights)
-{
- const complex_t a_plus = weights(0, 0);
- const complex_t b_plus = weights(1, 0);
- const complex_t a_min = weights(0, 1);
- const complex_t b_min = weights(1, 1);
-
- for (auto& term : coeff) {
- Eigen::Vector4cd w_plus = a_plus * term.m_w_plus + b_plus * term.m_w_min;
- Eigen::Vector4cd w_min = a_min * term.m_w_plus + b_min * term.m_w_min;
- term.m_w_plus = std::move(w_plus);
- term.m_w_min = std::move(w_min);
- }
+ return {result, denom};
}
namespace
diff --git a/Core/Multilayer/SpecularMagneticStrategy.h b/Core/Multilayer/SpecularMagneticStrategy.h
index 2ae617a4fe0..c93828d19c5 100644
--- a/Core/Multilayer/SpecularMagneticStrategy.h
+++ b/Core/Multilayer/SpecularMagneticStrategy.h
@@ -33,9 +33,9 @@ class Slice;
class BA_CORE_API_ SpecularMagneticStrategy : public ISpecularStrategy
{
public:
- using coefficient_type = MatrixRTCoefficients_v2;
+ using coefficient_type = MatrixRTCoefficients_v2;
using coefficient_pointer_type = std::unique_ptr<const coefficient_type>;
- using coeffs_t = std::vector<coefficient_pointer_type>;
+ using coeffs_t = std::vector<coefficient_pointer_type>;
//! Computes refraction angle reflection/transmission coefficients
//! for given sliced multilayer and wavevector k
@@ -61,18 +61,16 @@ private:
//! 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)
- static void propagateBackwards(std::vector<MatrixRTCoefficients_v2>& coeff,
- const std::vector<Slice>& slices);
+ //! 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 Eigen::Matrix2cd findNormalizationCoefficients(const MatrixRTCoefficients_v2& coeff);
-
- //! makes a linear combination of boundary conditions with using the given weights for each
- //! coefficient in the vector.
- static void propagateForwards(std::vector<MatrixRTCoefficients_v2>& coeff,
- const Eigen::Matrix2cd& weights);
+ static std::pair<Eigen::Matrix2cd, complex_t>
+ findNormalizationCoefficients(const MatrixRTCoefficients_v2& coeff);
};
#endif // BORNAGAIN_CORE_MULTILAYER_SPECULARMAGNETICSTRATEGY_H
--
GitLab