diff --git a/Core/Algorithms/SpecularMagnetic.cpp b/Core/Algorithms/SpecularMagnetic.cpp
index f71a7f80c663c8b9579a87caffdd36e19e7b3858..ff5689d091dfecec098a9ddd1819e8b6a887d1d8 100644
--- a/Core/Algorithms/SpecularMagnetic.cpp
+++ b/Core/Algorithms/SpecularMagnetic.cpp
@@ -57,6 +57,7 @@ void SpecularMagnetic::calculateEigenvalues(const MultiLayer& sample,
     }
 }
 
+// todo: avoid overflows (see SpecularMatrix.cpp)
 void SpecularMagnetic::calculateTransferAndBoundary(const MultiLayer& sample,
         const kvector_t k, MultiLayerCoeff_t& coeff)
 {
diff --git a/Core/Algorithms/SpecularMatrix.cpp b/Core/Algorithms/SpecularMatrix.cpp
index d6b33d70320ac857d71554e1a1d57e39ad484e01..cd93e873fb8ef0d0b25e7f77a6009a4087bec22a 100644
--- a/Core/Algorithms/SpecularMatrix.cpp
+++ b/Core/Algorithms/SpecularMatrix.cpp
@@ -25,6 +25,8 @@ namespace {
 }
 
 void setZeroBelow(SpecularMatrix::MultiLayerCoeff_t& coeff, size_t current_layer);
+bool calculateUpFromLayer(SpecularMatrix::MultiLayerCoeff_t& coeff, const MultiLayer& sample,
+                          const kvector_t k, size_t layer_index);
 
 // Computes refraction angles and transmission/reflection coefficients
 // for given coherent wave propagation in a multilayer.
@@ -66,7 +68,42 @@ void SpecularMatrix::execute(const MultiLayer& sample, const kvector_t k,
     }
     // Calculate transmission/refraction coefficients t_r for each layer,
     // from bottom to top.
-    for (int i=N-2; i>=0; --i) {
+    size_t start_layer_index = N-2;
+    while (!calculateUpFromLayer(coeff, sample, k, start_layer_index)) {
+        setZeroBelow(coeff, start_layer_index);
+        coeff[start_layer_index].t_r(0) = 1.0;
+        coeff[start_layer_index].t_r(1) = 0.0;
+        --start_layer_index;
+    }
+
+    // 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();
+    }
+}
+
+bool calculateUpFromLayer(SpecularMatrix::MultiLayerCoeff_t& coeff, const MultiLayer& sample,
+                          const kvector_t k, size_t layer_index)
+{
+    for (int i=layer_index; i>=0; --i) {
         complex_t roughness_factor = 1;
         if (sample.getLayerInterface(i)->getRoughness()) {
             double sigma = sample.getLayerBottomInterface(i)->getRoughness()->getSigma();
@@ -93,31 +130,8 @@ void SpecularMatrix::execute(const MultiLayer& sample, const kvector_t k,
                     (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);
+            return false;
         }
     }
-    // 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();
-    }
+    return true;
 }