Index: sasmodels/core.py
===================================================================
--- sasmodels/core.py	(revision 01c8d9ebbcedbff594f32ae8a488e9031320e3ba)
+++ sasmodels/core.py	(revision c036ddbd6dd028e64cda5692721d394c9845adaf)
@@ -140,5 +140,4 @@
     used with functions in generate to build the docs or extract model info.
     """
-
     if "+" in model_string:
         parts = [load_model_info(part)
@@ -290,4 +289,5 @@
     else:
         numpy_dtype = np.dtype(dtype)
+
     # Make sure that the type is supported by opencl, otherwise use dll
     if platform == "ocl":
Index: sasmodels/data.py
===================================================================
--- sasmodels/data.py	(revision 01c8d9ebbcedbff594f32ae8a488e9031320e3ba)
+++ sasmodels/data.py	(revision c036ddbd6dd028e64cda5692721d394c9845adaf)
@@ -522,5 +522,4 @@
                   else 'log')
         plt.xscale(xscale)
-        
         plt.xlabel("$q$/A$^{-1}$")
         plt.yscale(yscale)
Index: sasmodels/generate.py
===================================================================
--- sasmodels/generate.py	(revision 01c8d9ebbcedbff594f32ae8a488e9031320e3ba)
+++ sasmodels/generate.py	(revision c036ddbd6dd028e64cda5692721d394c9845adaf)
@@ -671,5 +671,4 @@
 
 
-# type in IQXY pattern could be single, float, double, long double, ...
 _IQXY_PATTERN = re.compile(r"(^|\s)double\s+I(?P<mode>q(ac|abc|xy))\s*[(]",
                            flags=re.MULTILINE)
@@ -701,7 +700,13 @@
     return 'qa'
 
-# type in IQXY pattern could be single, float, double, long double, ...
+# Note: not presently used.  Need to know whether Fq is available before
+# trying to compile the source.  Leave the code here in case we decide that
+# define have_Fq for each form factor is too tedious and error prone.
 _FQ_PATTERN = re.compile(r"(^|\s)void\s+Fq[(]", flags=re.MULTILINE)
-def has_Fq(source):
+def contains_Fq(source):
+    # type: (List[str]) -> bool
+    """
+    Return True if C source defines "void Fq(".
+    """
     for code in source:
         m = _FQ_PATTERN.search(code)
@@ -739,8 +744,10 @@
     # for computing volume even if we allow non-disperse volume parameters.
     partable = model_info.parameters
+
     # Load templates and user code
     kernel_header = load_template('kernel_header.c')
     kernel_code = load_template('kernel_iq.c')
     user_code = [(f, open(f).read()) for f in model_sources(model_info)]
+
     # Build initial sources
     source = []
@@ -755,4 +762,5 @@
     q, qx, qy, qab, qa, qb, qc \
         = [Parameter(name=v) for v in 'q qx qy qab qa qb qc'.split()]
+
     # Generate form_volume function, etc. from body only
     if isinstance(model_info.form_volume, str):
@@ -804,12 +812,13 @@
     source.append(call_volume)
     model_refs = _call_pars("_v.", partable.iq_parameters)
-    #create varaible BETA to turn on and off beta approximation
-    BETA = has_Fq(source)
-    if not BETA:
+
+    if model_info.have_Fq:
+        pars = ",".join(["_q", "&_F1", "&_F2",] + model_refs)
+        call_iq = "#define CALL_FQ(_q, _F1, _F2, _v) Fq(%s)" % pars
+        clear_iq = "#undef CALL_FQ"
+    else:
         pars = ",".join(["_q"] + model_refs)
         call_iq = "#define CALL_IQ(_q, _v) Iq(%s)" % pars
-    else:
-        pars = ",".join(["_q", "&_F1", "&_F2",] + model_refs)
-        call_iq = "#define CALL_IQ(_q, _F1, _F2, _v) Fq(%s)" % pars
+        clear_iq = "#undef CALL_IQ"
     if xy_mode == 'qabc':
         pars = ",".join(["_qa", "_qb", "_qc"] + model_refs)
@@ -820,8 +829,16 @@
         call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqac(%s)" % pars
         clear_iqxy = "#undef CALL_IQ_AC"
-    elif xy_mode == 'qa':
+    elif xy_mode == 'qa' and not model_info.have_Fq:
         pars = ",".join(["_qa"] + model_refs)
         call_iqxy = "#define CALL_IQ_A(_qa,_v) Iq(%s)" % pars
         clear_iqxy = "#undef CALL_IQ_A"
+    elif xy_mode == 'qa' and model_info.have_Fq:
+        pars = ",".join(["_qa", "&_F1", "&_F2",] + model_refs)
+        # Note: uses rare C construction (expr1, expr2) which computes
+        # expr1 then expr2 and evaluates to expr2.  This allows us to
+        # leave it looking like a function even though it is returning
+        # its values by reference.
+        call_iqxy = "#define CALL_FQ_A(_qa,_F1,_F2,_v) (Fq(%s),_F2)" % pars
+        clear_iqxy = "#undef CALL_FQ_A"
     elif xy_mode == 'qxy':
         orientation_refs = _call_pars("_v.", partable.orientation_parameters)
@@ -840,5 +857,4 @@
                if p.type == 'sld']
     # Fill in definitions for numbers of parameters
-    source.append("#define BETA %d" %(1 if BETA else 0))
     source.append("#define MAX_PD %s"%partable.max_pd)
     source.append("#define NUM_PARS %d"%partable.npars)
@@ -848,6 +864,6 @@
     source.append("#define PROJECTION %d"%PROJECTION)
     # TODO: allow mixed python/opencl kernels?
-    ocl = _kernels(kernel_code, call_iq, call_iqxy, clear_iqxy, model_info.name)
-    dll = _kernels(kernel_code, call_iq, call_iqxy, clear_iqxy, model_info.name)
+    ocl = _kernels(kernel_code, call_iq, clear_iq, call_iqxy, clear_iqxy, model_info.name)
+    dll = _kernels(kernel_code, call_iq, clear_iq, call_iqxy, clear_iqxy, model_info.name)
 
     result = {
@@ -855,9 +871,8 @@
         'opencl': '\n'.join(source+ocl[0]+ocl[1]+ocl[2]),
     }
-    #print(result['dll'])
     return result
 
 
-def _kernels(kernel, call_iq, call_iqxy, clear_iqxy, name):
+def _kernels(kernel, call_iq, clear_iq, call_iqxy, clear_iqxy, name):
     # type: ([str,str], str, str, str) -> List[str]
     code = kernel[0]
@@ -869,5 +884,5 @@
         '#line 1 "%s Iq"' % path,
         code,
-        "#undef CALL_IQ",
+        clear_iq,
         "#undef KERNEL_NAME",
         ]
@@ -1075,5 +1090,5 @@
         model_info = make_model_info(kernel_module)
         source = make_source(model_info)
-        #print(source['dll'])
+        print(source['dll'])
 
 
Index: sasmodels/kernel_iq.c
===================================================================
--- sasmodels/kernel_iq.c	(revision 01c8d9ebbcedbff594f32ae8a488e9031320e3ba)
+++ sasmodels/kernel_iq.c	(revision c036ddbd6dd028e64cda5692721d394c9845adaf)
@@ -29,4 +29,6 @@
 //  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
@@ -85,5 +87,5 @@
   out_spin = clip(out_spin, 0.0, 1.0);
   // Previous version of this function took the square root of the weights,
-  // under the assumption that 
+  // under the assumption that
   //
   //     w*I(q, rho1, rho2, ...) = I(q, sqrt(w)*rho1, sqrt(w)*rho2, ...)
@@ -272,4 +274,5 @@
 
 // ==================== KERNEL CODE ========================
+#define COMPUTE_F1_F2 defined(CALL_FQ)
 kernel
 void KERNEL_NAME(
@@ -338,30 +341,34 @@
   // The code differs slightly between opencl and dll since opencl is only
   // seeing one q value (stored in the variable "this_result") while the dll
-  // version must loop over all q
-  #if BETA
-  double *beta_result = &(result[nq+1]); // = result + nq + 1
-  double weight_norm = (pd_start == 0 ? 0.0 : result[2*nq+1]);
-  #endif
+  // version must loop over all q.
   #ifdef USE_OPENCL
-    double pd_norm = (pd_start == 0 ? 0.0 : result[nq]);
-    double this_result = (pd_start == 0 ? 0.0 : result[q_index]);
-    #if BETA
-      double this_beta_result = (pd_start == 0 ? 0.0 : result[nq + q_index];
+    #if COMPUTE_F1_F2
+      double pd_norm = (pd_start == 0 ? 0.0 : result[nq]);
+      double weight_norm = (pd_start == 0 ? 0.0 : result[2*nq+1]);
+      double this_F2 = (pd_start == 0 ? 0.0 : result[q_index]);
+      double this_F1 = (pd_start == 0 ? 0.0 : result[q_index+nq+1]);
+    #else
+      double pd_norm = (pd_start == 0 ? 0.0 : result[nq]);
+      double this_result = (pd_start == 0 ? 0.0 : result[q_index]);
+    #endif
   #else // !USE_OPENCL
-    double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
+    #if COMPUTE_F1_F2
+      double pd_norm = (pd_start == 0 ? 0.0 : result[nq]);
+      double weight_norm = (pd_start == 0 ? 0.0 : result[2*nq+1]);
+    #else
+      double pd_norm = (pd_start == 0 ? 0.0 : result[nq]);
+    #endif
     if (pd_start == 0) {
       #ifdef USE_OPENMP
       #pragma omp parallel for
       #endif
-      for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0;
-      #if BETA
-      for (int q_index=0; q_index < nq; q_index++) beta_result[q_index] = 0.0;
+      #if COMPUTE_F1_F2
+          for (int q_index=0; q_index < 2*nq+2; q_index++) result[q_index] = 0.0;
+      #else
+          for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0;
       #endif
     }
     //if (q_index==0) printf("start %d %g %g\n", pd_start, pd_norm, result[0]);
 #endif // !USE_OPENCL
-
-
-
 
 
@@ -454,21 +461,29 @@
 
 
-#if defined(CALL_IQ)
-  // unoriented 1D
+#if defined(CALL_FQ) // COMPUTE_F1_F2 is true
+  // unoriented 1D returning <F> and <F^2>
   double qk;
-  #if BETA == 0
-    #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)
-
-  // unoriented 1D Beta
-  #elif BETA == 1
-    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_IQ(qk,F1,F2,local_values.table)
-  #endif
+  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>
+  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)
 
 #elif defined(CALL_IQ_A)
@@ -504,4 +519,5 @@
   #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
@@ -517,5 +533,6 @@
 // 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)  // no orientation
+#if defined(CALL_IQ) || defined(CALL_IQ_A) || defined(CALL_FQ) || defined(CALL_FQ_A)
+  // no orientation
   #define APPLY_PROJECTION() const double weight=weight0
 #elif defined(CALL_IQ_XY) // pass orientation to the model
@@ -668,9 +685,9 @@
     if (weight > cutoff) {
       pd_norm += weight * CALL_VOLUME(local_values.table);
+      #if COMPUTE_F1_F2
+      weight_norm += weight;
+      #endif
       BUILD_ROTATION();
-#if BETA
-    if (weight > cutoff) {
-      weight_norm += weight;}
-#endif
+
 #ifndef USE_OPENCL
       // DLL needs to explicitly loop over the q values.
@@ -718,8 +735,6 @@
           }
         #else  // !MAGNETIC
-          #if defined(CALL_IQ) && BETA == 1
-            CALL_KERNEL();
-            const double scatteringF1 = F1;
-            const double scatteringF2 = F2;
+          #if COMPUTE_F1_F2
+            CALL_KERNEL(); // sets F1 and F2 by reference
           #else
             const double scattering = CALL_KERNEL();
@@ -729,16 +744,15 @@
 
         #ifdef USE_OPENCL
-          #if defined(CALL_IQ)&& BETA == 1
-             this_result += weight * scatteringF2;
-             this_beta_result += weight * scatteringF1;
-            #else
-              this_result += weight * scattering;
+          #if COMPUTE_F1_F2
+            this_F1 += weight * F1;
+            this_F2 += weight * F2;
+          #else
+            this_result += weight * scattering;
           #endif
         #else // !USE_OPENCL
-          #if defined(CALL_IQ)&& BETA == 1
-            result[q_index] += weight * scatteringF2;
-            beta_result[q_index] += weight * scatteringF1;
-            #endif
-            #else
+          #if COMPUTE_F1_F2
+            result[q_index] += weight * F2;
+            result[q_index+nq+1] += weight * F1;
+          #else
             result[q_index] += weight * scattering;
           #endif
@@ -767,16 +781,23 @@
 // Remember the current result and the updated norm.
 #ifdef USE_OPENCL
-  result[q_index] = this_result;
-  if (q_index == 0) result[nq] = pd_norm;
-  #if BETA
-  beta_result[q_index] = this_beta_result;
+  #if COMPUTE_F1_F2
+    result[q_index] = this_F2;
+    result[q_index+nq+1] = this_F1;
+    if (q_index == 0) {
+      result[nq] = pd_norm;
+      result[2*nq+1] = weight_norm;
+    }
+  #else
+    result[q_index] = this_result;
+    if (q_index == 0) result[nq] = pd_norm;
   #endif
-  if (q_index == 0) result[nq] = pd_norm;
 
 //if (q_index == 0) printf("res: %g/%g\n", result[0], pd_norm);
 #else // !USE_OPENCL
-  result[nq] = pd_norm;
-  #if BETA
-  result[2*nq+1] = weight_norm;
+  #if COMPUTE_F1_F2
+    result[nq] = pd_norm;
+    result[2*nq+1] = weight_norm;
+  #else
+    result[nq] = pd_norm;
   #endif
 //printf("res: %g/%g\n", result[0], pd_norm);
@@ -784,4 +805,5 @@
 
 // ** clear the macros in preparation for the next kernel **
+#undef COMPUTE_F1_F2
 #undef PD_INIT
 #undef PD_OPEN
Index: sasmodels/kernelcl.py
===================================================================
--- sasmodels/kernelcl.py	(revision 01c8d9ebbcedbff594f32ae8a488e9031320e3ba)
+++ sasmodels/kernelcl.py	(revision c036ddbd6dd028e64cda5692721d394c9845adaf)
@@ -420,6 +420,4 @@
         if self.program is None:
             compile_program = environment().compile_program
-            with open('model.c','w') as fid:
-                print(self.source['opencl'], file=fid)
             timestamp = generate.ocl_timestamp(self.info)
             self.program = compile_program(
@@ -540,6 +538,8 @@
         self.dtype = dtype
         self.dim = '2d' if q_input.is_2d else '1d'
-        # plus three for the normalization values
-        self.result = np.empty(2*q_input.nq+2,dtype)
+        # leave room for f1/f2 results in case we need to compute beta for 1d models
+        num_returns = 1 if self.dim == '2d' else 2  #
+        # plus 1 for the normalization value
+        self.result = np.empty((q_input.nq+1)*num_returns, dtype)
 
         # Inputs and outputs for each kernel call
@@ -549,5 +549,5 @@
 
         self.result_b = cl.Buffer(self.queue.context, mf.READ_WRITE,
-                                  q_input.global_size[0] * dtype.itemsize)
+                                  q_input.global_size[0] * num_returns * dtype.itemsize)
         self.q_input = q_input # allocated by GpuInput above
 
@@ -557,7 +557,31 @@
                      else np.float16 if dtype == generate.F16
                      else np.float32)  # will never get here, so use np.float32
-    __call__= Iq
 
     def Iq(self, call_details, values, cutoff, magnetic):
+        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray
+        self._call_kernel(call_details, values, cutoff, magnetic)
+        #print("returned",self.q_input.q, self.result)
+        pd_norm = self.result[self.q_input.nq]
+        scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0)
+        background = values[1]
+        #print("scale",scale,background)
+        return scale*self.result[:self.q_input.nq] + background
+    __call__ = Iq
+
+    def beta(self, call_details, values, cutoff, magnetic):
+        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray
+        if self.dim == '2d':
+            raise NotImplementedError("beta not yet supported for 2D")
+        self._call_kernel(call_details, values, cutoff, magnetic)
+        w_norm = self.result[2*self.q_input.nq + 1]
+        pd_norm = self.result[self.q_input.nq]
+        if w_norm == 0.:
+            w_norm = 1.
+        F2 = self.result[:self.q_input.nq]/w_norm
+        F1 = self.result[self.q_input.nq+1:2*self.q_input.nq+1]/w_norm
+        volume_avg = pd_norm/w_norm
+        return F1, F2, volume_avg
+
+    def _call_kernel(self, call_details, values, cutoff, magnetic):
         # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray
         context = self.queue.context
@@ -575,5 +599,5 @@
         ]
         #print("Calling OpenCL")
-        call_details.show(values)
+        #call_details.show(values)
         #Call kernel and retrieve results
         wait_for = None
@@ -601,14 +625,4 @@
                 v.release()
 
-        pd_norm = self.result[self.q_input.nq]
-        scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0)
-        background = values[1]
-        #print("scale",scale,values[0],self.result[self.q_input.nq],background)
-        return scale*self.result[:self.q_input.nq] + background
-        # return self.result[:self.q_input.nq]
-     #NEEDS TO BE FINISHED FOR OPENCL
-     def beta():
-         return 0
-
     def release(self):
         # type: () -> None
Index: sasmodels/kerneldll.py
===================================================================
--- sasmodels/kerneldll.py	(revision 01c8d9ebbcedbff594f32ae8a488e9031320e3ba)
+++ sasmodels/kerneldll.py	(revision c036ddbd6dd028e64cda5692721d394c9845adaf)
@@ -258,6 +258,6 @@
         # Note: if there is a syntax error then compile raises an error
         # and the source file will not be deleted.
-        #os.unlink(filename)
-        print("saving compiled file in %r"%filename)
+        os.unlink(filename)
+        #print("saving compiled file in %r"%filename)
     return dll
 
@@ -377,5 +377,8 @@
         self.dtype = q_input.dtype
         self.dim = '2d' if q_input.is_2d else '1d'
-        self.result = np.empty(2*q_input.nq+2, q_input.dtype)
+        # leave room for f1/f2 results in case we need to compute beta for 1d models
+        num_returns = 1 if self.dim == '2d' else 2  #
+        # plus 1 for the normalization value
+        self.result = np.empty((q_input.nq+1)*num_returns, self.dtype)
         self.real = (np.float32 if self.q_input.dtype == generate.F32
                      else np.float64 if self.q_input.dtype == generate.F64
@@ -417,7 +420,5 @@
             self.real(cutoff), # cutoff
         ]
-        #print(self.beta_result.ctypes.data)
-#        print("Calling DLL line 397")
-#        print("values", values)
+        #print("Calling DLL")
         #call_details.show(values)
         step = 100
Index: sasmodels/modelinfo.py
===================================================================
--- sasmodels/modelinfo.py	(revision 01c8d9ebbcedbff594f32ae8a488e9031320e3ba)
+++ sasmodels/modelinfo.py	(revision c036ddbd6dd028e64cda5692721d394c9845adaf)
@@ -417,4 +417,5 @@
     # scale and background are implicit parameters
     COMMON = [Parameter(*p) for p in COMMON_PARAMETERS]
+
     def __init__(self, parameters):
         # type: (List[Parameter]) -> None
@@ -440,4 +441,5 @@
         self.form_volume_parameters = [p for p in self.kernel_parameters
                                        if p.type == 'volume']
+
         # Theta offset
         offset = 0
@@ -786,4 +788,6 @@
     info.category = getattr(kernel_module, 'category', None)
     info.structure_factor = getattr(kernel_module, 'structure_factor', False)
+    # TODO: find Fq by inspection
+    info.have_Fq = getattr(kernel_module, 'have_Fq', False)
     info.profile_axes = getattr(kernel_module, 'profile_axes', ['x', 'y'])
     info.source = getattr(kernel_module, 'source', [])
@@ -909,4 +913,7 @@
     #: provided in the file.
     structure_factor = None # type: bool
+    #: True if the model defines an Fq function with signature
+    #: void Fq(double q, double *F1, double *F2, ...)
+    have_Fq = False
     #: List of C source files used to define the model.  The source files
     #: should define the *Iq* function, and possibly *Iqac* or *Iqabc* if the
Index: sasmodels/models/core_shell_sphere.py
===================================================================
--- sasmodels/models/core_shell_sphere.py	(revision dc76240c1dc8dcc0235d0432f3a909bea6e0c6a2)
+++ sasmodels/models/core_shell_sphere.py	(revision c036ddbd6dd028e64cda5692721d394c9845adaf)
@@ -58,6 +58,6 @@
 title = "Form factor for a monodisperse spherical particle with particle with a core-shell structure."
 description = """
-    F^2(q) = 3/V_s [V_c (sld_core-sld_shell) (sin(q*radius)-q*radius*cos(q*radius))/(q*radius)^3
-                   + V_s (sld_shell-sld_solvent) (sin(q*r_s)-q*r_s*cos(q*r_s))/(q*r_s)^3]
+    F(q) = [V_c (sld_core-sld_shell) 3 (sin(q*radius)-q*radius*cos(q*radius))/(q*radius)^3
+            + V_s (sld_shell-sld_solvent) 3 (sin(q*r_s)-q*r_s*cos(q*r_s))/(q*r_s)^3]
 
             V_s: Volume of the sphere shell
Index: sasmodels/models/ellipsoid.c
===================================================================
--- sasmodels/models/ellipsoid.c	(revision 01c8d9ebbcedbff594f32ae8a488e9031320e3ba)
+++ sasmodels/models/ellipsoid.c	(revision c036ddbd6dd028e64cda5692721d394c9845adaf)
@@ -5,4 +5,5 @@
 }
 
+/* Fq overrides Iq
 static  double
 Iq(double q,
@@ -20,5 +21,5 @@
     //     i(h) = int_0^1 Phi^2(h a sqrt(1 + u^2(v^2-1)) du
     const double v_square_minus_one = square(radius_polar/radius_equatorial) - 1.0;
-  
+
     // translate a point in [-1,1] to a point in [0, 1]
     // const double u = GAUSS_Z[i]*(upper-lower)/2 + (upper+lower)/2;
@@ -36,5 +37,6 @@
     const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial);
     return 1.0e-4 * s * s * form;
-} 
+}
+*/
 
 static void
@@ -72,12 +74,7 @@
     const double form_avg = total_F1*zm;
     const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial);
+    *F1 = 1e-2 * s * form_avg;
     *F2 = 1e-4 * s * s * form_squared_avg;
-    *F1 = 1e-2 * s * form_avg;
 }
-
-
-
-
-
 
 static double
Index: sasmodels/models/ellipsoid.py
===================================================================
--- sasmodels/models/ellipsoid.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
+++ sasmodels/models/ellipsoid.py	(revision c036ddbd6dd028e64cda5692721d394c9845adaf)
@@ -161,5 +161,8 @@
              ]
 
+
 source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "ellipsoid.c"]
+
+have_Fq = True
 
 def ER(radius_polar, radius_equatorial):
Index: sasmodels/models/sphere.c
===================================================================
--- sasmodels/models/sphere.c	(revision 01c8d9ebbcedbff594f32ae8a488e9031320e3ba)
+++ sasmodels/models/sphere.c	(revision c036ddbd6dd028e64cda5692721d394c9845adaf)
@@ -4,4 +4,5 @@
 }
 
+// TODO: remove Iq when Iqxy can use Fq directly
 static double Iq(double q, double sld, double sld_solvent, double radius)
 {
@@ -11,8 +12,8 @@
 static void Fq(double q, double *F1,double *F2, double sld, double solvent_sld, double radius)
 {
-    const double fq = sphere_volume(radius) * sas_3j1x_x(q*radius);
+    const double fq = sas_3j1x_x(q*radius);
     const double contrast = (sld - solvent_sld);
-    const double form = 1e-2*contrast * fq;
-    *F1= form;
+    const double form = 1e-2 * contrast * sphere_volume(radius) * fq;
+    *F1 = form;
     *F2 = form*form;
 }
Index: sasmodels/models/sphere.py
===================================================================
--- sasmodels/models/sphere.py	(revision 01c8d9ebbcedbff594f32ae8a488e9031320e3ba)
+++ sasmodels/models/sphere.py	(revision c036ddbd6dd028e64cda5692721d394c9845adaf)
@@ -67,6 +67,24 @@
              ]
 
-source = ["lib/sas_3j1x_x.c", "lib/sphere_form.c", "sphere.c"]
+source = ["lib/sas_3j1x_x.c", "lib/sphere_form.c"]
 
+c_code = """
+static double form_volume(double radius)
+{
+    return sphere_volume(radius);
+}
+
+static void Fq(double q, double *F1,double *F2, double sld, double solvent_sld, double radius)
+{
+    const double fq = sas_3j1x_x(q*radius);
+    const double contrast = (sld - solvent_sld);
+    const double form = 1e-2 * contrast * sphere_volume(radius) * fq;
+    *F1 = form;
+    *F2 = form*form;
+}
+"""
+
+# TODO: figure this out by inspection
+have_Fq = True
 
 def ER(radius):
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
 
