Provide stabilized computation of magnetic Fresnel computation

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

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