diff --git a/Resample/Specular/TransitionMagneticNevot.cpp b/Resample/Specular/TransitionMagneticNevot.cpp
index 40f87f2c13de3c74cab8d05937dc947822f05e22..253fa42a90282ea14b704e3dc1cfc6557cc83dae 100644
--- a/Resample/Specular/TransitionMagneticNevot.cpp
+++ b/Resample/Specular/TransitionMagneticNevot.cpp
@@ -32,30 +32,23 @@ std::pair<SpinMatrix, SpinMatrix> Compute::refractionMatrixBlocksNevot(const Mat
         complex_t beta_b = coeff_b.k_eigen_up() - coeff_b.k_eigen_dn();
 
         const complex_t alpha = alpha_b + sign * alpha_a;
-        C3 b_p_vec = beta_b * coeff_b.field() + sign * beta_a * coeff_a.field();
+        C3 b = beta_b * coeff_b.field() + sign * beta_a * coeff_a.field();
 
         auto square = [](auto& v) { return v.x() * v.x() + v.y() * v.y() + v.z() * v.z(); };
-        complex_t beta = std::sqrt(square(b_p_vec));
+        complex_t beta = std::sqrt(square(b));
         if (std::abs(beta) < std::numeric_limits<double>::epsilon() * 10.) {
             const complex_t alpha_pp = -(alpha * alpha) * sigma * sigma / 8.;
             return SpinMatrix(std::exp(alpha_pp), 0, 0, std::exp(alpha_pp));
         }
 
-        b_p_vec /= beta;
+        b /= beta;
 
         const complex_t alpha_pp = -(alpha * alpha + beta * beta) * sigma * sigma / 8.;
         const complex_t beta_pp = -alpha * beta * sigma * sigma / 4.;
+        SpinMatrix Q(b.z() + 1., b.x() - I * b.y(), b.x() + I * b.y(), -1. - b.z());
+        const SpinMatrix M(std::exp(beta_pp), 0, 0, std::exp(-beta_pp));
 
-        const complex_t factor1 = std::sqrt(2. * (1. + b_p_vec.z()));
-        const complex_t factor2 = std::sqrt(2. * (1. - b_p_vec.z()));
-        SpinMatrix Q((b_p_vec.z() + 1.) / factor1, (b_p_vec.z() - 1.) / factor2,
-                     (b_p_vec.x() + I * b_p_vec.y()) / factor1,
-                     (b_p_vec.x() + I * b_p_vec.y()) / factor2);
-
-        const SpinMatrix exp1(std::exp(alpha_pp), 0, 0, std::exp(alpha_pp));
-        const SpinMatrix exp2(std::exp(beta_pp), 0, 0, std::exp(-beta_pp));
-
-        return exp1 * Q * exp2 * Q.adjoint();
+        return std::exp(alpha_pp) * Q * M * Q.adjoint() / (2. * (1. + b.z()));
     };
 
     const auto kk = SpinMatrix(coeff_a.computeInverseKappa() * coeff_b.computeKappa());
diff --git a/Tests/ReferenceData/Suite/MagneticRotationReflectivityMM_NevotCroce.int.gz b/Tests/ReferenceData/Suite/MagneticRotationReflectivityMM_NevotCroce.int.gz
index 4d3236e9370b7936f65e1b77fe7cc90d6224009e..a6dcaed5e0737f9b3f8bb44ded38d0724c04115c 100644
Binary files a/Tests/ReferenceData/Suite/MagneticRotationReflectivityMM_NevotCroce.int.gz and b/Tests/ReferenceData/Suite/MagneticRotationReflectivityMM_NevotCroce.int.gz differ
diff --git a/Tests/ReferenceData/Suite/MagneticRotationReflectivityMP_NevotCroce.int.gz b/Tests/ReferenceData/Suite/MagneticRotationReflectivityMP_NevotCroce.int.gz
index eee57da80d11bb0981ba878efa2f2fec4b8b1468..f73940f5eba338439af61e18e6841d44c459cd11 100644
Binary files a/Tests/ReferenceData/Suite/MagneticRotationReflectivityMP_NevotCroce.int.gz and b/Tests/ReferenceData/Suite/MagneticRotationReflectivityMP_NevotCroce.int.gz differ
diff --git a/Tests/ReferenceData/Suite/MagneticRotationReflectivityPM_NevotCroce.int.gz b/Tests/ReferenceData/Suite/MagneticRotationReflectivityPM_NevotCroce.int.gz
index 37b3e16e07fbeba77d49ac5f419ba658814048e9..06280d7c212cdbef2f32ff85305155ba49641622 100644
Binary files a/Tests/ReferenceData/Suite/MagneticRotationReflectivityPM_NevotCroce.int.gz and b/Tests/ReferenceData/Suite/MagneticRotationReflectivityPM_NevotCroce.int.gz differ
diff --git a/Tests/ReferenceData/Suite/MagneticRotationReflectivityPP_NevotCroce.int.gz b/Tests/ReferenceData/Suite/MagneticRotationReflectivityPP_NevotCroce.int.gz
index cc1ef6e76bb3b39161872469ae103ff6190ffd76..4555407243b6fa5883fa3cef87d664ee77537523 100644
Binary files a/Tests/ReferenceData/Suite/MagneticRotationReflectivityPP_NevotCroce.int.gz and b/Tests/ReferenceData/Suite/MagneticRotationReflectivityPP_NevotCroce.int.gz differ