Index: sasmodels/product.py
===================================================================
--- sasmodels/product.py	(revision 01c8d9ebbcedbff594f32ae8a488e9031320e3ba)
+++ sasmodels/product.py	(revision c036ddbd6dd028e64cda5692721d394c9845adaf)
@@ -37,4 +37,6 @@
 ER_ID = "radius_effective"
 VF_ID = "volfraction"
+BETA_DEFINITION = ("beta_mode", "", 0, [["P*S"],["P*(1+beta*(S-1))"]], "",
+                   "Structure factor dispersion calculation mode")
 
 # TODO: core_shell_sphere model has suppressed the volume ratio calculation
@@ -74,6 +76,6 @@
     translate_name = dict((old.id, new.id) for old, new
                           in zip(s_pars.kernel_parameters[1:], s_list))
-    beta_parameter = Parameter("beta_mode", "", 0, [["P*S"],["P*(1+beta*(S-1))"], "", "Structure factor dispersion calculation mode"])
-    combined_pars = p_pars.kernel_parameters + s_list + [beta_parameter]
+    beta = [Parameter(*BETA_DEFINITION)] if p_info.have_Fq else []
+    combined_pars = p_pars.kernel_parameters + s_list + beta
     parameters = ParameterTable(combined_pars)
     parameters.max_pd = p_pars.max_pd + s_pars.max_pd
@@ -152,5 +154,5 @@
         #: Structure factor modelling interaction between particles.
         self.S = S
-        
+
         #: Model precision. This is not really relevant, since it is the
         #: individual P and S models that control the effective dtype,
@@ -170,5 +172,5 @@
         # in opencl; or both in opencl, but one in single precision and the
         # other in double precision).
-        
+
         p_kernel = self.P.make_kernel(q_vectors)
         s_kernel = self.S.make_kernel(q_vectors)
@@ -196,4 +198,5 @@
         # type: (CallDetails, np.ndarray, float, bool) -> np.ndarray
         p_info, s_info = self.info.composition[1]
+
         # if there are magnetic parameters, they will only be on the
         # form factor P, not the structure factor S.
@@ -207,4 +210,5 @@
         nweights = call_details.num_weights
         weights = values[nvalues:nvalues + 2*nweights]
+
         # Construct the calling parameters for P.
         p_npars = p_info.parameters.npars
@@ -212,4 +216,5 @@
         p_offset = call_details.offset[:p_npars]
         p_details = make_details(p_info, p_length, p_offset, nweights)
+
         # Set p scale to the volume fraction in s, which is the first of the
         # 'S' parameters in the parameter list, or 2+np in 0-origin.
@@ -219,8 +224,10 @@
         p_values.append([0.]*spacer)
         p_values = np.hstack(p_values).astype(self.p_kernel.dtype)
+
         # Call ER and VR for P since these are needed for S.
         p_er, p_vr = calc_er_vr(p_info, p_details, p_values)
         s_vr = (volfrac/p_vr if p_vr != 0. else volfrac)
         #print("volfrac:%g p_er:%g p_vr:%g s_vr:%g"%(volfrac,p_er,p_vr,s_vr))
+
         # Construct the calling parameters for S.
         # The  effective radius is not in the combined parameter list, so
@@ -252,14 +259,47 @@
         s_values.append([0.]*spacer)
         s_values = np.hstack(s_values).astype(self.s_kernel.dtype)
+
         # beta mode is the first parameter after the structure factor pars
         beta_index = 2+p_npars+s_npars
         beta_mode = values[beta_index]
+
         # Call the kernels
-        if beta_mode: # beta:
+        s_result = self.s_kernel.Iq(s_details, s_values, cutoff, False)
+        scale, background = values[0], values[1]
+        if beta_mode:
             F1, F2, volume_avg = self.p_kernel.beta(p_details, p_values, cutoff, magnetic)
+            combined_scale = scale*volfrac/volume_avg
+            # Define lazy results based on intermediate values.
+            # The return value for the calculation should be an ordered
+            # dictionary containing any result the user might want to see
+            # at the end of the calculation, including scalars, strings,
+            # and plottable data.  Don't want to build this structure during
+            # fits, only when displaying the final result (or a one-off
+            # computation which asks for it).
+            # Do not use the current hack of storing the intermediate values
+            # in self.results since that leads to awkward threading issues.
+            # Instead return the function along with the bundle of inputs.
+            # P and Q may themselves have intermediate results they want to
+            # include, such as A and B if P = A + B.  Might use this mechanism
+            # to return the computed effective radius as well.
+            #def lazy_results(Q, S, F1, F2, scale):
+            #    """
+            #    beta = F1**2 / F2  # what about divide by zero errors?
+            #    return {
+            #        'P' : Data1D(Q, scale*F2),
+            #        'beta': Data1D(Q, beta),
+            #        'S' : Data1D(Q, S),
+            #        'Seff': Data1D(Q, 1 + beta*(S-1)),
+            #        'I' : Data1D(Q, scale*(F2 + (F1**2)*(S-1)) + background),
+            #    }
+            #lazy_pars = s_result, F1, F2, combined_scale
+            self.results = [F2, s_result]
+            final_result = combined_scale*(F2 + (F1**2)*(s_result - 1)) + background
         else:
             p_result = self.p_kernel.Iq(p_details, p_values, cutoff, magnetic)
-        s_result = self.s_kernel.Iq(s_details, s_values, cutoff, False)
-        #print("p_npars",p_npars,s_npars,p_er,s_vr,values[2+p_npars+1:2+p_npars+s_npars])
+            # remember the parts for plotting later
+            self.results = [p_result, s_result]
+            final_result = scale*(p_result*s_result) + background
+
         #call_details.show(values)
         #print("values", values)
@@ -272,14 +312,4 @@
         #plt.subplot(212); plt.loglog(self.s_kernel.q_input.q, s_result, '-')
         #plt.figure()
-        if beta_mode:#beta
-            beta_factor = F1**2/F2
-            Sq_eff = 1+beta_factor*(s_result - 1)
-            self.results = [F2, Sq_eff,F1,s_result]
-            final_result = volfrac*values[0]*(F2 + (F1**2)*(s_result - 1))/volume_avg+values[1]
-            #final_result =  volfrac * values[0] * F2 * Sq_eff / volume_avg + values[1]
-        else:
-            # remember the parts for plotting later
-            self.results = [p_result, s_result]
-            final_result = values[0]*(p_result*s_result) + values[1]
         return final_result
 
