Fix complex overflows in SpecularMatrix

Core/Algorithms/ScalarSpecularInfoMap.cpp

@@ -49,9 +49,5 @@ const ScalarRTCoefficients *ScalarSpecularInfoMap::getCoefficients(kvector_t kve
     SpecularMatrix::MultiLayerCoeff_t coeffs;
     SpecularMatrix::execute(*mp_multilayer, kvec, coeffs);
     auto layer_coeffs = coeffs[m_layer];
-    if (std::isnan(layer_coeffs.getScalarT().real())) {
-        throw RuntimeErrorException("ScalarSpecularInfoMap::getCoefficients() -> "
-                                    "Error! Amplitude is NaN");
-    }
     return new ScalarRTCoefficients(layer_coeffs);
 }

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
@@ -90,7 +92,7 @@ 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::isnan(std::abs(coeff[i].t_r(0)))) {
+        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);
@@ -98,16 +100,24 @@ void SpecularMatrix::execute(const MultiLayer& sample, const kvector_t k,
     }
     // 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 SpecularMatrix::setZeroBelow(MultiLayerCoeff_t& coeff, size_t current_layer)
+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/Algorithms/SpecularMatrix.h

@@ -37,8 +37,6 @@ public:
     //! Computes refraction angles and transmission/reflection coefficients
     //! for given coherent wave propagation in a multilayer.
     static void execute(const MultiLayer& sample, const kvector_t k, MultiLayerCoeff_t& coeff);
-private:
-    void setZeroBelow(MultiLayerCoeff_t& coeff, size_t current_layer);
 };
 
 #endif /* SPECULARMATRIX_H_ */