Changeset a4280bd in sasmodels for sasmodels/kernel_iq.cl

Timestamp:

Jul 25, 2016 11:54:30 PM (8 years ago)

Author:

Paul Kienzle <pkienzle@…>

Branches:

master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

2f5c6d4

Parents:

2c74c11

Message:

restructure magnetic models to use less code

File:

• sasmodels/kernel_iq.cl (modified) (5 diffs)

sasmodels/kernel_iq.cl

@@ -34 +34 @@
 
 
-#ifdef MAGNETIC
+#if defined(MAGNETIC) && NUM_MAGNETIC>0
 
 // Return value restricted between low and high
@@ -48 +48 @@
 //     ud * (m_sigma_y + 1j*m_sigma_z);
 //     du * (m_sigma_y - 1j*m_sigma_z);
-static void spins(double in_spin, double out_spin,
-    double *uu, double *dd, double *ud, double *du)
+static void set_spins(double in_spin, double out_spin, double spins[4])
 {
   in_spin = clip(in_spin, 0.0, 1.0);
   out_spin = clip(out_spin, 0.0, 1.0);
-  *uu = sqrt(sqrt(in_spin * out_spin));
-  *dd = sqrt(sqrt((1.0-in_spin) * (1.0-out_spin)));
-  *ud = sqrt(sqrt(in_spin * (1.0-out_spin)));
-  *du = sqrt(sqrt((1.0-in_spin) * out_spin));
+  spins[0] = sqrt(sqrt((1.0-in_spin) * (1.0-out_spin))); // dd
+  spins[1] = sqrt(sqrt((1.0-in_spin) * out_spin));       // du
+  spins[2] = sqrt(sqrt(in_spin * (1.0-out_spin)));       // ud
+  spins[3] = sqrt(sqrt(in_spin * out_spin));             // uu
+}
+
+static double mag_sld(double qx, double qy, double p,
+                       double mx, double my, double sld)
+{
+    const double perp = qy*mx - qx*my;
+    return sld + perp*p;
 }
 
@@ -84 +90 @@
 
   // Fill in the initial variables
-  for (int i=0; i < NPARS; i++) {
+  for (int i=0; i < NUM_PARS; i++) {
     pvec[i] = values[2+i];
 //if (q_index==0) printf("p%d = %g\n",i, pvec[i]);
   }
 
-#ifdef MAGNETIC
+#if defined(MAGNETIC) && NUM_MAGNETIC>0
   // Location of the sld parameters in the parameter pvec.
   // These parameters are updated with the effective sld due to magnetism.
+  #if NUM_MAGNETIC > 3
   const int32_t slds[] = { MAGNETIC_PARS };
-
-  const double up_frac_i = values[NPARS+2];
-  const double up_frac_f = values[NPARS+3];
-  const double up_angle = values[NPARS+4];
-  #define MX(_k) (values[NPARS+5+3*_k])
-  #define MY(_k) (values[NPARS+6+3*_k])
-  #define MZ(_k) (values[NPARS+7+3*_k])
+  #endif
 
   // TODO: could precompute these outside of the kernel.
   // Interpret polarization cross section.
-  double uu, dd, ud, du;
+  //     up_frac_i = values[NUM_PARS+2];
+  //     up_frac_f = values[NUM_PARS+3];
+  //     up_angle = values[NUM_PARS+4];
+  double spins[4];
   double cos_mspin, sin_mspin;
-  spins(up_frac_i, up_frac_f, &uu, &dd, &ud, &du);
-  SINCOS(-up_angle*M_PI_180, sin_mspin, cos_mspin);
+  set_spins(values[NUM_PARS+2], values[NUM_PARS+3], spins);
+  SINCOS(-values[NUM_PARS+4]*M_PI_180, sin_mspin, cos_mspin);
 #endif // MAGNETIC
 
@@ -222 +226 @@
 #endif
 
-//if (q_index == 0) {printf("step:%d of %d, pars:",step,pd_stop); for (int i=0; i < NPARS; i++) printf("p%d=%g ",i, pvec[i]); printf("\n"); }
+//if (q_index == 0) {printf("step:%d of %d, pars:",step,pd_stop); for (int i=0; i < NUM_PARS; i++) printf("p%d=%g ",i, pvec[i]); printf("\n"); }
 //if (q_index == 0) printf("sphcor: %g\n", spherical_correction);
 
@@ -237 +241 @@
         pd_norm += weight * CALL_VOLUME(local_values);
 
-#ifdef MAGNETIC
-          const double qx = q[2*q_index];
-          const double qy = q[2*q_index+1];
-          const double qsq = qx*qx + qy*qy;
-
-          // Constant across orientation, polydispersity for given qx, qy
-          double px, py, pz;
-          if (qsq > 1.e-16) {
-            px = (qy*cos_mspin + qx*sin_mspin)/qsq;
-            py = (qy*sin_mspin - qx*cos_mspin)/qsq;
-            pz = 1.0;
-          } else {
-            px = py = pz = 0.0;
+#if defined(MAGNETIC) && NUM_MAGNETIC > 0
+        const double qx = q[2*q_index];
+        const double qy = q[2*q_index+1];
+        const double qsq = qx*qx + qy*qy;
+
+        // Constant across orientation, polydispersity for given qx, qy
+        double scattering = 0.0;
+        // TODO: what is the magnetic scattering at q=0
+        if (qsq > 1.e-16) {
+          double p[4];  // spin_i, spin_f
+          p[0] = (qy*cos_mspin + qx*sin_mspin)/qsq;
+          p[3] = -p[0];
+          p[1] = p[2] = (qy*sin_mspin - qx*cos_mspin)/qsq;
+
+          for (int index=0; index<4; index++) {
+            const double xs = spins[index];
+            if (xs > 1.e-8) {
+              const int spin_flip = (index==1) || (index==2);
+              const double pk = p[index];
+              for (int axis=0; axis<=spin_flip; axis++) {
+                #define M1 NUM_PARS+5
+                #define M2 NUM_PARS+8
+                #define M3 NUM_PARS+13
+                #define SLD(_M_offset, _sld_offset) \
+                    pvec[_sld_offset] = xs * (axis \
+                    ? (index==1 ? -values[_M_offset+2] : values[_M_offset+2]) \
+                    : mag_sld(qx, qy, pk, values[_M_offset], values[_M_offset+1], \
+                              (spin_flip ? 0.0 : values[_sld_offset+2])))
+                #if NUM_MAGNETIC==1
+                    SLD(M1, MAGNETIC_PAR1);
+                #elif NUM_MAGNETIC==2
+                    SLD(M1, MAGNETIC_PAR1);
+                    SLD(M2, MAGNETIC_PAR2);
+                #elif NUM_MAGNETIC==3
+                    SLD(M1, MAGNETIC_PAR1);
+                    SLD(M2, MAGNETIC_PAR2);
+                    SLD(M3, MAGNETIC_PAR3);
+                #else
+                for (int sk=0; sk<NUM_MAGNETIC; sk++) {
+                    SLD(M1+3*sk, slds[sk]);
+                }
+                #endif
+                scattering += CALL_IQ(q, q_index, local_values);
+              }
+            }
           }
-
-          double scattering = 0.0;
-          if (uu > 1.e-8) {
-            for (int sk=0; sk<NUM_MAGNETIC; sk++) {
-                const double perp = (qy*MX(sk) - qx*MY(sk));
-                pvec[slds[sk]] = (values[slds[sk]+2] - perp*px)*uu;
-            }
-            scattering += CALL_IQ(q, q_index, local_values);
-          }
-
-          if (dd > 1.e-8){
-            for (int sk=0; sk<NUM_MAGNETIC; sk++) {
-                const double perp = (qy*MX(sk) - qx*MY(sk));
-                pvec[slds[sk]] = (values[slds[sk]+2] + perp*px)*dd;
-            }
-            scattering += CALL_IQ(q, q_index, local_values);
-          }
-          if (ud > 1.e-8){
-            for (int sk=0; sk<NUM_MAGNETIC; sk++) {
-                const double perp = (qy*MX(sk) - qx*MY(sk));
-                pvec[slds[sk]] = perp*py*ud;
-            }
-            scattering += CALL_IQ(q, q_index, local_values);
-            for (int sk=0; sk<NUM_MAGNETIC; sk++) {
-                pvec[slds[sk]] = MZ(sk)*pz*ud;
-            }
-            scattering += CALL_IQ(q, q_index, local_values);
-          }
-          if (du > 1.e-8) {
-            for (int sk=0; sk<NUM_MAGNETIC; sk++) {
-                const double perp = (qy*MX(sk) - qx*MY(sk));
-                pvec[slds[sk]] = perp*py*du;
-            }
-            scattering += CALL_IQ(q, q_index, local_values);
-            for (int sk=0; sk<NUM_MAGNETIC; sk++) {
-                pvec[slds[sk]] = -MZ(sk)*pz*du;
-            }
-            scattering += CALL_IQ(q, q_index, local_values);
-          }
+        }
 #else  // !MAGNETIC
         const double scattering = CALL_IQ(q, q_index, local_values);