Changeset 885753a in sasmodels

Timestamp:

Apr 6, 2018 10:35:22 AM (7 years ago)

Author:

Paul Kienzle <pkienzle@…>

Branches:

master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

dc6f601

Parents:

5770493

Message:

correct magnetic calculations; fractional polarization no longer matches 3.1.2

Location:

sasmodels

Files:

• details.py (modified) (3 diffs)
• kernel_iq.c (modified) (5 diffs)

sasmodels/details.py

@@ -243 +243 @@
     offset = np.cumsum(np.hstack((0, length)))
     call_details = make_details(kernel.info, length, offset[:-1], offset[-1])
-    # Pad value array to a 32 value boundaryd
+    # Pad value array to a 32 value boundary
     data_len = nvalues + 2*sum(len(v) for v in dispersity)
     extra = (32 - data_len%32)%32
@@ -250 +250 @@
     is_magnetic = convert_magnetism(kernel.info.parameters, data)
     #call_details.show()
+    #print("data", data)
     return call_details, data, is_magnetic
 
@@ -296 +297 @@
     mag = values[parameters.nvalues-3*parameters.nmagnetic:parameters.nvalues]
     mag = mag.reshape(-1, 3)
-    scale = mag[:, 0]
-    if np.any(scale):
+    if np.any(mag[:, 0] != 0.0):
+        M0 = mag[:, 0].copy()
         theta, phi = radians(mag[:, 1]), radians(mag[:, 2])
-        cos_theta = cos(theta)
-        mag[:, 0] = scale*cos_theta*cos(phi)  # mx
-        mag[:, 1] = scale*sin(theta)  # my
-        mag[:, 2] = -scale*cos_theta*sin(phi)  # mz
+        mag[:, 0] = +M0*cos(theta)*cos(phi)  # mx
+        mag[:, 1] = +M0*sin(theta) # my
+        mag[:, 2] = -M0*cos(theta)*sin(phi)  # mz
         return True
     else:

sasmodels/kernel_iq.c

@@ -79 +79 @@
 //     du * (m_sigma_y + 1j*m_sigma_z);
 // weights for spin crosssections: dd du real, ud real, uu, du imag, ud imag
-static void set_spin_weights(double in_spin, double out_spin, double spins[4])
+static void set_spin_weights(double in_spin, double out_spin, double weight[6])
 {
   in_spin = clip(in_spin, 0.0, 1.0);
   out_spin = clip(out_spin, 0.0, 1.0);
-  spins[0] = sqrt(sqrt((1.0-in_spin) * (1.0-out_spin))); // dd
-  spins[1] = sqrt(sqrt((1.0-in_spin) * out_spin));       // du real
-  spins[2] = sqrt(sqrt(in_spin * (1.0-out_spin)));       // ud real
-  spins[3] = sqrt(sqrt(in_spin * out_spin));             // uu
-  spins[4] = spins[1]; // du imag
-  spins[5] = spins[2]; // ud imag
+  // Note: sasview 3.1 scaled all slds by sqrt(weight) and assumed that
+  //     w*I(q, rho1, rho2, ...) = I(q, sqrt(w)*rho1, sqrt(w)*rho2, ...)
+  // which is likely to be the case for simple models.
+  weight[0] = ((1.0-in_spin) * (1.0-out_spin)); // dd
+  weight[1] = ((1.0-in_spin) * out_spin);       // du.real
+  weight[2] = (in_spin * (1.0-out_spin));       // ud.real
+  weight[3] = (in_spin * out_spin);             // uu
+  weight[4] = weight[1]; // du.imag
+  weight[5] = weight[2]; // ud.imag
 }
 
 // Compute the magnetic sld
 static double mag_sld(
-  const unsigned int xs, // 0=dd, 1=du real, 2=ud real, 3=uu, 4=du imag, 5=up imag
+  const unsigned int xs, // 0=dd, 1=du.real, 2=ud.real, 3=uu, 4=du.imag, 5=ud.imag
   const double qx, const double qy,
   const double px, const double py,
@@ -106 +109 @@
       case 0: // uu => sld - D M_perpx
           return sld - px*perp;
-      case 1: // ud real => -D M_perpy
+      case 1: // ud.real => -D M_perpy
           return py*perp;
-      case 2: // du real => -D M_perpy
+      case 2: // du.real => -D M_perpy
           return py*perp;
-      case 3: // dd real => sld + D M_perpx
+      case 3: // dd => sld + D M_perpx
           return sld + px*perp;
     }
   } else {
     if (xs== 4) {
-      return -mz;  // ud imag => -D M_perpz
+      return -mz;  // du.imag => +D M_perpz
     } else { // index == 5
-      return mz;   // du imag => D M_perpz
+      return +mz;  // ud.imag => -D M_perpz
     }
   }
@@ -312 +315 @@
   //     up_angle = values[NUM_PARS+4];
   // TODO: could precompute more magnetism parameters before calling the kernel.
-  double spins[8];  // uu, ud real, du real, dd, ud imag, du imag, fill, fill
+  double xs_weights[8];  // uu, ud real, du real, dd, ud imag, du imag, fill, fill
   double cos_mspin, sin_mspin;
-  set_spin_weights(values[NUM_PARS+2], values[NUM_PARS+3], spins);
+  set_spin_weights(values[NUM_PARS+2], values[NUM_PARS+3], xs_weights);
   SINCOS(-values[NUM_PARS+4]*M_PI_180, sin_mspin, cos_mspin);
 #endif // MAGNETIC
@@ -661 +664 @@
             // loop over uu, ud real, du real, dd, ud imag, du imag
             for (unsigned int xs=0; xs<6; xs++) {
-              const double xs_weight = spins[xs];
+              const double xs_weight = xs_weights[xs];
               if (xs_weight > 1.e-8) {
                 // Since the cross section weight is significant, set the slds
@@ -674 +677 @@
                   local_values.vector[sld_index] =
                     mag_sld(xs, qx, qy, px, py, values[sld_index+2], mx, my, mz);
+//if (q_index==0) printf("%d: (qx,qy)=(%g,%g) xs=%d sld%d=%g p=(%g,%g) m=(%g,%g,%g)\n",
+//q_index, qx, qy, xs, sk, local_values.vector[sld_index], px, py, mx, my, mz);
                 }
                 scattering += xs_weight * CALL_KERNEL();