Changes in / [9c490bc:19bf626] in sasmodels

File:

• sasmodels/kernel_iq.c (modified) (10 diffs)

sasmodels/kernel_iq.c

@@ -40 +40 @@
            double sld_solvent
 
-   PARAMETER_CALL(var) is the declaration of a call to the kernel.
+   CALL_IQ(q, nq, i, pars) is the declaration of a call to the kernel.
 
        Cylinder:
 
-           #define PARAMETER_CALL(var) \
+           #define CALL_IQ(q, nq, i, var) \
+           Iq(q[i], \
            var.length, \
            var.radius, \
            var.sld, \
-           var.sld_solvent
+           var.sld_solvent)
 
        Multi-shell cylinder:
-           #define PARAMETER_CALL(var) \
+           #define CALL_IQ(q, nq, i, var) \
+           Iq(q[i], \
            var.num_shells, \
            var.length, \
            var.radius, \
            var.sld, \
-           var.sld_solvent
+           var.sld_solvent)
+
+   CALL_VOLUME(var) is similar, but for calling the form volume.
 
    INVALID is a test for model parameters in the correct range
@@ -104 +108 @@
         par_offset = {2, 3, 303, 313}  // parameter offsets
         par_coord = {0, 3, 2, 1} // bitmap of parameter dependencies
+        fast_coord_index = {5, 3, x, x}
         fast_coord_count = 2  // two parameters vary with *length* distribution
-        fast_coord_index = {5, 3, x, x, x, x}
+        theta_var = -1   // no spherical correction
+        fast_theta = 0   // spherical correction angle is not pd 1
     }
 
@@ -197 +203 @@
 Scale and background cannot be coordinated with other polydisperse parameters
 
-Cutoff paramater is basically used to restrict the region where integration
-is peformed i.e. polydispersity hypercude is limitted spheres.
-
+Oriented objects in 2-D need a spherical correction on the angular variation
+in order to preserve the 'surface area' of the weight distribution.
+
+TODO: cutoff
 */
 
@@ -213 +220 @@
     int par_offset[NPARS];  // offset of par values in the par & weight vector
     int par_coord[NPARS];   // polydispersity coordination bitvector
+    int fast_coord_index[NPARS]; // index of the fast coordination parameters
     int fast_coord_count;   // number of parameters coordinated with pd 1
-    int fast_coord_index[NPARS]; // index of the fast coordination parameters
+    int theta_var;          // id of spherical correction variable
+    int fast_theta;         // true if spherical correction depends on pd 1
 } ProblemDetails;
 
@@ -317 +326 @@
   double weight = 1.0;        // set to 1 in case there are no weights
   double vol_weight = 1.0;    // set to 1 in case there are no vol weights
+  double spherical_correction = 1.0;  // correction for latitude variation
 
   // Loop over the weights then loop over q, accumulating values
@@ -322 +332 @@
     // check if indices need to be updated
     if (pd_index[0] >= pd_length[0]) {
+
+      // RESET INDICES
       pd_index[0] = loop_index%pd_length[0];
       partial_weight = 1.0;
@@ -331 +343 @@
         if (pd_isvol[k]) partial_volweight *= wi;
       }
-      weight = partial_weight * weights[pd_offset[0]+pd_index[0]]
       for (int k=0; k < NPARS; k++) {
         int coord = par_coord[k];
@@ -346 +357 @@
         parvec[k] = pars[this_offset];
       }
+      weight = partial_weight * weights[pd_offset[0]+pd_index[0]]
+      if (theta_var >= 0) {
+        spherical_correction = fabs(cos(M_PI_180*parvec[theta_var]));
+      }
+      if (!fast_theta) weight *= spherical_correction;
+
     } else {
+
+      // INCREMENT INDICES
       pd_index[0] += 1;
       const double wi = weights[pd_offset[0]+pd_index[0]];
@@ -355 +374 @@
             = pars[offset[fast_coord_index[k]] + pd_index[0]];
       }
-    }
+      if (fast_theta) weight *= fabs(cos(M_PI_180*parvec[theta_var]));
+
+    }
+
     #ifdef INVALID
     if (INVALID(local_pars)) continue;
     #endif
+
     if (weight > cutoff) {
       norm += weight;
-      vol += vol_weight * form_volume(VOLUME_PARAMETERS);
+      vol += vol_weight * CALL_VOLUME(local_pars);
       norm_vol += vol_weight;
 
@@ -369 +392 @@
       for (int i=0; i < nq; i++) {
       {
-        const double scattering = IQ_FUNC(IQ_PARS, IQ_PARAMETERS);
+        const double scattering = CALL_IQ(q, nq, i, local_pars);
         //const double scattering = Iq(q[i], IQ_PARAMETERS);
         result[i] += weight*scattering;