Changes in sasmodels/kernel_iq.c [c57ee9e:70530778] in sasmodels

File:

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

sasmodels/kernel_iq.c

@@ -27 +27 @@
 //      parameters in the parameter table.
 //  CALL_VOLUME(table) : call the form volume function
-//  CALL_EFFECTIVE_RADIUS(type, table) : call the R_eff function
 //  CALL_IQ(q, table) : call the Iq function for 1D calcs.
 //  CALL_IQ_A(q, table) : call the Iq function with |q| for 2D data.
-//  CALL_FQ(q, F1, F2, table) : call the Fq function for 1D calcs.
-//  CALL_FQ_A(q, F1, F2, table) : call the Iq function with |q| for 2D data.
 //  CALL_IQ_AC(qa, qc, table) : call the Iqxy function for symmetric shapes
 //  CALL_IQ_ABC(qa, qc, table) : call the Iqxy function for asymmetric shapes
@@ -39 +36 @@
 //  PROJECTION : equirectangular=1, sinusoidal=2
 //      see explore/jitter.py for definitions.
-
 
 #ifndef _PAR_BLOCK_ // protected block so we can include this code twice.
@@ -192 +188 @@
     QACRotation *rotation,
     double qx, double qy,
-    double *qa_out, double *qc_out)
+    double *qab_out, double *qc_out)
 {
+    // Indirect calculation of qab, from qab^2 = |q|^2 - qc^2
     const double dqc = rotation->R31*qx + rotation->R32*qy;
-    // Indirect calculation of qab, from qab^2 = |q|^2 - qc^2
-    const double dqa = sqrt(-dqc*dqc + qx*qx + qy*qy);
-
-    *qa_out = dqa;
+    const double dqab_sq = -dqc*dqc + qx*qx + qy*qy;
+    //*qab_out = sqrt(fabs(dqab_sq));
+    *qab_out = dqab_sq > 0.0 ? sqrt(dqab_sq) : 0.0;
     *qc_out = dqc;
 }
@@ -274 +270 @@
 #endif // _QABC_SECTION
 
+
 // ==================== KERNEL CODE ========================
-#define COMPUTE_F1_F2 defined(CALL_FQ)
+
 kernel
 void KERNEL_NAME(
@@ -285 +282 @@
     global const double *q, // nq q values, with padding to boundary
     global double *result,  // nq+1 return values, again with padding
-    const double cutoff,     // cutoff in the dispersity weight product
-    int32_t effective_radius_type // which effective radius to compute
+    const double cutoff     // cutoff in the dispersity weight product
     )
 {
@@ -345 +341 @@
   // version must loop over all q.
   #ifdef USE_OPENCL
-    #if COMPUTE_F1_F2
-      double weight_norm = (pd_start == 0 ? 0.0 : result[2*nq]);
-      double weighted_volume = (pd_start == 0 ? 0.0 : result[2*nq+1]);
-      double weighted_radius = (pd_start == 0 ? 0.0 : result[2*nq+2]);
-      double this_F2 = (pd_start == 0 ? 0.0 : result[2*q_index+0]);
-      double this_F1 = (pd_start == 0 ? 0.0 : result[2*q_index+1]);
-    #else
-      double weight_norm = (pd_start == 0 ? 0.0 : result[nq]);
-      double weighted_volume = (pd_start == 0 ? 0.0 : result[nq+1]);
-      double weighted_radius = (pd_start == 0 ? 0.0 : result[nq+2]);
-      double this_result = (pd_start == 0 ? 0.0 : result[q_index]);
-    #endif
+    double pd_norm = (pd_start == 0 ? 0.0 : result[nq]);
+    double this_result = (pd_start == 0 ? 0.0 : result[q_index]);
   #else // !USE_OPENCL
-    #if COMPUTE_F1_F2
-      double weight_norm = (pd_start == 0 ? 0.0 : result[2*nq]);
-      double weighted_volume = (pd_start == 0 ? 0.0 : result[2*nq+1]);
-      double weighted_radius = (pd_start == 0 ? 0.0 : result[2*nq+2]);
-    #else
-      double weight_norm = (pd_start == 0 ? 0.0 : result[nq]);
-      double weighted_volume = (pd_start == 0 ? 0.0 : result[nq+1]);
-      double weighted_radius = (pd_start == 0 ? 0.0 : result[nq+2]);
-    #endif
+    double pd_norm = (pd_start == 0 ? 0.0 : result[nq]);
     if (pd_start == 0) {
       #ifdef USE_OPENMP
       #pragma omp parallel for
       #endif
-      #if COMPUTE_F1_F2
-          // 2*nq for F^2,F pairs
-          for (int q_index=0; q_index < 2*nq; q_index++) result[q_index] = 0.0;
-      #else
-          for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0;
-      #endif
+      for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0;
     }
-    //if (q_index==0) printf("start %d %g %g\n", pd_start, weighted_volume, result[0]);
+    //if (q_index==0) printf("start %d %g %g\n", pd_start, pd_norm, result[0]);
 #endif // !USE_OPENCL
 
@@ -469 +442 @@
 // inner loop and defines the macros that use them.
 
-
-#if defined(CALL_FQ) // COMPUTE_F1_F2 is true
-  // unoriented 1D returning <F> and <F^2>
-  double qk;
-  double F1, F2;
-  #define FETCH_Q() do { qk = q[q_index]; } while (0)
-  #define BUILD_ROTATION() do {} while(0)
-  #define APPLY_ROTATION() do {} while(0)
-  #define CALL_KERNEL() CALL_FQ(qk,F1,F2,local_values.table)
-
-#elif defined(CALL_FQ_A)
-  // unoriented 2D return <F> and <F^2>
-  double qx, qy;
-  double F1, F2;
-  #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0)
-  #define BUILD_ROTATION() do {} while(0)
-  #define APPLY_ROTATION() do {} while(0)
-  #define CALL_KERNEL() CALL_FQ_A(sqrt(qx*qx+qy*qy),F1,F2,local_values.table)
-
-#elif defined(CALL_IQ)
-  // unoriented 1D return <F^2>
+#if defined(CALL_IQ)
+  // unoriented 1D
   double qk;
   #define FETCH_Q() do { qk = q[q_index]; } while (0)
   #define BUILD_ROTATION() do {} while(0)
   #define APPLY_ROTATION() do {} while(0)
-  #define CALL_KERNEL() CALL_IQ(qk,local_values.table)
+  #define CALL_KERNEL() CALL_IQ(qk, local_values.table)
 
 #elif defined(CALL_IQ_A)
@@ -528 +482 @@
   #define APPLY_ROTATION() qabc_apply(&rotation, qx, qy, &qa, &qb, &qc)
   #define CALL_KERNEL() CALL_IQ_ABC(qa, qb, qc, local_values.table)
-
 #elif defined(CALL_IQ_XY)
   // direct call to qx,qy calculator
@@ -542 +495 @@
 // 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) || defined(CALL_FQ) || defined(CALL_FQ_A)
-  // no orientation
+#if defined(CALL_IQ) || defined(CALL_IQ_A)  // no orientation
   #define APPLY_PROJECTION() const double weight=weight0
 #elif defined(CALL_IQ_XY) // pass orientation to the model
@@ -693 +645 @@
     // Note: weight==0 must always be excluded
     if (weight > cutoff) {
-      weighted_volume += weight * CALL_VOLUME(local_values.table);
-      #if COMPUTE_F1_F2
-      weight_norm += weight;
-      #endif
-      if (effective_radius_type != 0) {
-        weighted_radius += weight * CALL_EFFECTIVE_RADIUS(effective_radius_type, local_values.table);
-      }
+      pd_norm += weight * CALL_VOLUME(local_values.table);
       BUILD_ROTATION();
 
@@ -747 +693 @@
           }
         #else  // !MAGNETIC
-          #if COMPUTE_F1_F2
-            CALL_KERNEL(); // sets F1 and F2 by reference
-          #else
-            const double scattering = CALL_KERNEL();
-          #endif
+          const double scattering = CALL_KERNEL();
         #endif // !MAGNETIC
 //printf("q_index:%d %g %g %g %g\n", q_index, scattering, weight0);
 
         #ifdef USE_OPENCL
-          #if COMPUTE_F1_F2
-            this_F2 += weight * F2;
-            this_F1 += weight * F1;
-          #else
-            this_result += weight * scattering;
-          #endif
+          this_result += weight * scattering;
         #else // !USE_OPENCL
-          #if COMPUTE_F1_F2
-            result[2*q_index+0] += weight * F2;
-            result[2*q_index+1] += weight * F1;
-          #else
-            result[q_index] += weight * scattering;
-          #endif
+          result[q_index] += weight * scattering;
         #endif // !USE_OPENCL
       }
     }
   }
+
 // close nested loops
 ++step;
@@ -793 +726 @@
 // Remember the current result and the updated norm.
 #ifdef USE_OPENCL
-  #if COMPUTE_F1_F2
-    result[2*q_index+0] = this_F2;
-    result[2*q_index+1] = this_F1;
-    if (q_index == 0) {
-      result[2*nq+0] = weight_norm;
-      result[2*nq+1] = weighted_volume;
-      result[2*nq+2] = weighted_radius;
-    }
-  #else
-    result[q_index] = this_result;
-    if (q_index == 0) {
-      result[nq+0] = weight_norm;
-      result[nq+1] = weighted_volume;
-      result[nq+2] = weighted_radius;
-    }
-  #endif
-
-//if (q_index == 0) printf("res: %g/%g\n", result[0], weigthed_volume);
+  result[q_index] = this_result;
+  if (q_index == 0) result[nq] = pd_norm;
+//if (q_index == 0) printf("res: %g/%g\n", result[0], pd_norm);
 #else // !USE_OPENCL
-  #if COMPUTE_F1_F2
-    result[2*nq] = weight_norm;
-    result[2*nq+1] = weighted_volume;
-    result[2*nq+2] = weighted_radius;
-  #else
-    result[nq] = weight_norm;
-    result[nq+1] = weighted_volume;
-    result[nq+2] = weighted_radius;
-  #endif
-//printf("res: %g/%g\n", result[0], weighted_volume);
+  result[nq] = pd_norm;
+//printf("res: %g/%g\n", result[0], pd_norm);
 #endif // !USE_OPENCL
 
 // ** clear the macros in preparation for the next kernel **
-#undef COMPUTE_F1_F2
 #undef PD_INIT
 #undef PD_OPEN