diff --git a/Resample/Specular/ComputeFluxMagnetic.cpp b/Resample/Specular/ComputeFluxMagnetic.cpp
index c213c08911d83abb473a00f106927a9f0daaece1..6aeb955e416720c2a8c1dbc0e0d56d44fd16dedc 100644
--- a/Resample/Specular/ComputeFluxMagnetic.cpp
+++ b/Resample/Specular/ComputeFluxMagnetic.cpp
@@ -21,6 +21,7 @@
 #include "Resample/Specular/TransitionMagneticNevot.h"
 #include "Resample/Specular/TransitionMagneticTanh.h"
 #include "Sample/Interface/LayerRoughness.h"
+#include <algorithm>
 
 namespace {
 
@@ -79,9 +80,8 @@ void calculateUpwards(std::vector<MatrixFlux>& coeff, const SliceStack& slices,
             coeff[i - 1].computeDeltaMatrix(slices[i - 1].thicknessOr0());
 
         // compute the rotation matrix
-        SpinMatrix S, Si;
-        Si = mp + mm * coeff[i].m_R;
-        S << Si(1, 1), -Si(0, 1), -Si(1, 0), Si(0, 0);
+        SpinMatrix Si = mp + mm * coeff[i].m_R;
+        SpinMatrix S =MakeSpinMatrix( 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;
 
@@ -154,10 +154,10 @@ std::vector<MatrixFlux> computeFlux(const SliceStack& slices, const std::vector<
 
     if (kzs[0] == 0.) {
         result[0].m_T = SpinMatrix::Identity();
-        result[0].m_R = -SpinMatrix::Identity();
+        result[0].m_R = MakeSpinMatrix(-1, 0, 0, -1);
         for (size_t i = 1; i < N; ++i) {
-            result[i].m_T.setZero();
-            result[i].m_R.setZero();
+            result[i].m_T = SpinMatrix::Zero();
+            result[i].m_R = SpinMatrix::Zero();
         }
         return result;
     }
@@ -230,9 +230,8 @@ SpinMatrix Compute::SpecularMagnetic::topLayerR(const SliceStack& slices,
         const SpinMatrix delta = c_i.computeDeltaMatrix(slices[i - 1].thicknessOr0());
 
         // compute the rotation matrix
-        SpinMatrix S, Si;
-        Si = mp + mm * c_i1.m_R;
-        S << Si(1, 1), -Si(0, 1), -Si(1, 0), Si(0, 0);
+        SpinMatrix Si = mp + mm * c_i1.m_R;
+        SpinMatrix S =MakeSpinMatrix( 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;
         S /= norm;
diff --git a/Resample/Specular/ComputeFluxScalar.cpp b/Resample/Specular/ComputeFluxScalar.cpp
index 0492be025b59abee9ed2d69a1011cf79ca8d010e..8b45ce0df9f403917e6f69d467efa10b4ad8f019 100644
--- a/Resample/Specular/ComputeFluxScalar.cpp
+++ b/Resample/Specular/ComputeFluxScalar.cpp
@@ -80,7 +80,7 @@ std::vector<Spinor> computeTR(const SliceStack& slices, const std::vector<comple
     if (kz[X[0]] == 0.0) { // If kz in layer 0 is zero, R0 = -T0 and all others equal to 0
         TR[X[0]] = {1.0, -1.0};
         for (size_t i = 1; i < N; ++i)
-            TR[X[i]].setZero();
+            TR[X[i]] = Spinor(0, 0);
         return TR;
     }