Merge branch 'nanbug' into develop

Core/Algorithms/ScalarSpecularInfoMap.cpp

@@ -28,26 +28,26 @@ ScalarSpecularInfoMap *ScalarSpecularInfoMap::clone() const
     return new ScalarSpecularInfoMap(mp_multilayer, m_layer);
 }
 
-//! \todo Can we avoid the code duplication in the two following functions ?
-
 const ScalarRTCoefficients *ScalarSpecularInfoMap::getOutCoefficients(
         double alpha_f, double phi_f, double wavelength) const
 {
-    (void)phi_f;
-    SpecularMatrix::MultiLayerCoeff_t coeffs;
-    // phi has no effect on R,T, so just pass zero:
+    (void)phi_f; // phi has no effect on R,T, so just pass zero:
     kvector_t kvec = Geometry::vecOfLambdaAlphaPhi(wavelength, -alpha_f, 0.0);
-    SpecularMatrix::execute(*mp_multilayer, kvec, coeffs);
-    return new ScalarRTCoefficients(coeffs[m_layer]);
+    return getCoefficients(kvec);
 }
 
 const ScalarRTCoefficients *ScalarSpecularInfoMap::getInCoefficients(
         double alpha_i, double phi_i, double wavelength) const
 {
-    (void)phi_i;
-    SpecularMatrix::MultiLayerCoeff_t coeffs;
-    // phi has no effect on R,T, so just pass zero:
+    (void)phi_i; // phi has no effect on R,T, so just pass zero:
     kvector_t kvec = Geometry::vecOfLambdaAlphaPhi(wavelength, alpha_i, 0.0);
+    return getCoefficients(kvec);
+}
+
+const ScalarRTCoefficients *ScalarSpecularInfoMap::getCoefficients(kvector_t kvec) const
+{
+    SpecularMatrix::MultiLayerCoeff_t coeffs;
     SpecularMatrix::execute(*mp_multilayer, kvec, coeffs);
-    return new ScalarRTCoefficients(coeffs[m_layer]);
+    auto layer_coeffs = coeffs[m_layer];
+    return new ScalarRTCoefficients(layer_coeffs);
 }

Core/Algorithms/ScalarSpecularInfoMap.h

@@ -45,6 +45,7 @@ public:
 private:
     const MultiLayer *mp_multilayer;
     const int m_layer;
+    const ScalarRTCoefficients *getCoefficients(kvector_t kvec) const;
 };
 
 #endif /* SCALARSPECULARINFOMAP_H_ */

Core/Algorithms/SpecularMatrix.cpp

@@ -24,6 +24,8 @@ namespace {
     const complex_t imag_unit = complex_t(0.0, 1.0);
 }
 
+void setZeroBelow(SpecularMatrix::MultiLayerCoeff_t& coeff, size_t current_layer);
+
 // Computes refraction angles and transmission/reflection coefficients
 // for given coherent wave propagation in a multilayer.
 // k : length: wavenumber in vacuum, direction: defined in layer 0
@@ -59,9 +61,7 @@ void SpecularMatrix::execute(const MultiLayer& sample, const kvector_t k,
     if (coeff[0].lambda == 0.0) {
         coeff[0].t_r(0) = 1.0;
         coeff[0].t_r(1) = -1.0;
-        for (size_t i=1; i<N; ++i) {
-            coeff[i].t_r.setZero();
-        }
+        setZeroBelow(coeff, 0);
         return;
     }
     // Calculate transmission/refraction coefficients t_r for each layer,
@@ -92,10 +92,32 @@ void SpecularMatrix::execute(const MultiLayer& sample, const kvector_t k,
                     (lambda_rough-lambda_below)*coeff[i+1].t_r(0) +
                     (lambda_rough+lambda_below)*coeff[i+1].t_r(1) )/2.0/lambda *
                     std::exp( ikd*lambda);
+        if (std::isinf(std::abs(coeff[i].getScalarT()))) {
+            coeff[i].t_r(0) = 1.0;
+            coeff[i].t_r(1) = 0.0;
+            setZeroBelow(coeff, i);
+        }
     }
     // Normalize to incoming downward travelling amplitude = 1.
     complex_t T0 = coeff[0].getScalarT();
+    // This trick is used to avoid overflows in the intermediate steps of
+    // complex division:
+    double tmax = std::max(std::abs(T0.real()), std::abs(T0.imag()));
+    for (size_t i=0; i<N; ++i) {
+        coeff[i].t_r(0) /= tmax;
+        coeff[i].t_r(1) /= tmax;
+    }
+    // Now the real normalization to T_0 proceeds
+    T0 = coeff[0].getScalarT();
     for (size_t i=0; i<N; ++i) {
         coeff[i].t_r = coeff[i].t_r/T0;
     }
 }
+
+void setZeroBelow(SpecularMatrix::MultiLayerCoeff_t& coeff, size_t current_layer)
+{
+    size_t N = coeff.size();
+    for (size_t i=current_layer+1; i<N; ++i) {
+        coeff[i].t_r.setZero();
+    }
+}

Core/FormFactors/FormFactorDWBA.cpp

@@ -47,7 +47,8 @@ void FormFactorDWBA::setSpecularInfo(const ILayerRTCoefficients *p_in_coeffs,
 complex_t FormFactorDWBA::evaluate(const WavevectorInfo& wavevectors) const
 {
     calculateTerms(wavevectors);
-    return m_term_S + m_term_RS + m_term_SR + m_term_RSR;
+    complex_t sum_terms = m_term_S + m_term_RS + m_term_SR + m_term_RSR;
+    return sum_terms;
 }
 
 void FormFactorDWBA::calculateTerms(const WavevectorInfo& wavevectors) const
@@ -71,14 +72,16 @@ void FormFactorDWBA::calculateTerms(const WavevectorInfo& wavevectors) const
     WavevectorInfo k_TR(k_i_T, k_f_R, wavelength);
     WavevectorInfo k_RR(k_i_R, k_f_R, wavelength);
 
+    // Get the four R,T coefficients
+    complex_t T_in = mp_in_coeffs->getScalarT();
+    complex_t R_in = mp_in_coeffs->getScalarR();
+    complex_t T_out = mp_out_coeffs->getScalarT();
+    complex_t R_out = mp_out_coeffs->getScalarR();
+
     // The four different scattering contributions; S stands for scattering
     // off the particle, R for reflection off the layer interface
-    m_term_S = mp_in_coeffs->getScalarT()*mp_form_factor->evaluate(k_TT)
-            * mp_out_coeffs->getScalarT();
-    m_term_RS = mp_in_coeffs->getScalarR()*mp_form_factor->evaluate(k_RT)
-            * mp_out_coeffs->getScalarT();
-    m_term_SR = mp_in_coeffs->getScalarT()*mp_form_factor->evaluate(k_TR)
-            * mp_out_coeffs->getScalarR();
-    m_term_RSR = mp_in_coeffs->getScalarR()*mp_form_factor->evaluate(k_RR)
-            * mp_out_coeffs->getScalarR();
+    m_term_S = T_in * mp_form_factor->evaluate(k_TT) * T_out;
+    m_term_RS = R_in * mp_form_factor->evaluate(k_RT) * T_out;
+    m_term_SR = T_in * mp_form_factor->evaluate(k_TR) * R_out;
+    m_term_RSR = R_in * mp_form_factor->evaluate(k_RR) * R_out;
 }