Write scalar kernel as proper Parrat iteration

Core/Multilayer/SpecularScalarNCStrategy.cpp

@@ -15,7 +15,8 @@
 #include "Core/Multilayer/SpecularScalarNCStrategy.h"
 #include <Eigen/Dense>
 
-Eigen::Vector2cd SpecularScalarNCStrategy::transition(complex_t kzi, complex_t kzi1, double sigma,
+std::pair<complex_t, complex_t>
+SpecularScalarNCStrategy::transition(complex_t kzi, complex_t kzi1, double sigma,
                                                       double thickness,
                                                       const Eigen::Vector2cd& t_r1) const
 {
@@ -31,8 +32,8 @@ Eigen::Vector2cd SpecularScalarNCStrategy::transition(complex_t kzi, complex_t k
     const complex_t a00 = 0.5 * (1. + kz_ratio) * roughness_diff;
     const complex_t a01 = 0.5 * (1. - kz_ratio) * roughness_sum;
 
-    Eigen::Vector2cd result;
-    result << (a00 * t_r1(0) + a01 * t_r1(1)) / phase_shift,
-        (a01 * t_r1(0) + a00 * t_r1(1)) * phase_shift;
-    return result;
+//    Eigen::Vector2cd result;
+//    result << (a00 * t_r1(0) + a01 * t_r1(1)) / phase_shift,
+//        (a01 * t_r1(0) + a00 * t_r1(1)) * phase_shift;
+    return {a00, a01};
 }

Core/Multilayer/SpecularScalarNCStrategy.h

@@ -34,7 +34,7 @@ private:
     //! reflection coefficients.
     //! Implementation follows A. Gibaud and G. Vignaud, in X-ray and Neutron Reflectivity, edited
     //! by J. Daillant and A. Gibaud, volume 770 of Lecture Notes in Physics (2009)
-    virtual Eigen::Vector2cd transition(complex_t kzi, complex_t kzi1, double sigma,
+    virtual std::pair<complex_t, complex_t> transition(complex_t kzi, complex_t kzi1, double sigma,
                                         double thickness,
                                         const Eigen::Vector2cd& t_r1) const override;
 };

Core/Multilayer/SpecularScalarStrategy.cpp

@@ -94,21 +94,18 @@ void SpecularScalarStrategy::calculateUpFromLayer(std::vector<ScalarRTCoefficien
         if (const auto roughness = GetBottomRoughness(slices, i))
             sigma = roughness->getSigma();
 
-        coeff[i].t_r = transition(kz[i], kz[i + 1], sigma, slices[i].thickness(), coeff[i + 1].t_r);
+        auto mpmm = transition(kz[i], kz[i + 1], sigma, slices[i].thickness(), coeff[i + 1].t_r);
+        auto mp = std::get<0>(mpmm);
+        auto mm = std::get<1>(mpmm);
 
-        if (std::isinf(std::norm(coeff[i].t_r(0))) || std::isnan(std::norm(coeff[i].t_r(0)))) {
-            coeff[i].t_r(0) = 1.0;
-            coeff[i].t_r(1) = 0.0;
+        auto delta = exp_I(kz[i] * slices[i].thickness());
 
-            setZeroBelow(coeff, i);
-        }
+        complex_t S = mp + mm * coeff[i+1].t_r(1);
+        S = 1. / S * delta;
+        factors[i] = S;
 
-        // normalize the t-coefficient in each step of the computation
-        // this avoids an overflow in the backwards propagation
-        // the used factors are stored in order to correct the amplitudes
-        // in forward direction at the end of the computation
-        factors[i] = coeff[i].t_r(0);
-        coeff[i].t_r = coeff[i].t_r / coeff[i].t_r(0);
+        coeff[i].t_r(0) = 1;
+        coeff[i].t_r(1) = delta * (mm + mp * coeff[i+1].t_r(1)) * S;
     }
 
     // now correct all amplitudes by dividing the with the remaining factors in forward direction
@@ -121,7 +118,8 @@ void SpecularScalarStrategy::calculateUpFromLayer(std::vector<ScalarRTCoefficien
             setZeroBelow(coeff, j - 1);
             break;
         }
-        coeff[j].t_r = coeff[j].t_r / dumpingFactor;
+        coeff[j].t_r(0) = dumpingFactor;
+        coeff[j].t_r(1) *= dumpingFactor;
     }
 }
 

Core/Multilayer/SpecularScalarStrategy.h

@@ -47,7 +47,7 @@ public:
                                                 const std::vector<complex_t>& kz) const override;
 
 private:
-    virtual Eigen::Vector2cd transition(complex_t kzi, complex_t kzi1, double sigma,
+    virtual std::pair<complex_t, complex_t> transition(complex_t kzi, complex_t kzi1, double sigma,
                                         double thickness, const Eigen::Vector2cd& t_r1) const = 0;
 
     std::vector<ScalarRTCoefficients> computeTR(const std::vector<Slice>& slices,

Core/Multilayer/SpecularScalarTanhStrategy.cpp

@@ -22,7 +22,8 @@ namespace
 const double pi2_15 = std::pow(M_PI_2, 1.5);
 }
 
-Eigen::Vector2cd SpecularScalarTanhStrategy::transition(complex_t kzi, complex_t kzi1, double sigma,
+std::pair<complex_t, complex_t>
+SpecularScalarTanhStrategy::transition(complex_t kzi, complex_t kzi1, double sigma,
                                                         double thickness,
                                                         const Eigen::Vector2cd& t_r1) const
 {
@@ -39,8 +40,8 @@ Eigen::Vector2cd SpecularScalarTanhStrategy::transition(complex_t kzi, complex_t
     const complex_t a00 = 0.5 * (inv_roughness + kz_ratio);
     const complex_t a01 = 0.5 * (inv_roughness - kz_ratio);
 
-    Eigen::Vector2cd result;
-    result << (a00 * t_r1(0) + a01 * t_r1(1)) / phase_shift,
-        (a01 * t_r1(0) + a00 * t_r1(1)) * phase_shift;
-    return result;
+//    Eigen::Vector2cd result;
+//    result << (a00 * t_r1(0) + a01 * t_r1(1)) / phase_shift,
+//        (a01 * t_r1(0) + a00 * t_r1(1)) * phase_shift;
+    return {a00, a01};
 }

Core/Multilayer/SpecularScalarTanhStrategy.h

@@ -30,7 +30,7 @@ class BA_CORE_API_ SpecularScalarTanhStrategy : public SpecularScalarStrategy
 {
 private:
     //! Roughness is modelled by tanh profile [see e.g. Phys. Rev. B, vol. 47 (8), p. 4385 (1993)].
-    virtual Eigen::Vector2cd transition(complex_t kzi, complex_t kzi1, double sigma,
+    virtual std::pair<complex_t, complex_t> transition(complex_t kzi, complex_t kzi1, double sigma,
                                         double thickness,
                                         const Eigen::Vector2cd& t_r1) const override;
 };