Index: sasmodels/kernel_iq.c
===================================================================
--- sasmodels/kernel_iq.c	(revision ba32cddebeff04b66d828ca58b1843244d4425fb)
+++ sasmodels/kernel_iq.c	(revision 5ff1b037170e74a37b38e653c5baeaa8cf54e08c)
@@ -21,10 +21,11 @@
     int32_t pd_offset[MAX_PD];  // offset of pd weights in the value & weight vector
     int32_t pd_stride[MAX_PD];  // stride to move to the next index at this level
-    int32_t pd_isvol[MAX_PD];   // True if parameter is a volume weighting parameter
 #endif // MAX_PD > 0
-    int32_t par_offset[NPARS];  // offset of par values in the value & weight vector
-    int32_t par_coord[NPARS];   // polydispersity coordination bitvector
-    int32_t fast_coord_pars[NPARS]; // ids of the fast coordination parameters
-    int32_t fast_coord_count;   // number of parameters coordinated with pd 1
+    int32_t par_offset[NPARS];  // offset of par value blocks in the value & weight vector
+    int32_t par_coord[NPARS];   // ids of the coordination parameters
+    int32_t pd_coord[NPARS];    // polydispersity coordination bitvector
+    int32_t num_active;         // number of non-trivial pd loops
+    int32_t total_pd;           // total number of voxels in hypercube
+    int32_t num_coord;          // number of coordinated parameters
     int32_t theta_par;          // id of spherical correction variable
 } ProblemDetails;
@@ -54,19 +55,40 @@
   double *pvec = (double *)(&local_values);  // Alias named parameters with a vector
 
-  // Monodisperse computation
-  if (pd_stop == 1) {
-    // Shouldn't need to copy!!
-    for (int k=0; k < NPARS; k++) {
-      pvec[k] = values[k+2];  // skip scale and background
-    }
-
-    const double volume = CALL_VOLUME(local_values);
+  // Fill in the initial variables
+  #ifdef USE_OPENMP
+  #pragma omp parallel for
+  #endif
+  for (int k=0; k < NPARS; k++) {
+    pvec[k] = values[problem->par_offset[k]];
+  }
+
+  // If it is the first round initialize the result to zero, otherwise
+  // assume that the previous result has been passed back.
+  // Note: doing this even in the monodisperse case in order to handle the
+  // rare case where the model parameters are invalid and zero is returned.
+  // So slightly increased cost for slightly smaller code size.
+  if (pd_start == 0) {
     #ifdef USE_OPENMP
     #pragma omp parallel for
     #endif
+    for (int i=0; i < nq+1; i++) {
+      result[i] = 0.0;
+    }
+  }
+
+  // Monodisperse computation
+  if (problem->num_active == 0) {
+    #ifdef INVALID
+    if (INVALID(local_values)) { return; }
+    #endif
+
+    const double norm = CALL_VOLUME(local_values);
+    #ifdef USE_OPENMP
+    #pragma omp parallel for
+    #endif
+    result[nq] = norm; // Total volume normalization
     for (int i=0; i < nq; i++) {
       double scattering = CALL_IQ(q, i, local_values);
-      if (volume != 0.0) scattering /= volume;
-      result[i] = values[0]*scattering + values[1];
+      result[i] = values[0]*scattering/norm + values[1];
     }
     return;
@@ -74,110 +96,89 @@
 
 #if MAX_PD > 0
-  //printf("Entering polydispersity\n");
-
+  //printf("Entering polydispersity from %d to %d\n", pd_start, pd_stop);
   // Since we are no longer looping over the entire polydispersity hypercube
-  // for each q, we need to track the normalization values for each q in a
-  // separate work vector.
-  double norm;   // contains sum over weights
-  double vol; // contains sum over volume
-  double norm_vol; // contains weights over volume
-
-  // Initialize the results to zero
-  if (pd_start == 0) {
-    norm_vol = 0.0;
-    norm = 0.0;
-    vol = 0.0;
-
-    #ifdef USE_OPENMP
-    #pragma omp parallel for
-    #endif
-    for (int i=0; i < nq; i++) {
-      result[i] = 0.0;
-    }
-  } else {
-    //Pulling values from previous segment
-    norm = result[nq];
-    vol = result[nq+1];
-    norm_vol = result[nq+2];
-  }
-
-  // Location in the polydispersity hypercube, one index per dimension.
-  local int pd_index[MAX_PD];
-
-  // polydispersity loop index positions
-  local int offset[NPARS];  // NPARS excludes scale/background
-
-  // Trigger the reset behaviour that happens at the end the fast loop
-  // by setting the initial index >= weight vector length.
-  pd_index[0] = problem->pd_length[0];
-
+  // for each q, we need to track the normalization values between calls.
+  double norm = 0.0;
 
   // need product of weights at every Iq calc, so keep product of
   // weights from the outer loops so that weight = partial_weight * fast_weight
   double partial_weight = NAN; // product of weight w4*w3*w2 but not w1
-  double partial_volweight = NAN;
-  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
+  double spherical_correction = 1.0;  // cosine correction for latitude variation
+
+  // Location in the polydispersity hypercube, one index per dimension.
+  local int pd_index[MAX_PD];
+
+  // Location of the coordinated parameters in their own sub-cubes.
+  local int offset[NPARS];
+
+  // Trigger the reset behaviour that happens at the end the fast loop
+  // by setting the initial index >= weight vector length.
+  const int fast_length = problem->pd_length[0];
+  pd_index[0] = fast_length;
 
   // Loop over the weights then loop over q, accumulating values
   for (int loop_index=pd_start; loop_index < pd_stop; loop_index++) {
     // check if indices need to be updated
-    if (pd_index[0] >= problem->pd_length[0]) {
-
-      // RESET INDICES
-      pd_index[0] = loop_index%problem->pd_length[0];
+    if (pd_index[0] == fast_length) {
+      //printf("should be here with %d active\n", problem->num_active);
+
+      // Compute position in polydispersity hypercube
+      for (int k=0; k < problem->num_active; k++) {
+        pd_index[k] = (loop_index/problem->pd_stride[k])%problem->pd_length[k];
+        //printf("pd_index[%d] = %d\n",k,pd_index[k]);
+      }
+
+      // Compute partial weights
       partial_weight = 1.0;
-      partial_volweight = 1.0;
-      for (int k=1; k < MAX_PD; k++) {
-        pd_index[k] = (loop_index%problem->pd_length[k])/problem->pd_stride[k];
-        const double wi = weights[problem->pd_offset[k]+pd_index[k]];
+      //printf("partial weight %d: ", loop_index);
+      for (int k=1; k < problem->num_active; k++) {
+        double wi = weights[problem->pd_offset[k] + pd_index[k]];
+        //printf("pd[%d]=par[%d]=%g ", k, problem->pd_par[k], wi);
         partial_weight *= wi;
-        if (problem->pd_isvol[k]) partial_volweight *= wi;
-      }
+      }
+      //printf("\n");
+
+      // Update parameter offsets in weight vector
       //printf("slow %d: ", loop_index);
-      for (int k=0; k < NPARS; k++) {
-        int coord = problem->par_coord[k];
-        int this_offset = problem->par_offset[k];
+      for (int k=0; k < problem->num_coord; k++) {
+        int par = problem->par_coord[k];
+        int coord = problem->pd_coord[k];
+        int this_offset = problem->par_offset[par];
         int block_size = 1;
-        for (int bit=0; bit < MAX_PD && coord != 0; bit++) {
+        for (int bit=0; coord != 0; bit++) {
           if (coord&1) {
               this_offset += block_size * pd_index[bit];
               block_size *= problem->pd_length[bit];
           }
-          coord /= 2;
+          coord >>= 1;
         }
-        offset[k] = this_offset;
-        pvec[k] = values[this_offset];
-        //printf("p[%d]=v[%d]=%g ", k, offset[k], pvec[k]);
+        offset[par] = this_offset;
+        pvec[par] = values[this_offset];
+        //printf("par[%d]=v[%d]=%g \n", k, offset[k], pvec[k]);
+        // if theta is not coordinated with fast index, precompute spherical correction
+        if (par == problem->theta_par && !(problem->par_coord[k]&1)) {
+          spherical_correction = fmax(fabs(cos(M_PI_180*pvec[problem->theta_par])), 1e-6);
+        }
       }
       //printf("\n");
-      weight = partial_weight * weights[problem->pd_offset[0]+pd_index[0]];
-      if (problem->theta_par >= 0) {
-        spherical_correction = fabs(cos(M_PI_180*pvec[problem->theta_par]));
-      }
-      if (problem->theta_par == problem->pd_par[0]) {
-        weight *= spherical_correction;
-      }
-      pd_index[0] += 1;
-
-    } else {
-
-      // INCREMENT INDICES
-      const double wi = weights[problem->pd_offset[0]+pd_index[0]];
-      weight = partial_weight*wi;
-      if (problem->pd_isvol[0]) vol_weight *= wi;
-      //printf("fast %d: ", loop_index);
-      for (int k=0; k < problem->fast_coord_count; k++) {
-        const int pindex = problem->fast_coord_pars[k];
-        pvec[pindex] = values[++offset[pindex]];
-        //printf("p[%d]=v[%d]=%g ", pindex, offset[pindex], pvec[pindex]);
-      }
-      //printf("\n");
-      if (problem->theta_par == problem->pd_par[0]) {
-        weight *= fabs(cos(M_PI_180*pvec[problem->theta_par]));
-      }
-      pd_index[0] += 1;
-    }
+    }
+
+    // Increment fast index
+    const double wi = weights[problem->pd_offset[0] + pd_index[0]++];
+    double weight = partial_weight*wi;
+    //printf("fast %d: ", loop_index);
+    for (int k=0; k < problem->num_coord; k++) {
+      if (problem->pd_coord[k]&1) {
+        const int par = problem->par_coord[k];
+        pvec[par] = values[offset[par]++];
+        //printf("p[%d]=v[%d]=%g ", par, offset[par]-1, pvec[par]);
+        // if theta is coordinated with fast index, compute spherical correction each time
+        if (par == problem->theta_par) {
+          spherical_correction = fmax(fabs(cos(M_PI_180*pvec[problem->theta_par])), 1e-6);
+        }
+      }
+    }
+    //printf("\n");
+
     #ifdef INVALID
     if (INVALID(local_values)) continue;
@@ -187,7 +188,8 @@
     // Note: weight==0 must always be excluded
     if (weight > cutoff) {
-      norm += weight;
-      vol += vol_weight * CALL_VOLUME(local_values);
-      norm_vol += vol_weight;
+      // spherical correction has some nasty effects when theta is +90 or -90
+      // where it becomes zero.  If the entirety of the correction
+      weight *= spherical_correction;
+      norm += weight * CALL_VOLUME(local_values);
 
       #ifdef USE_OPENMP
@@ -202,17 +204,12 @@
 
   // Make normalization available for the next round
-  result[nq] = norm;
-  result[nq+1] = vol;
-  result[nq+2] = norm_vol;
+  result[nq] += norm;
 
   // End of the PD loop we can normalize
-  if (pd_stop >= problem->pd_stride[MAX_PD-1]) {
+  if (pd_stop >= problem->total_pd) {
     #ifdef USE_OPENMP
     #pragma omp parallel for
     #endif
     for (int i=0; i < nq; i++) {
-      if (vol*norm_vol != 0.0) {
-        result[i] *= norm_vol/vol;
-      }
       result[i] = values[0]*result[i]/norm + values[1];
     }
