Changeset ff10479 in sasmodels for sasmodels/kernel_iq.c

Timestamp:

Nov 2, 2017 4:27:16 PM (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:

9d9fcbd

Parents:

cf8b630

git-author:

Paul Kienzle <pkienzle@…> (11/01/17 18:16:45)

git-committer:

Paul Kienzle <pkienzle@…> (11/02/17 16:27:16)

Message:

set projection used in theory to that selected in jitter.py

File:

• sasmodels/kernel_iq.c (modified) (8 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
+//  PROJECTION : equirectangular=1, sinusoidal=2
+//      see explore/jitter.py for definitions.
 
 #ifndef _PAR_BLOCK_ // protected block so we can include this code twice.
@@ -457 +449 @@
 // capture the differences in one spot so the rest of the code is easier
 // to read.
+
 #if defined(CALL_IQ)
   // 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)
@@ -469 +461 @@
   // 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)
@@ -481 +472 @@
   double qa, qc;
   QACRotation rotation;
-  // Grab the "view" angles (theta, phi) from the initial parameter table.
-  // These are separate from the "jitter" angles (dtheta, dphi) which are
-  // stored with the dispersity values and copied to the local parameter
-  // table as we go through the mesh.
-  const double theta = values[details->theta_par+2];
-  const double phi = values[details->theta_par+3];
-  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
+  // theta, phi, dtheta, dphi are defined below in projection to avoid repeated code.
   #define BUILD_ROTATION() qac_rotation(&rotation, theta, phi, dtheta, dphi);
   #define APPLY_ROTATION() qac_apply(rotation, qx, qy, &qa, &qc)
@@ -514 +483 @@
   double qa, qb, qc;
   QABCRotation rotation;
+  // theta, phi, dtheta, dphi are defined below in projection to avoid repeated code.
+  // psi and dpsi are only for IQ_ABC, so they are processed here.
+  const double psi = values[details->theta_par+4];
+  #define BUILD_ROTATION() qabc_rotation(&rotation, theta, phi, psi, dtheta, dphi, local_values.table.psi)
+  #define APPLY_ROTATION() qabc_apply(rotation, qx, qy, &qa, &qb, &qc)
+  #define CALL_KERNEL() CALL_IQ_ABC(qa, qb, qc, local_values.table)
+#endif
+
+// Doing jitter projection code outside the previous if block so that we don't
+// need to repeat the identical logic in the IQ_AC and IQ_ABC branches.  This
+// will become more important if we implement more projections, or more
+// complicated projections.
+#if defined(CALL_IQ) || defined(CALL_IQ_A)
+  #define APPLY_PROJECTION() const double weight=weight0
+#else // !spherosymmetric projection
   // Grab the "view" angles (theta, phi, psi) from the initial parameter table.
-  // These are separate from the "jitter" angles (dtheta, dphi, dpsi) which are
-  // stored with the dispersity values and copied to the local parameter
-  // table as we go through the mesh.
   const double theta = values[details->theta_par+2];
   const double phi = values[details->theta_par+3];
-  const double psi = values[details->theta_par+4];
-  double dtheta, dphi, dpsi, weight;
-  #if USE_SCALED_DPHI
-    #define SCALE_DPHI_AND_SET_WEIGHT() do { \
+  // The "jitter" angles (dtheta, dphi, dpsi) are stored with the
+  // dispersity values and copied to the local parameter table as
+  // we go through the mesh.
+  double dtheta, dphi, weight;
+  #if PROJECTION == 1
+    #define APPLY_PROJECTION() 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)
+  #elif PROJECTION == 2
+    #define APPLY_PROJECTION() do { \
+      dtheta = local_values.table.theta; \
+      dphi = local_values.table.phi; \
       weight = weight0; \
-      if (dtheta != 90.0) dphi /= fabs(cos(dtheta*M_PI_180)); \
+      if (dtheta != 90.0) dphi /= 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)
-
-#endif
+#endif // !spherosymmetric projection
+
 
 // ====== construct the loops =======
@@ -596 +574 @@
 
 //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 =======
@@ -603 +580 @@
   #endif
   {
+     APPLY_PROJECTION();
+
     // Accumulate I(q)
     // Note: weight==0 must always be excluded
@@ -697 +676 @@
 #undef PD_OPEN
 #undef PD_CLOSE
-#undef SCALE_DPHI_AND_SET_WEIGHT
 #undef FETCH_Q
+#undef APPLY_PROJECTION
 #undef BUILD_ROTATION
 #undef APPLY_ROTATION