diff --git a/Resample/Specular/ComputeFluxMagnetic.cpp b/Resample/Specular/ComputeFluxMagnetic.cpp
index 020339c653ebd09bbee26f67e40660aa5028c0a3..2e76fc737a8a365c05e36c419bc867eedd42e17d 100644
--- a/Resample/Specular/ComputeFluxMagnetic.cpp
+++ b/Resample/Specular/ComputeFluxMagnetic.cpp
@@ -67,34 +67,34 @@ void calculateUpwards(std::vector<MatrixFlux>& coeff, const SliceStack& slices,
     coeff.back().m_T = Eigen::Matrix2cd::Identity();
     coeff.back().m_R = Eigen::Matrix2cd::Zero();
 
-    for (signed long i = N - 2; i >= 0; --i) {
+    for (size_t i = N - 1; i > 0; --i) {
         double sigma = 0.;
         if (const auto roughness = slices.bottomRoughness(i))
             sigma = roughness->getSigma();
 
         // compute the 2x2 submatrices in the write-up denoted as P+, P- and delta
-        const auto [mp, mm] = computeBackwardsSubmatrices(coeff[i], coeff[i + 1], sigma, r_model);
+        const auto [mp, mm] = computeBackwardsSubmatrices(coeff[i-1], coeff[i], sigma, r_model);
 
-        const Eigen::Matrix2cd delta = coeff[i].computeDeltaMatrix(slices[i].thicknessOr0());
+        const Eigen::Matrix2cd delta = coeff[i-1].computeDeltaMatrix(slices[i-1].thicknessOr0());
 
         // compute the rotation matrix
         Eigen::Matrix2cd S, Si;
-        Si = mp + mm * coeff[i + 1].m_R;
+        Si = mp + mm * coeff[i].m_R;
         S << Si(1, 1), -Si(0, 1), -Si(1, 0), Si(0, 0);
         const complex_t norm = S(0, 0) * S(1, 1) - S(0, 1) * S(1, 0);
         S = S * delta;
 
         // store the rotation matrix and normalization constant in order to rotate
         // the coefficients for all lower slices at the end of the computation
-        SMatrices[i] = S;
-        Normalization[i] = norm;
+        SMatrices[i-1] = S;
+        Normalization[i-1] = norm;
 
         // compute the reflection matrix and
         // rotate the polarization such that we have pure incoming states (T = I)
         S /= norm;
 
         // T is always equal to the identity at this point, no need to store
-        coeff[i].m_R = delta * (mm + mp * coeff[i + 1].m_R) * S;
+        coeff[i-1].m_R = delta * (mm + mp * coeff[i].m_R) * S;
     }
 
     // now correct all amplitudes in forward direction by dividing with the remaining
@@ -216,17 +216,17 @@ Eigen::Matrix2cd Compute::SpecularMagnetic::topLayerR(const SliceStack& slices,
     // bottom boundary condition
     c_i1.m_R = Eigen::Matrix2cd::Zero();
 
-    for (signed long i = N - 2; i >= 0; --i) {
-        auto c_i = createCoeff(i);
+    for (size_t i = N - 1; i > 0; --i) {
+        auto c_i = createCoeff(i-1);
 
         double sigma = 0.;
-        if (const auto roughness = slices.bottomRoughness(i))
+        if (const auto roughness = slices.bottomRoughness(i-1))
             sigma = roughness->getSigma();
 
         // compute the 2x2 submatrices in the write-up denoted as P+, P- and delta
         const auto [mp, mm] = computeBackwardsSubmatrices(c_i, c_i1, sigma, r_model);
 
-        const Eigen::Matrix2cd delta = c_i.computeDeltaMatrix(slices[i].thicknessOr0());
+        const Eigen::Matrix2cd delta = c_i.computeDeltaMatrix(slices[i-1].thicknessOr0());
 
         // compute the rotation matrix
         Eigen::Matrix2cd S, Si;
diff --git a/Resample/Specular/ComputeFluxScalar.cpp b/Resample/Specular/ComputeFluxScalar.cpp
index 8726a0dd44575750f5df6eb45ca61c551e09b5a6..599ce48eec336464c030a3fc2ad9916ce420deed 100644
--- a/Resample/Specular/ComputeFluxScalar.cpp
+++ b/Resample/Specular/ComputeFluxScalar.cpp
@@ -88,9 +88,9 @@ std::vector<Eigen::Vector2cd> computeTR(const SliceStack& slices, const std::vec
     // Calculate transmission/refraction coefficients t_r for each layer, from bottom to top.
     TR[X[N - 1]] = {1.0, 0.0};
     std::vector<complex_t> factors(N - 1);
-    for (int i = N - 2; i >= 0; i--) {
-        size_t jthis = X[i];
-        size_t jlast = X[i + 1];
+    for (size_t i = N - 1; i > 0; i--) {
+        size_t jthis = X[i - 1];
+        size_t jlast = X[i];
         const auto roughness = slices.bottomRoughness(jthis); // TODO verify
         const double sigma = roughness ? roughness->getSigma() : 0.;
 
@@ -99,7 +99,7 @@ std::vector<Eigen::Vector2cd> computeTR(const SliceStack& slices, const std::vec
         const complex_t delta = exp_I(kz[jthis] * slices[jthis].thicknessOr0());
 
         complex_t S = delta / (mp + mm * TR[jlast](1));
-        factors[i] = S;
+        factors[i-1] = S;
         TR[jthis](1) = delta * (mm + mp * TR[jlast](1)) * S;
     }
 
@@ -151,14 +151,14 @@ complex_t Compute::SpecularScalar::topLayerR(const SliceStack& slices,
 
     complex_t R_i1 = 0.;
 
-    for (int i = N - 2; i >= 0; i--) {
+    for (size_t i = N - 1; i > 0; i--) {
         double sigma = 0.0;
-        if (const auto roughness = slices.bottomRoughness(i))
+        if (const auto roughness = slices.bottomRoughness(i-1))
             sigma = roughness->getSigma();
 
-        const auto [mp, mm] = transition(kz[i], kz[i + 1], sigma, r_model);
+        const auto [mp, mm] = transition(kz[i-1], kz[i], sigma, r_model);
 
-        const complex_t delta = exp_I(kz[i] * slices[i].thicknessOr0());
+        const complex_t delta = exp_I(kz[i-1] * slices[i-1].thicknessOr0());
 
         R_i1 = pow(delta, 2) * (mm + mp * R_i1) / (mp + mm * R_i1);
     }