Changeset ff10479 in sasmodels

Timestamp:

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

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

Files:

• explore/jitter.py (modified) (6 diffs)
• sasmodels/generate.py (modified) (2 diffs)
• sasmodels/kernel_iq.c (modified) (8 diffs)

explore/jitter.py

@@ -162 +162 @@
 
 PROJECTIONS = [
-    'equirectangular', 'azimuthal_equidistance', 'sinusoidal', 'guyou',
+    # in order of PROJECTION number; do not change without updating the
+    # constants in kernel_iq.c
+    'equirectangular', 'sinusoidal', 'guyou', 'azimuthal_equidistance',
 ]
 def draw_mesh(ax, view, jitter, radius=1.2, n=11, dist='gaussian',
@@ -239 +241 @@
         raise ValueError("expected dist to be gaussian, rectangle or uniform")
 
-    if projection == 'equirectangular':
+    if projection == 'equirectangular':  #define PROJECTION 1
         def rotate(theta_i, phi_j):
             return Rx(phi_j)*Ry(theta_i)
         def weight(theta_i, phi_j, wi, wj):
             return wi*wj*abs(cos(radians(theta_i)))
-    elif projection == 'azimuthal_equidistance':
+    elif projection == 'sinusoidal':  #define PROJECTION 2
+        def rotate(theta_i, phi_j):
+            latitude = theta_i
+            scale = cos(radians(latitude))
+            longitude = phi_j/scale if abs(phi_j) < abs(scale)*180 else 0
+            #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude))
+            return Rx(longitude)*Ry(latitude)
+        def weight(theta_i, phi_j, wi, wj):
+            latitude = theta_i
+            scale = cos(radians(latitude))
+            w = 1 if abs(phi_j) < abs(scale)*180 else 0
+            return w*wi*wj
+    elif projection == 'guyou':  #define PROJECTION 3  (eventually?)
+        def rotate(theta_i, phi_j):
+            from guyou import guyou_invert
+            #latitude, longitude = guyou_invert([theta_i], [phi_j])
+            longitude, latitude = guyou_invert([phi_j], [theta_i])
+            return Rx(longitude[0])*Ry(latitude[0])
+        def weight(theta_i, phi_j, wi, wj):
+            return wi*wj
+    elif projection == 'azimuthal_equidistance':  # Note: Rz Ry, not Rx Ry
         def rotate(theta_i, phi_j):
             latitude = sqrt(theta_i**2 + phi_j**2)
@@ -282 +304 @@
             w = sin(radians(latitude))/latitude if latitude != 0 else 1
             return w*wi*wj if latitude < 180 else 0
-    elif projection == 'sinusoidal':
-        def rotate(theta_i, phi_j):
-            latitude = theta_i
-            scale = cos(radians(latitude))
-            longitude = phi_j/scale if abs(phi_j) < abs(scale)*180 else 0
-            #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude))
-            return Rx(longitude)*Ry(latitude)
-        def weight(theta_i, phi_j, wi, wj):
-            latitude = theta_i
-            scale = cos(radians(latitude))
-            w = 1 if abs(phi_j) < abs(scale)*180 else 0
-            return w*wi*wj
     elif projection == 'azimuthal_equal_area':
         def rotate(theta_i, phi_j):
@@ -305 +315 @@
             w = sin(radians(latitude))/latitude if latitude != 0 else 1
             return w*wi*wj if latitude < 180 else 0
-    elif projection == 'guyou':
-        def rotate(theta_i, phi_j):
-            from guyou import guyou_invert
-            #latitude, longitude = guyou_invert([theta_i], [phi_j])
-            longitude, latitude = guyou_invert([phi_j], [theta_i])
-            return Rx(longitude[0])*Ry(latitude[0])
-        def weight(theta_i, phi_j, wi, wj):
-            return wi*wj
     else:
         raise ValueError("unknown projection %r"%projection)
@@ -490 +492 @@
         vmin, vmax = calculator.limits
     else:
-        vmin, vmax = clipped_range(Iqxy, portion=0.95, mode='top')
+        vmax = Iqxy.max()
+        vmin = vmax*10**-7
+        #vmin, vmax = clipped_range(Iqxy, portion=portion, mode='top')
     #print("range",(vmin,vmax))
     #qx, qy = np.meshgrid(qx, qy)
@@ -605 +609 @@
     *model_name* is one of the models available in :func:`select_model`.
     """
+    # projection number according to 1-order position in list, but
+    # only 1 and 2 are implemented so far.
+    from sasmodels import generate
+    generate.PROJECTION = PROJECTIONS.index(projection) + 1
+    if generate.PROJECTION > 2:
+        print("*** PROJECTION %s not implemented in scattering function ***"%projection)
+        generate.PROJECTION = 2
+
     # set up calculator
     calculator, size = select_calculator(model_name, n=150, size=size)

sasmodels/generate.py

@@ -179 +179 @@
 except ImportError:
     pass
+
+# jitter projection to use in the kernel code.  See explore/jitter.py
+# for details.  To change it from a program, set generate.PROJECTION.
+PROJECTION = 1
 
 def get_data_path(external_dir, target_file):
@@ -764 +768 @@
     source.append("#define NUM_MAGNETIC %d" % partable.nmagnetic)
     source.append("#define MAGNETIC_PARS %s"%",".join(str(k) for k in magpars))
+    source.append("#define PROJECTION %d"%PROJECTION)
 
     # TODO: allow mixed python/opencl kernels?

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