Index: sasmodels/product.py
===================================================================
--- sasmodels/product.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
+++ sasmodels/product.py	(revision 01c8d9ebbcedbff594f32ae8a488e9031320e3ba)
@@ -16,5 +16,5 @@
 import numpy as np  # type: ignore
 
-from .modelinfo import ParameterTable, ModelInfo
+from .modelinfo import ParameterTable, ModelInfo, Parameter
 from .kernel import KernelModel, Kernel
 from .details import make_details, dispersion_mesh
@@ -74,5 +74,6 @@
     translate_name = dict((old.id, new.id) for old, new
                           in zip(s_pars.kernel_parameters[1:], s_list))
-    combined_pars = p_pars.kernel_parameters + 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]
     parameters = ParameterTable(combined_pars)
     parameters.max_pd = p_pars.max_pd + s_pars.max_pd
@@ -151,4 +152,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,
@@ -168,4 +170,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)
@@ -193,5 +196,4 @@
         # 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.
@@ -205,5 +207,4 @@
         nweights = call_details.num_weights
         weights = values[nvalues:nvalues + 2*nweights]
-
         # Construct the calling parameters for P.
         p_npars = p_info.parameters.npars
@@ -218,10 +219,8 @@
         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
@@ -253,9 +252,13 @@
         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
-        p_result = self.p_kernel(p_details, p_values, cutoff, magnetic)
-        s_result = self.s_kernel(s_details, s_values, cutoff, False)
-
+        if beta_mode: # beta:
+            F1, F2, volume_avg = self.p_kernel.beta(p_details, p_values, cutoff, magnetic)
+        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])
         #call_details.show(values)
@@ -265,14 +268,19 @@
         #s_details.show(s_values)
         #print("=>", s_result)
-
-        # remember the parts for plotting later
-        self.results = [p_result, s_result]
-
         #import pylab as plt
         #plt.subplot(211); plt.loglog(self.p_kernel.q_input.q, p_result, '-')
         #plt.subplot(212); plt.loglog(self.s_kernel.q_input.q, s_result, '-')
         #plt.figure()
-
-        return values[0]*(p_result*s_result) + values[1]
+        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
 
     def release(self):
