Changeset 767dca8 in sasmodels for sasmodels/kernel_iq.c

Timestamp:

Oct 27, 2017 10:00:07 AM (6 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:

da5536f

Parents:

a9bc435

Message:

use alternative jitter, which maintains constant arc length for the phi distribution

File:

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

sasmodels/kernel_iq.c

@@ -32 +32 @@
 //  INVALID(table) : test if the current point is feesible to calculate.  This
 //      will be defined in the kernel definition file.
+
+// Set SCALED_PHI to 1 if dphi represents distance rather than longitude.
+// That is, in order to make an equally spaced mesh on a curved surface,
+// you will need to scale longitude displacement by cos(latitude) to account
+// for the reduction in distance per degree latitude as you move away from
+// the equator.  Points on the mesh with scaled dphi values more than 180
+// degrees are not included in the dispersity calculation. This should make
+// the integral over the entire surface work by sampling fewer points near
+// the poles.
+# define USE_SCALED_DPHI 1
 
 #ifndef _PAR_BLOCK_ // protected block so we can include this code twice.
@@ -450 +460 @@
   // unoriented 1D
   double qk;
+  #define SCALE_DPHI_AND_SET_WEIGHT() const double weight=weight0
   #define FETCH_Q() do { qk = q[q_index]; } while (0)
   #define BUILD_ROTATION() do {} while(0)
@@ -458 +469 @@
   // unoriented 2D
   double qx, qy;
+  #define SCALE_DPHI_AND_SET_WEIGHT() const double weight=weight0
   #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0)
   #define BUILD_ROTATION() do {} while(0)
@@ -475 +487 @@
   const double theta = values[details->theta_par+2];
   const double phi = values[details->theta_par+3];
-  #define BUILD_ROTATION() qac_rotation(&rotation, \
-    theta, phi, local_values.table.theta, local_values.table.phi)
+  double dtheta, dphi, weight;
+  #if USE_SCALED_DPHI
+    #define SCALE_DPHI_AND_SET_WEIGHT() do { \
+      dtheta = local_values.table.theta; \
+      dphi = local_values.table.phi; \
+      weight = weight0; \
+      if (dtheta != 90.0) dphi /= fabs(cos(dtheta*M_PI_180)); \
+      else if (dphi != 0.0) weight = 0.; \
+      if (fabs(dphi) >= 180.) weight = 0.; \
+    } while (0)
+  #else
+    #define SCALE_DPHI_AND_SET_WEIGHT() do { \
+      dtheta = local_values.table.theta; \
+      dphi = local_values.table.phi; \
+      weight = fabs(cos(dtheta*M_PI_180)) * weight0; \
+    } while (0)
+  #endif
+  #define BUILD_ROTATION() qac_rotation(&rotation, theta, phi, dtheta, dphi);
   #define APPLY_ROTATION() qac_apply(rotation, qx, qy, &qa, &qc)
   #define CALL_KERNEL() CALL_IQ_AC(qa, qc, local_values.table)
@@ -493 +521 @@
   const double phi = values[details->theta_par+3];
   const double psi = values[details->theta_par+4];
-  #define BUILD_ROTATION() qabc_rotation(&rotation, \
-    theta, phi, psi, local_values.table.theta, \
-    local_values.table.phi, local_values.table.psi)
+  double dtheta, dphi, dpsi, weight;
+  #if USE_SCALED_DPHI
+    #define SCALE_DPHI_AND_SET_WEIGHT() do { \
+      dtheta = local_values.table.theta; \
+      dphi = local_values.table.phi; \
+      dpsi = local_values.table.psi; \
+      weight = weight0; \
+      if (dtheta != 90.0) dphi /= fabs(cos(dtheta*M_PI_180)); \
+      else if (dphi != 0.0) weight = 0.; \
+      if (fabs(dphi) >= 180.) weight = 0.; \
+    } while (0)
+  #else
+    #define SCALE_DPHI_AND_SET_WEIGHT() do { \
+      dtheta = local_values.table.theta; \
+      dphi = local_values.table.phi; \
+      dpsi = local_values.table.psi; \
+      weight = fabs(cos(dtheta*M_PI_180)) * weight0; \
+    } while (0)
+  #endif
+  #define BUILD_ROTATION() qabc_rotation(&rotation, theta, phi, psi, dtheta, dphi, dpsi)
   #define APPLY_ROTATION() qabc_apply(rotation, qx, qy, &qa, &qb, &qc)
   #define CALL_KERNEL() CALL_IQ_ABC(qa, qb, qc, local_values.table)
@@ -551 +596 @@
 
 //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, local_values.vector[i]); printf("\n");}
+  SCALE_DPHI_AND_SET_WEIGHT();
 
   // ====== loop body =======
@@ -559 +605 @@
     // Accumulate I(q)
     // Note: weight==0 must always be excluded
-    if (weight0 > cutoff) {
-      pd_norm += weight0 * CALL_VOLUME(local_values.table);
+    if (weight > cutoff) {
+      pd_norm += weight * CALL_VOLUME(local_values.table);
       BUILD_ROTATION();
 
@@ -608 +654 @@
           const double scattering = CALL_KERNEL();
         #endif // !MAGNETIC
-//printf("q_index:%d %g %g %g %g\n", q_index, scattering, weight, weight0);
+//printf("q_index:%d %g %g %g %g\n", q_index, scattering, weight0);
 
         #ifdef USE_OPENCL
-          this_result += weight0 * scattering;
+          this_result += weight * scattering;
         #else // !USE_OPENCL
-          result[q_index] += weight0 * scattering;
+          result[q_index] += weight * scattering;
         #endif // !USE_OPENCL
       }
@@ -651 +697 @@
 #undef PD_OPEN
 #undef PD_CLOSE
+#undef SCALE_DPHI_AND_SET_WEIGHT
 #undef FETCH_Q
 #undef BUILD_ROTATION