Index: sasmodels/details.py
===================================================================
--- sasmodels/details.py	(revision 40a87fa5096610e24c3ead10299f609e99f093f2)
+++ sasmodels/details.py	(revision bde38b559fcdc618ebafdbb69cf7af6120ea8466)
@@ -71,9 +71,9 @@
         #   pd_offset[max_pd]  offset of pd values in parameter array
         #   pd_stride[max_pd]  index of pd value in loop = n//stride[k]
-        #   pd_prod            total length of pd loop
-        #   pd_sum             total length of the weight vector
+        #   num_eval           total length of pd loop
+        #   num_weights        total length of the weight vector
         #   num_active         number of pd params
         #   theta_par          parameter number for theta parameter
-        self.buffer = np.zeros(4*max_pd + 4, 'i4')
+        self.buffer = np.empty(4*max_pd + 4, 'i4')
 
         # generate views on different parts of the array
@@ -86,4 +86,15 @@
         self.theta_par = parameters.theta_offset
 
+        # offset and length are for all parameters, not just pd parameters
+        # They are not sent to the kernel function, though they could be.
+        # They are used by the composite models (sum and product) to
+        # figure out offsets into the combined value list.
+        self.offset = None  # type: np.ndarray
+        self.length = None  # type: np.ndarray
+
+        # keep hold of ifno show() so we can break a values vector
+        # into the individual components
+        self.info = model_info
+
     @property
     def pd_par(self):
@@ -107,20 +118,20 @@
 
     @property
-    def pd_prod(self):
+    def num_eval(self):
         """Total size of the pd mesh"""
         return self.buffer[-4]
 
-    @pd_prod.setter
-    def pd_prod(self, v):
+    @num_eval.setter
+    def num_eval(self, v):
         """Total size of the pd mesh"""
         self.buffer[-4] = v
 
     @property
-    def pd_sum(self):
+    def num_weights(self):
         """Total length of all the weight vectors"""
         return self.buffer[-3]
 
-    @pd_sum.setter
-    def pd_sum(self, v):
+    @num_weights.setter
+    def num_weights(self, v):
         """Total length of all the weight vectors"""
         self.buffer[-3] = v
@@ -146,70 +157,110 @@
         self.buffer[-1] = v
 
-    def show(self):
+    def show(self, values=None):
         """Print the polydispersity call details to the console"""
-        print("num_active", self.num_active)
-        print("pd_prod", self.pd_prod)
-        print("pd_sum", self.pd_sum)
-        print("theta par", self.theta_par)
-        print("pd_par", self.pd_par)
-        print("pd_length", self.pd_length)
-        print("pd_offset", self.pd_offset)
-        print("pd_stride", self.pd_stride)
-
-
-def mono_details(model_info):
-    # type: (ModelInfo) -> CallDetails
-    """
-    Return a :class:`CallDetails` object for a monodisperse calculation
-    of the model defined by *model_info*.
-    """
-    call_details = CallDetails(model_info)
-    call_details.pd_prod = 1
-    call_details.pd_sum = model_info.parameters.nvalues
-    call_details.pd_par[:] = np.arange(0, model_info.parameters.max_pd)
-    call_details.pd_length[:] = 1
-    call_details.pd_offset[:] = np.arange(0, model_info.parameters.max_pd)
-    call_details.pd_stride[:] = 1
-    return call_details
-
-
-def poly_details(model_info, weights):
+        print("===== %s details ===="%self.info.name)
+        print("num_active:%d  num_eval:%d  num_weights:%d  theta=%d"
+              % (self.num_active, self.num_eval, self.num_weights, self.theta_par))
+        if self.pd_par.size:
+            print("pd_par", self.pd_par)
+            print("pd_length", self.pd_length)
+            print("pd_offset", self.pd_offset)
+            print("pd_stride", self.pd_stride)
+        if values is not None:
+            nvalues = self.info.parameters.nvalues
+            print("scale, background", values[:2])
+            print("val", values[2:nvalues])
+            print("pd", values[nvalues:nvalues+self.num_weights])
+            print("wt", values[nvalues+self.num_weights:nvalues+2*self.num_weights])
+            print("offsets", self.offset)
+
+
+def make_details(model_info, length, offset, num_weights):
     """
     Return a :class:`CallDetails` object for a polydisperse calculation
-    of the model defined by *model_info* for the given set of *weights*.
-    *weights* is a list of vectors, one for each parameter.  Monodisperse
-    parameters should provide a weight vector containing [1.0].
-    """
-    # type: (ModelInfo) -> CallDetails
-    #print("weights",weights)
-    #weights = weights[2:] # Skip scale and background
-
-    # Decreasing list of polydispersity lengths
-    #print([p.id for p in model_info.parameters.call_parameters])
-    pd_length = np.array([len(w)
-                          for w in weights[2:2+model_info.parameters.npars]])
-    num_active = np.sum(pd_length > 1)
+    of the model defined by *model_info*.  Polydispersity is defined by
+    the *length* of the polydispersity distribution for each parameter
+    and the *offset* of the distribution in the polydispersity array.
+    Monodisperse parameters should use a polydispersity length of one
+    with weight 1.0. *num_weights* is the total length of the polydispersity
+    array.
+    """
+    # type: (ModelInfo, np.ndarray, np.ndarray, int) -> CallDetails
+    #pars = model_info.parameters.call_parameters[2:model_info.parameters.npars+2]
+    #print(", ".join(str(i)+"-"+p.id for i,p in enumerate(pars)))
+    #print("len:",length)
+    #print("off:",offset)
+
+    # Check that we arn't using too many polydispersity loops
+    num_active = np.sum(length > 1)
     max_pd = model_info.parameters.max_pd
     if num_active > max_pd:
         raise ValueError("Too many polydisperse parameters")
 
-    pd_offset = np.cumsum(np.hstack((0, pd_length)))
-    #print(", ".join(str(i)+"-"+p.id for i,p in enumerate(model_info.parameters.call_parameters)))
-    #print("len:",pd_length)
-    #print("off:",pd_offset)
+    # Decreasing list of polydpersity lengths
     # Note: the reversing view, x[::-1], does not require a copy
-    idx = np.argsort(pd_length)[::-1][:max_pd]
-    pd_stride = np.cumprod(np.hstack((1, pd_length[idx])))
+    idx = np.argsort(length)[::-1][:max_pd]
+    pd_stride = np.cumprod(np.hstack((1, length[idx])))
 
     call_details = CallDetails(model_info)
     call_details.pd_par[:max_pd] = idx
-    call_details.pd_length[:max_pd] = pd_length[idx]
-    call_details.pd_offset[:max_pd] = pd_offset[idx]
+    call_details.pd_length[:max_pd] = length[idx]
+    call_details.pd_offset[:max_pd] = offset[idx]
     call_details.pd_stride[:max_pd] = pd_stride[:-1]
-    call_details.pd_prod = pd_stride[-1]
-    call_details.pd_sum = sum(len(w) for w in weights)
+    call_details.num_eval = pd_stride[-1]
+    call_details.num_weights = num_weights
     call_details.num_active = num_active
+    call_details.length = length
+    call_details.offset = offset
     #call_details.show()
     return call_details
+
+
+ZEROS = tuple([0.]*31)
+def make_kernel_args(kernel, pairs):
+    # type: (Kernel, Tuple[List[np.ndarray], List[np.ndarray]]) -> Tuple[CallDetails, np.ndarray, bool]
+    """
+    Converts (value, weight) pairs into parameters for the kernel call.
+
+    Returns a CallDetails object indicating the polydispersity, a data object
+    containing the different values, and the magnetic flag indicating whether
+    any magnetic magnitudes are non-zero. Magnetic vectors (M0, phi, theta) are
+    converted to rectangular coordinates (mx, my, mz).
+    """
+    npars = kernel.info.parameters.npars
+    nvalues = kernel.info.parameters.nvalues
+    scalars = [v[0][0] for v in pairs]
+    values, weights = zip(*pairs[2:npars+2]) if npars else ((),())
+    length = np.array([len(w) for w in weights])
+    offset = np.cumsum(np.hstack((0, length)))
+    call_details = make_details(kernel.info, length, offset[:-1], offset[-1])
+    # Pad value array to a 32 value boundaryd
+    data_len = nvalues + 2*sum(len(v) for v in values)
+    extra = (32 - data_len%32)%32
+    data = np.hstack((scalars,) + values + weights + ZEROS[:extra])
+    data = data.astype(kernel.dtype)
+    is_magnetic = convert_magnetism(kernel.info.parameters, data)
+    #call_details.show()
+    return call_details, data, is_magnetic
+
+
+def convert_magnetism(parameters, values):
+    """
+    Convert magnetism values from polar to rectangular coordinates.
+
+    Returns True if any magnetism is present.
+    """
+    mag = values[parameters.nvalues-3*parameters.nmagnetic:parameters.nvalues]
+    mag = mag.reshape(-1, 3)
+    scale = mag[:,0]
+    if np.any(scale):
+        theta, phi = mag[:, 1]*pi/180., mag[:, 2]*pi/180.
+        cos_theta = cos(theta)
+        mag[:, 0] = scale*cos_theta*cos(phi)  # mx
+        mag[:, 1] = scale*sin(theta)  # my
+        mag[:, 2] = -scale*cos_theta*sin(phi)  # mz
+        return True
+    else:
+        return False
 
 
@@ -239,53 +290,2 @@
         value = pars
     return value, weight
-
-
-
-def build_details(kernel, pairs):
-    # type: (Kernel, Tuple[List[np.ndarray], List[np.ndarray]]) -> Tuple[CallDetails, np.ndarray, bool]
-    """
-    Converts (value, weight) pairs into parameters for the kernel call.
-
-    Returns a CallDetails object indicating the polydispersity, a data object
-    containing the different values, and the magnetic flag indicating whether
-    any magnetic magnitudes are non-zero. Magnetic vectors (M0, phi, theta) are
-    converted to rectangular coordinates (mx, my, mz).
-    """
-    values, weights = zip(*pairs)
-    scalars = [v[0] for v in values]
-    if all(len(w) == 1 for w in weights):
-        call_details = mono_details(kernel.info)
-        # Pad value array to a 32 value boundary
-        data_len = 3*len(scalars)
-        extra = ((data_len+31)//32)*32 - data_len
-        data = np.array(scalars+scalars+[1.]*len(scalars)+[0.]*extra,
-                        dtype=kernel.dtype)
-    else:
-        call_details = poly_details(kernel.info, weights)
-        # Pad value array to a 32 value boundary
-        data_len = len(scalars) + 2*sum(len(v) for v in values)
-        extra = ((data_len+31)//32)*32 - data_len
-        data = np.hstack(scalars+list(values)+list(weights)+[0.]*extra)
-        data = data.astype(kernel.dtype)
-    is_magnetic = convert_magnetism(kernel.info.parameters, data)
-    #call_details.show()
-    return call_details, data, is_magnetic
-
-def convert_magnetism(parameters, values):
-    """
-    Convert magnetism values from polar to rectangular coordinates.
-
-    Returns True if any magnetism is present.
-    """
-    mag = values[parameters.nvalues-3*parameters.nmagnetic:parameters.nvalues]
-    mag = mag.reshape(-1, 3)
-    scale = mag[:,0]
-    if np.any(scale):
-        theta, phi = mag[:, 1]*pi/180., mag[:, 2]*pi/180.
-        cos_theta = cos(theta)
-        mag[:, 0] = scale*cos_theta*cos(phi)  # mx
-        mag[:, 1] = scale*sin(theta)  # my
-        mag[:, 2] = -scale*cos_theta*sin(phi)  # mz
-        return True
-    else:
-        return False
Index: sasmodels/direct_model.py
===================================================================
--- sasmodels/direct_model.py	(revision 40a87fa5096610e24c3ead10299f609e99f093f2)
+++ sasmodels/direct_model.py	(revision bde38b559fcdc618ebafdbb69cf7af6120ea8466)
@@ -30,5 +30,5 @@
 from . import resolution
 from . import resolution2d
-from .details import build_details, dispersion_mesh
+from .details import make_kernel_args, dispersion_mesh
 
 try:
@@ -70,5 +70,5 @@
                 for p in parameters.call_parameters]
 
-    call_details, values, is_magnetic = build_details(calculator, vw_pairs)
+    call_details, values, is_magnetic = make_kernel_args(calculator, vw_pairs)
     #print("values:", values)
     return calculator(call_details, values, cutoff, is_magnetic)
Index: sasmodels/kernel.py
===================================================================
--- sasmodels/kernel.py	(revision 9eb36329a5763c7ea05bd24c85b4810c08c2c2e2)
+++ sasmodels/kernel.py	(revision bde38b559fcdc618ebafdbb69cf7af6120ea8466)
@@ -39,4 +39,5 @@
     info = None  # type: ModelInfo
     results = None # type: List[np.ndarray]
+    dtype = None  # type: np.dtype
 
     def __call__(self, call_details, values, cutoff, magnetic):
Index: sasmodels/kernel_iq.c
===================================================================
--- sasmodels/kernel_iq.c	(revision 0f00d95c97536e141fcc2852f5d81a3553c3dc70)
+++ sasmodels/kernel_iq.c	(revision bde38b559fcdc618ebafdbb69cf7af6120ea8466)
@@ -22,12 +22,17 @@
     int32_t pd_stride[MAX_PD];  // stride to move to the next index at this level
 #endif // MAX_PD > 0
-    int32_t pd_prod;            // total number of voxels in hypercube
-    int32_t pd_sum;             // total length of the weights vector
+    int32_t num_eval;           // total number of voxels in hypercube
+    int32_t num_weights;        // total length of the weights vector
     int32_t num_active;         // number of non-trivial pd loops
     int32_t theta_par;          // id of spherical correction variable
 } ProblemDetails;
 
+// Intel HD 4000 needs private arrays to be a multiple of 4 long
 typedef struct {
     PARAMETER_TABLE
+} ParameterTable;
+typedef union {
+    ParameterTable table;
+    double vector[4*((NUM_PARS+3)/4)];
 } ParameterBlock;
 #endif // _PAR_BLOCK_
@@ -79,12 +84,10 @@
     )
 {
-
   // Storage for the current parameter values.  These will be updated as we
-  // walk the polydispersity cube.  local_values will be aliased to pvec.
+  // walk the polydispersity cube.
   ParameterBlock local_values;
-  double *pvec = (double *)&local_values;
 
 #if defined(MAGNETIC) && NUM_MAGNETIC>0
-  // Location of the sld parameters in the parameter pvec.
+  // Location of the sld parameters in the parameter vector.
   // These parameters are updated with the effective sld due to magnetism.
   #if NUM_MAGNETIC > 3
@@ -110,25 +113,22 @@
   #endif
   for (int i=0; i < NUM_PARS; i++) {
-    pvec[i] = values[2+i];
-//printf("p%d = %g\n",i, pvec[i]);
-  }
-
-  double pd_norm;
-//printf("start: %d %d\n",pd_start, pd_stop);
+    local_values.vector[i] = values[2+i];
+//printf("p%d = %g\n",i, local_values.vector[i]);
+  }
+//printf("NUM_VALUES:%d  NUM_PARS:%d  MAX_PD:%d\n", NUM_VALUES, NUM_PARS, MAX_PD);
+//printf("start:%d  stop:%d\n", pd_start, pd_stop);
+
+  double pd_norm = (pd_start == 0 ? 0.0 : result[nq]);
   if (pd_start == 0) {
-    pd_norm = 0.0;
     #ifdef USE_OPENMP
     #pragma omp parallel for
     #endif
     for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0;
-//printf("initializing %d\n", nq);
-  } else {
-    pd_norm = result[nq];
   }
 //printf("start %d %g %g\n", pd_start, pd_norm, result[0]);
 
 #if MAX_PD>0
-  global const double *pd_value = values + NUM_VALUES + 2;
-  global const double *pd_weight = pd_value + details->pd_sum;
+  global const double *pd_value = values + NUM_VALUES;
+  global const double *pd_weight = pd_value + details->num_weights;
 #endif
 
@@ -169,5 +169,5 @@
   global const double *v0 = pd_value + details->pd_offset[0];
   global const double *w0 = pd_weight + details->pd_offset[0];
-//printf("w0:%p, values:%p, diff:%d, %d\n",w0,values,(w0-values),NUM_VALUES);
+//printf("w0:%p, values:%p, diff:%ld, %d\n",w0,values,(w0-values), NUM_VALUES);
 #endif
 
@@ -186,11 +186,10 @@
   int step = pd_start;
 
-
 #if MAX_PD>4
   const double weight5 = 1.0;
   while (i4 < n4) {
-    pvec[p4] = v4[i4];
+    local_values.vector[p4] = v4[i4];
     double weight4 = w4[i4] * weight5;
-//printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 4, p4, i4, n4, pvec[p4], weight4);
+//printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 4, p4, i4, n4, local_values.vector[p4], weight4);
 #elif MAX_PD>3
     const double weight4 = 1.0;
@@ -198,7 +197,7 @@
 #if MAX_PD>3
   while (i3 < n3) {
-    pvec[p3] = v3[i3];
+    local_values.vector[p3] = v3[i3];
     double weight3 = w3[i3] * weight4;
-//printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 3, p3, i3, n3, pvec[p3], weight3);
+//printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 3, p3, i3, n3, local_values.vector[p3], weight3);
 #elif MAX_PD>2
     const double weight3 = 1.0;
@@ -206,7 +205,7 @@
 #if MAX_PD>2
   while (i2 < n2) {
-    pvec[p2] = v2[i2];
+    local_values.vector[p2] = v2[i2];
     double weight2 = w2[i2] * weight3;
-//printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 2, p2, i2, n2, pvec[p2], weight2);
+//printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 2, p2, i2, n2, local_values.vector[p2], weight2);
 #elif MAX_PD>1
     const double weight2 = 1.0;
@@ -214,7 +213,7 @@
 #if MAX_PD>1
   while (i1 < n1) {
-    pvec[p1] = v1[i1];
+    local_values.vector[p1] = v1[i1];
     double weight1 = w1[i1] * weight2;
-//printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 1, p1, i1, n1, pvec[p1], weight1);
+//printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 1, p1, i1, n1, local_values.vector[p1], weight1);
 #elif MAX_PD>0
     const double weight1 = 1.0;
@@ -222,12 +221,12 @@
 #if MAX_PD>0
   if (slow_theta) { // Theta is not in inner loop
-    spherical_correction = fmax(fabs(cos(M_PI_180*pvec[theta_par])), 1.e-6);
+    spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[theta_par])), 1.e-6);
   }
   while(i0 < n0) {
-    pvec[p0] = v0[i0];
+    local_values.vector[p0] = v0[i0];
     double weight0 = w0[i0] * weight1;
-//printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 0, p0, i0, n0, pvec[p0], weight0);
+//printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 0, p0, i0, n0, local_values.vector[p0], weight0);
     if (fast_theta) { // Theta is in inner loop
-      spherical_correction = fmax(fabs(cos(M_PI_180*pvec[p0])), 1.e-6);
+      spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[p0])), 1.e-6);
     }
 #else
@@ -235,9 +234,9 @@
 #endif
 
-//printf("step:%d of %d, pars:",step,pd_stop); for (int i=0; i < NUM_PARS; i++) printf("p%d=%g ",i, pvec[i]); printf("\n");
+//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");
 //printf("sphcor: %g\n", spherical_correction);
 
     #ifdef INVALID
-    if (!INVALID(local_values))
+    if (!INVALID(local_values.table))
     #endif
     {
@@ -248,5 +247,5 @@
         // would be problems looking at models with theta=90.
         const double weight = weight0 * spherical_correction;
-        pd_norm += weight * CALL_VOLUME(local_values);
+        pd_norm += weight * CALL_VOLUME(local_values.table);
 
         #ifdef USE_OPENMP
@@ -263,5 +262,5 @@
           // TODO: what is the magnetic scattering at q=0
           if (qsq > 1.e-16) {
-            double p[4];
+            double p[4];  // dd, du, ud, uu
             p[0] = (qy*cos_mspin + qx*sin_mspin)/qsq;
             p[3] = -p[0];
@@ -278,5 +277,5 @@
                   #define M3 NUM_PARS+13
                   #define SLD(_M_offset, _sld_offset) \
-                      pvec[_sld_offset] = xs * (axis \
+                      local_values.vector[_sld_offset] = xs * (axis \
                       ? (index==1 ? -values[_M_offset+2] : values[_M_offset+2]) \
                       : mag_sld(qx, qy, pk, values[_M_offset], values[_M_offset+1], \
@@ -296,5 +295,5 @@
                   }
                   #endif
-                  scattering += CALL_IQ(q, q_index, local_values);
+                  scattering += CALL_IQ(q, q_index, local_values.table);
                 }
               }
@@ -302,5 +301,5 @@
           }
 #else  // !MAGNETIC
-          const double scattering = CALL_IQ(q, q_index, local_values);
+          const double scattering = CALL_IQ(q, q_index, local_values.table);
 #endif // !MAGNETIC
 //printf("q_index:%d %g %g %g %g\n",q_index, scattering, weight, spherical_correction, weight0);
Index: sasmodels/kernel_iq.cl
===================================================================
--- sasmodels/kernel_iq.cl	(revision 4f1f876c5f6328417afdd2dfe88f52e2f99e444b)
+++ sasmodels/kernel_iq.cl	(revision bde38b559fcdc618ebafdbb69cf7af6120ea8466)
@@ -22,6 +22,6 @@
     int32_t pd_stride[MAX_PD];  // stride to move to the next index at this level
 #endif // MAX_PD > 0
-    int32_t pd_prod;            // total number of voxels in hypercube
-    int32_t pd_sum;             // total length of the weights vector
+    int32_t num_eval;           // total number of voxels in hypercube
+    int32_t num_weights;        // total length of the weights vector
     int32_t num_active;         // number of non-trivial pd loops
     int32_t theta_par;          // id of spherical correction variable
@@ -90,12 +90,6 @@
 
   // Storage for the current parameter values.  These will be updated as we
-  // walk the polydispersity cube.  local_values will be aliased to pvec.
+  // walk the polydispersity cube.
   ParameterBlock local_values;
-
-  // Fill in the initial variables
-  for (int i=0; i < NUM_PARS; i++) {
-    local_values.vector[i] = values[2+i];
-//if (q_index==0) printf("p%d = %g\n",i, local_values.vector[i]);
-  }
 
 #if defined(MAGNETIC) && NUM_MAGNETIC>0
@@ -117,16 +111,21 @@
 #endif // MAGNETIC
 
-  double pd_norm, this_result;
-  if (pd_start == 0) {
-    pd_norm = this_result = 0.0;
-  } else {
-    pd_norm = result[nq];
-    this_result = result[q_index];
-  }
+  // Fill in the initial variables
+  //   values[0] is scale
+  //   values[1] is background
+  for (int i=0; i < NUM_PARS; i++) {
+    local_values.vector[i] = values[2+i];
+//if (q_index==0) printf("p%d = %g\n",i, local_values.vector[i]);
+  }
+//if (q_index==0) printf("NUM_VALUES:%d  NUM_PARS:%d  MAX_PD:%d\n", NUM_VALUES, NUM_PARS, MAX_PD);
+//if (q_index==0) printf("start:%d stop:%d\n", pd_start, pd_stop);
+
+  double pd_norm = (pd_start == 0 ? 0.0 : result[nq]);
+  double this_result = (pd_start == 0 ? 0.0 : result[q_index]);
 //if (q_index==0) printf("start %d %g %g\n", pd_start, pd_norm, this_result);
 
 #if MAX_PD>0
-  global const double *pd_value = values + NUM_VALUES + 2;
-  global const double *pd_weight = pd_value + details->pd_sum;
+  global const double *pd_value = values + NUM_VALUES;
+  global const double *pd_weight = pd_value + details->num_weights;
 #endif
 
@@ -256,5 +255,5 @@
         // TODO: what is the magnetic scattering at q=0
         if (qsq > 1.e-16) {
-          double p[4];  // spin_i, spin_f
+          double p[4];  // dd, du, ud, uu
           p[0] = (qy*cos_mspin + qx*sin_mspin)/qsq;
           p[3] = -p[0];
Index: sasmodels/kernelcl.py
===================================================================
--- sasmodels/kernelcl.py	(revision 40a87fa5096610e24c3ead10299f609e99f093f2)
+++ sasmodels/kernelcl.py	(revision bde38b559fcdc618ebafdbb69cf7af6120ea8466)
@@ -552,10 +552,10 @@
         ]
         #print("Calling OpenCL")
-        #print("values",values)
+        #call_details.show(values)
         # Call kernel and retrieve results
         last_call = None
         step = 100
-        for start in range(0, call_details.pd_prod, step):
-            stop = min(start+step, call_details.pd_prod)
+        for start in range(0, call_details.num_eval, step):
+            stop = min(start + step, call_details.num_eval)
             #print("queuing",start,stop)
             args[1:3] = [np.int32(start), np.int32(stop)]
@@ -563,4 +563,5 @@
                                 None, *args, wait_for=last_call)]
         cl.enqueue_copy(self.queue, self.result, self.result_b)
+        #print("result", self.result)
 
         # Free buffers
@@ -570,5 +571,5 @@
         scale = values[0]/self.result[self.q_input.nq]
         background = values[1]
-        #print("scale",scale,background)
+        #print("scale",scale,values[0],self.result[self.q_input.nq])
         return scale*self.result[:self.q_input.nq] + background
         # return self.result[:self.q_input.nq]
Index: sasmodels/kerneldll.py
===================================================================
--- sasmodels/kerneldll.py	(revision 40a87fa5096610e24c3ead10299f609e99f093f2)
+++ sasmodels/kerneldll.py	(revision bde38b559fcdc618ebafdbb69cf7af6120ea8466)
@@ -380,10 +380,10 @@
             self.real(cutoff), # cutoff
         ]
-        #print("kerneldll values", values)
+        #print("Calling DLL")
+        #call_details.show(values)
         step = 100
-        for start in range(0, call_details.pd_prod, step):
-            stop = min(start+step, call_details.pd_prod)
+        for start in range(0, call_details.num_eval, step):
+            stop = min(start + step, call_details.num_eval)
             args[1:3] = [start, stop]
-            #print("calling DLL")
             kernel(*args) # type: ignore
 
Index: sasmodels/kernelpy.py
===================================================================
--- sasmodels/kernelpy.py	(revision 40a87fa5096610e24c3ead10299f609e99f093f2)
+++ sasmodels/kernelpy.py	(revision bde38b559fcdc618ebafdbb69cf7af6120ea8466)
@@ -100,4 +100,5 @@
     """
     def __init__(self, kernel, model_info, q_input):
+        # type: (callable, ModelInfo, List[np.ndarray]) -> None
         self.dtype = np.dtype('d')
         self.info = model_info
@@ -156,5 +157,9 @@
 
     def __call__(self, call_details, values, cutoff, magnetic):
-        assert isinstance(call_details, details.CallDetails)
+        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray
+        if magnetic:
+            raise NotImplementedError("Magnetism not implemented for pure python models")
+        #print("Calling python kernel")
+        #call_details.show(values)
         res = _loops(self._parameter_vector, self._form, self._volume,
                      self.q_input.nq, call_details, values, cutoff)
@@ -162,4 +167,5 @@
 
     def release(self):
+        # type: () -> None
         """
         Free resources associated with the kernel.
@@ -189,6 +195,6 @@
             return np.ones(nq, 'd')*background
 
-    pd_value = values[2+n_pars:2+n_pars + call_details.pd_sum]
-    pd_weight = values[2+n_pars + call_details.pd_sum:]
+    pd_value = values[2+n_pars:2+n_pars + call_details.num_weights]
+    pd_weight = values[2+n_pars + call_details.num_weights:]
 
     pd_norm = 0.0
@@ -209,5 +215,5 @@
 
     total = np.zeros(nq, 'd')
-    for loop_index in range(call_details.pd_prod):
+    for loop_index in range(call_details.num_eval):
         # update polydispersity parameter values
         if p0_index == p0_length:
Index: sasmodels/sasview_model.py
===================================================================
--- sasmodels/sasview_model.py	(revision fd19811b9fa206bba6e8f8d5edcc44a10d44aa41)
+++ sasmodels/sasview_model.py	(revision bde38b559fcdc618ebafdbb69cf7af6120ea8466)
@@ -24,5 +24,5 @@
 from . import weights
 from . import modelinfo
-from .details import build_details, dispersion_mesh
+from .details import make_kernel_args, dispersion_mesh
 
 try:
@@ -565,5 +565,5 @@
         parameters = self._model_info.parameters
         pairs = [self._get_weights(p) for p in parameters.call_parameters]
-        call_details, values, is_magnetic = build_details(calculator, pairs)
+        call_details, values, is_magnetic = make_kernel_args(calculator, pairs)
         #call_details.show()
         #print("pairs", pairs)
