remove approximation shortcut for magnetic_amplitude when rhoM < MINI…

…MAL_RHO_M

Users noticed an incorrect calculation of the spin asymmetry when the applied field was less than
the equivalent MINIMAL_RHO_M (when H < 0.00425 T).  The reason for the incorrectness is a shortcut that we
were taking to speed up calculations when rhoM in the fronting or backing medium is small, but
a) this shortcut only speeds up the calculations by a factor of 2
b) it seems to be noticeably wrong whenever it is applied, getting worse the higher the field is
so in the interest of correctness, we can remove the shortcut

refl1d/lib/c/magnetic.cc

@@ -12,7 +12,6 @@
 #define M_PI 3.141592653589793
 #endif
 
-#define MINIMAL_RHO_M 1e-2  // in units of 1e-6/A^2
 const double EPS = std::numeric_limits<double>::epsilon();
 const double B2SLD = 2.31604654;  // Scattering factor for B field 1e-6
 
@@ -437,35 +436,24 @@ magnetic_amplitude(const int layers,
 {
   Cplx dummy1,dummy2;
   int ip;
-  if (fabs(rhoM[0]) <= MINIMAL_RHO_M && fabs(rhoM[layers-1]) <= MINIMAL_RHO_M) {
-    ip = 1; // calculations for I+ and I- are the same in the fronting and backing.
-    #ifdef _OPENMP
-    #pragma omp parallel for
-    #endif
-    for (int i=0; i < points; i++) {
-      const int offset = layers*(rho_index != NULL?rho_index[i]:0);
-      Cr4xa(layers,d,sigma,ip,rho+offset,irho+offset,rhoM,u1,u3,
-            KZ[i],Ra[i],Rb[i],Rc[i],Rd[i]);
-    }
-  } else {
-    ip = 1; // plus polarization
-    #ifdef _OPENMP
-    #pragma omp parallel for
-    #endif
-    for (int i=0; i < points; i++) {
-      const int offset = layers*(rho_index != NULL?rho_index[i]:0);
-      Cr4xa(layers,d,sigma,ip,rho+offset,irho+offset,rhoM,u1,u3,
-            KZ[i],Ra[i],Rb[i],dummy1,dummy2);
-    }
-    ip = -1; // minus polarization
-    #ifdef _OPENMP
-    #pragma omp parallel for
-    #endif
-    for (int i=0; i < points; i++) {
-      const int offset = layers*(rho_index != NULL?rho_index[i]:0);
-      Cr4xa(layers,d,sigma,ip,rho+offset,irho+offset,rhoM,u1,u3,
-            KZ[i],dummy1,dummy2,Rc[i],Rd[i]);
-    }
+
+  ip = 1; // plus polarization
+  #ifdef _OPENMP
+  #pragma omp parallel for
+  #endif
+  for (int i=0; i < points; i++) {
+    const int offset = layers*(rho_index != NULL?rho_index[i]:0);
+    Cr4xa(layers,d,sigma,ip,rho+offset,irho+offset,rhoM,u1,u3,
+          KZ[i],Ra[i],Rb[i],dummy1,dummy2);
+  }
+  ip = -1; // minus polarization
+  #ifdef _OPENMP
+  #pragma omp parallel for
+  #endif
+  for (int i=0; i < points; i++) {
+    const int offset = layers*(rho_index != NULL?rho_index[i]:0);
+    Cr4xa(layers,d,sigma,ip,rho+offset,irho+offset,rhoM,u1,u3,
+          KZ[i],dummy1,dummy2,Rc[i],Rd[i]);
   }
 }
 

refl1d/lib/python/magnetic.py

@@ -5,7 +5,6 @@
 M_PI = pi
 PI4 = 4.0e-6 * pi
 B2SLD = 2.31604654  # Scattering factor for B field 1e-6
-MINIMAL_RHO_M = 1e-2  # in units of 1e-6/A^2
 
 prange = range
 
@@ -306,20 +305,16 @@ def magnetic_amplitude(d, sigma, rho, irho, rhoM, u1, u3, KZ, rho_index, Ra, Rb,
     # assert rho_index is None
     layers = len(d)
     points = len(KZ)
-    if fabs(rhoM[0]) <= MINIMAL_RHO_M and fabs(rhoM[layers - 1]) <= MINIMAL_RHO_M:
-        # calculations for I+ and I- are the same in the fronting and backing.
-        for i in prange(points):
-            Cr4xa(layers, d, sigma, 1.0, rho, irho, rhoM, u1, u3, KZ[i], i, Ra, Rb, Rc, Rd)
-    else:
-        # plus polarization must be before minus polarization because it
-        # fills in all R++, R+-, R-+, R--, but minus polarization only fills
-        # in R-+, R--.
-        for i in prange(points):
-            Cr4xa(layers, d, sigma, 1.0, rho, irho, rhoM, u1, u3, KZ[i], i, Ra, Rb, Rc, Rd)
-
-        # minus polarization
-        for i in prange(points):
-            Cr4xa(layers, d, sigma, -1.0, rho, irho, rhoM, u1, u3, KZ[i], i, Ra, Rb, Rc, Rd)
+
+    # plus polarization must be before minus polarization because it
+    # fills in all R++, R+-, R-+, R--, but minus polarization only fills
+    # in R-+, R--.
+    for i in prange(points):
+        Cr4xa(layers, d, sigma, 1.0, rho, irho, rhoM, u1, u3, KZ[i], i, Ra, Rb, Rc, Rd)
+
+    # minus polarization
+    for i in prange(points):
+        Cr4xa(layers, d, sigma, -1.0, rho, irho, rhoM, u1, u3, KZ[i], i, Ra, Rb, Rc, Rd)
 
 
 BASE_GUIDE_ANGLE = 270.0