Index: sasmodels/generate.py
===================================================================
--- sasmodels/generate.py	(revision c036ddbd6dd028e64cda5692721d394c9845adaf)
+++ sasmodels/generate.py	(revision 6e7ba14b7bd7549d79f388fe396c80945675b80b)
@@ -805,4 +805,5 @@
         refs = _call_pars("_v.", partable.form_volume_parameters)
         call_volume = "#define CALL_VOLUME(_v) form_volume(%s)"%(",".join(refs))
+        call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(mode, _v) effective_radius(mode, %s)"%(",".join(refs))
     else:
         # Model doesn't have volume.  We could make the kernel run a little
@@ -810,5 +811,7 @@
         # the ifdef's reduce readability more than is worthwhile.
         call_volume = "#define CALL_VOLUME(v) 1.0"
+        call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(mode, v) 0.0"
     source.append(call_volume)
+    source.append(call_effective_radius)
     model_refs = _call_pars("_v.", partable.iq_parameters)
 
Index: sasmodels/kernel.py
===================================================================
--- sasmodels/kernel.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
+++ sasmodels/kernel.py	(revision 6e7ba14b7bd7549d79f388fe396c80945675b80b)
@@ -41,9 +41,54 @@
     dtype = None  # type: np.dtype
 
-    def __call__(self, call_details, values, cutoff, magnetic):
+    def Iq(self, call_details, values, cutoff, magnetic):
         # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray
-        raise NotImplementedError("need to implement __call__")
+        Pq, Reff = self.Pq_Reff(call_details, values, cutoff, magnetic, effective_radius_type=0)
+        scale = values[0]
+        background = values[1]
+        return scale*Pq + background
+    __call__ = Iq
+
+    def Pq_Reff(self, call_details, values, cutoff, magnetic, effective_radius_type):
+        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray
+        self._call_kernel(call_details, values, cutoff, magnetic, effective_radius_type)
+        #print("returned",self.q_input.q, self.result)
+        nout = 2 if self.info.have_Fq else 1
+        total_weight = self.result[nout*self.q_input.nq + 0]
+        if total_weight == 0.:
+            total_weight = 1.
+        weighted_volume = self.result[nout*self.q_input.nq + 1]
+        weighted_radius = self.result[nout*self.q_input.nq + 2]
+        effective_radius = weighted_radius/total_weight
+        # compute I = scale*P + background
+        #           = scale*(sum(w*F^2)/sum w)/(sum (w*V)/sum w) + background
+        #           = scale/sum (w*V) * sum(w*F^2) + background
+        scale = values[0]/(weighted_volume if weighted_volume != 0.0 else 1.0)
+        background = values[1]
+        #print("scale",scale,background)
+        return self.result[0:nout*self.q_input.nq:nout], effective_radius
+
+    def beta(self, call_details, values, cutoff, magnetic, effective_radius_type):
+        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray
+        if self.dim == '2d':
+            raise NotImplementedError("beta not yet supported for 2D")
+        if not self.info.have_Fq:
+            raise NotImplementedError("beta not yet supported for "+self.info.id)
+        self._call_kernel(call_details, values, cutoff, magnetic, effective_radius_type)
+        total_weight = self.result[2*self.q_input.nq + 0]
+        if total_weight == 0.:
+            total_weight = 1.
+        weighted_volume = self.result[2*self.q_input.nq + 1]
+        weighted_radius = self.result[2*self.q_input.nq + 2]
+        volume_average = weighted_volume/total_weight
+        effective_radius = weighted_radius/total_weight
+        F2 = self.result[0:2*self.q_input.nq:2]/total_weight
+        F1 = self.result[1:2*self.q_input.nq:2]/total_weight
+        return F1, F2, volume_average, effective_radius
 
     def release(self):
         # type: () -> None
         pass
+
+    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type):
+        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray
+        raise NotImpmentedError()
Index: sasmodels/kernel_iq.c
===================================================================
--- sasmodels/kernel_iq.c	(revision c036ddbd6dd028e64cda5692721d394c9845adaf)
+++ sasmodels/kernel_iq.c	(revision 6e7ba14b7bd7549d79f388fe396c80945675b80b)
@@ -27,4 +27,5 @@
 //      parameters in the parameter table.
 //  CALL_VOLUME(table) : call the form volume function
+//  CALL_EFFECTIVE_RADIUS(type, table) : call the R_eff function
 //  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.
@@ -284,5 +285,6 @@
     global const double *q, // nq q values, with padding to boundary
     global double *result,  // nq+1 return values, again with padding
-    const double cutoff     // cutoff in the dispersity weight product
+    const double cutoff,     // cutoff in the dispersity weight product
+    int effective_radius_type // which effective radius to compute
     )
 {
@@ -344,18 +346,24 @@
   #ifdef USE_OPENCL
     #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]);
+      double weight_norm = (pd_start == 0 ? 0.0 : result[2*nq]);
+      double weighted_volume = (pd_start == 0 ? 0.0 : result[2*nq+1]);
+      double weighted_radius = (pd_start == 0 ? 0.0 : result[2*nq+2]);
+      double this_F2 = (pd_start == 0 ? 0.0 : result[2*q_index+0]);
+      double this_F1 = (pd_start == 0 ? 0.0 : result[2*q_index+1]);
     #else
-      double pd_norm = (pd_start == 0 ? 0.0 : result[nq]);
+      double weight_norm = (pd_start == 0 ? 0.0 : result[nq]);
+      double weighted_volume = (pd_start == 0 ? 0.0 : result[nq+1]);
+      double weighted_radius = (pd_start == 0 ? 0.0 : result[nq+2]);
       double this_result = (pd_start == 0 ? 0.0 : result[q_index]);
     #endif
   #else // !USE_OPENCL
     #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 weight_norm = (pd_start == 0 ? 0.0 : result[2*nq]);
+      double weighted_volume = (pd_start == 0 ? 0.0 : result[2*nq+1]);
+      double weighted_radius = (pd_start == 0 ? 0.0 : result[2*nq+2]);
     #else
-      double pd_norm = (pd_start == 0 ? 0.0 : result[nq]);
+      double weight_norm = (pd_start == 0 ? 0.0 : result[nq]);
+      double weighted_volume = (pd_start == 0 ? 0.0 : result[nq+1]);
+      double weighted_radius = (pd_start == 0 ? 0.0 : result[nq+2]);
     #endif
     if (pd_start == 0) {
@@ -364,10 +372,11 @@
       #endif
       #if COMPUTE_F1_F2
-          for (int q_index=0; q_index < 2*nq+2; q_index++) result[q_index] = 0.0;
+          // 2*nq for F^2,F pairs
+          for (int q_index=0; q_index < 2*nq; 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]);
+    //if (q_index==0) printf("start %d %g %g\n", pd_start, weighted_volume, result[0]);
 #endif // !USE_OPENCL
 
@@ -684,8 +693,11 @@
     // Note: weight==0 must always be excluded
     if (weight > cutoff) {
-      pd_norm += weight * CALL_VOLUME(local_values.table);
+      weighted_volume += weight * CALL_VOLUME(local_values.table);
       #if COMPUTE_F1_F2
       weight_norm += weight;
       #endif
+      if (effective_radius_type != 0) {
+        weighted_radius += weight * CALL_EFFECTIVE_RADIUS(effective_radius_type, local_values.table);
+      }
       BUILD_ROTATION();
 
@@ -745,6 +757,6 @@
         #ifdef USE_OPENCL
           #if COMPUTE_F1_F2
+            this_F2 += weight * F2;
             this_F1 += weight * F1;
-            this_F2 += weight * F2;
           #else
             this_result += weight * scattering;
@@ -752,6 +764,6 @@
         #else // !USE_OPENCL
           #if COMPUTE_F1_F2
-            result[q_index] += weight * F2;
-            result[q_index+nq+1] += weight * F1;
+            result[2*q_index+0] += weight * F2;
+            result[2*q_index+1] += weight * F1;
           #else
             result[q_index] += weight * scattering;
@@ -782,24 +794,32 @@
 #ifdef USE_OPENCL
   #if COMPUTE_F1_F2
-    result[q_index] = this_F2;
-    result[q_index+nq+1] = this_F1;
+    result[2*q_index+0] = this_F2;
+    result[2*q_index+1] = this_F1;
     if (q_index == 0) {
-      result[nq] = pd_norm;
-      result[2*nq+1] = weight_norm;
+      result[2*nq+0] = weight_norm;
+      result[2*nq+1] = weighted_volume;
+      result[2*nq+2] = weighted_radius;
     }
   #else
     result[q_index] = this_result;
-    if (q_index == 0) result[nq] = pd_norm;
+    if (q_index == 0) {
+      result[nq+0] = weight_norm;
+      result[nq+1] = weighted_volume;
+      result[nq+2] = weighted_radius;
+    }
   #endif
 
-//if (q_index == 0) printf("res: %g/%g\n", result[0], pd_norm);
+//if (q_index == 0) printf("res: %g/%g\n", result[0], weigthed_volume);
 #else // !USE_OPENCL
   #if COMPUTE_F1_F2
-    result[nq] = pd_norm;
-    result[2*nq+1] = weight_norm;
+    result[2*nq] = weight_norm;
+    result[2*nq+1] = weighted_volume;
+    result[2*nq+2] = weighted_radius;
   #else
-    result[nq] = pd_norm;
+    result[nq] = weight_norm;
+    result[nq+1] = weighted_volume;
+    result[nq+2] = weighted_radius;
   #endif
-//printf("res: %g/%g\n", result[0], pd_norm);
+//printf("res: %g/%g\n", result[0], weighted_volume);
 #endif // !USE_OPENCL
 
Index: sasmodels/kernelcl.py
===================================================================
--- sasmodels/kernelcl.py	(revision c036ddbd6dd028e64cda5692721d394c9845adaf)
+++ sasmodels/kernelcl.py	(revision 6e7ba14b7bd7549d79f388fe396c80945675b80b)
@@ -540,6 +540,6 @@
         # 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)
+        # plus 1 for the normalization value, plus another for R_eff
+        self.result = np.empty((q_input.nq+2)*num_returns, dtype)
 
         # Inputs and outputs for each kernel call
@@ -558,30 +558,5 @@
                      else np.float32)  # will never get here, so use np.float32
 
-    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):
+    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type):
         # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray
         context = self.queue.context
@@ -597,4 +572,5 @@
             details_b, values_b, self.q_input.q_b, self.result_b,
             self.real(cutoff),
+            np.uint32(effective_radius_type),
         ]
         #print("Calling OpenCL")
Index: sasmodels/kerneldll.py
===================================================================
--- sasmodels/kerneldll.py	(revision 2f8cbb9312dd1c060fcab5c78754e49ab71b15b7)
+++ sasmodels/kerneldll.py	(revision 6e7ba14b7bd7549d79f388fe396c80945675b80b)
@@ -315,5 +315,5 @@
 
         # int, int, int, int*, double*, double*, double*, double*, double
-        argtypes = [ct.c_int32]*3 + [ct.c_void_p]*4 + [float_type]
+        argtypes = [ct.c_int32]*3 + [ct.c_void_p]*4 + [float_type, ct.c_int32]
         names = [generate.kernel_name(self.info, variant)
                  for variant in ("Iq", "Iqxy", "Imagnetic")]
@@ -382,37 +382,13 @@
         self.dim = '2d' if q_input.is_2d else '1d'
         # 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)
+        nout = 2 if self.info.have_Fq else 1
+        # plus 3 weight, volume, radius
+        self.result = np.empty(q_input.nq*nout + 3, self.dtype)
         self.real = (np.float32 if self.q_input.dtype == generate.F32
                      else np.float64 if self.q_input.dtype == generate.F64
                      else np.float128)
 
-    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):
+    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type):
+        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray
         kernel = self.kernel[1 if magnetic else 0]
         args = [
@@ -425,4 +401,5 @@
             self.result.ctypes.data,   # results
             self.real(cutoff), # cutoff
+            effective_radius_type, # cutoff
         ]
         #print("Calling DLL")
Index: sasmodels/kernelpy.py
===================================================================
--- sasmodels/kernelpy.py	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
+++ sasmodels/kernelpy.py	(revision 6e7ba14b7bd7549d79f388fe396c80945675b80b)
@@ -12,4 +12,9 @@
 
 import numpy as np  # type: ignore
+from numpy import pi
+try:
+    from numpy import cbrt
+except ImportError:
+    def cbrt(x): return x ** (1.0/3.0)
 
 from .generate import F64
@@ -158,9 +163,16 @@
 
         # Generate a closure which calls the form_volume if it exists.
-        form_volume = model_info.form_volume
-        self._volume = ((lambda: form_volume(*volume_args)) if form_volume else
-                        (lambda: 1.0))
-
-    def __call__(self, call_details, values, cutoff, magnetic):
+        self._volume_args = volume_args
+        volume = model_info.form_volume
+        radius = model_info.effective_radius
+        self._volume = ((lambda: volume(*volume_args)) if volume
+                        else (lambda: 1.0))
+        self._radius = ((lambda mode: radius(mode, *volume_args)) if radius
+                        else (lambda mode: cbrt(0.75/pi*volume(*volume_args))) if volume
+                        else (lambda mode: 1.0))
+
+
+
+    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type):
         # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray
         if magnetic:
@@ -168,5 +180,7 @@
         #print("Calling python kernel")
         #call_details.show(values)
-        res = _loops(self._parameter_vector, self._form, self._volume,
+        radius = ((lambda: 0.0) if effective_radius_type == 0
+                  else (lambda: self._radius(effective_radius_type)))
+        res = _loops(self._parameter_vector, self._form, self._volume, radius,
                      self.q_input.nq, call_details, values, cutoff)
         return res
@@ -183,8 +197,9 @@
            form,          # type: Callable[[], np.ndarray]
            form_volume,   # type: Callable[[], float]
+           form_radius,   # type: Callable[[], float]
            nq,            # type: int
            call_details,  # type: details.CallDetails
            values,        # type: np.ndarray
-           cutoff         # type: float
+           cutoff,        # type: float
           ):
     # type: (...) -> None
@@ -209,5 +224,7 @@
     pd_weight = values[2+n_pars + call_details.num_weights:]
 
-    pd_norm = 0.0
+    weight_norm = 0.0
+    weighted_volume = 0.0
+    weighted_radius = 0.0
     partial_weight = np.NaN
     weight = np.NaN
@@ -246,9 +263,10 @@
             # update value and norm
             total += weight * Iq
-            pd_norm += weight * form_volume()
-
-    scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0)
-    background = values[1]
-    return scale*total + background
+            weight_norm += weight
+            weighted_volume += weight * form_volume()
+            weighted_radius += weight * form_radius()
+
+    result = np.hstack((total, weight_norm, weighted_volume, weighted_radius))
+    return result
 
 
Index: sasmodels/modelinfo.py
===================================================================
--- sasmodels/modelinfo.py	(revision c036ddbd6dd028e64cda5692721d394c9845adaf)
+++ sasmodels/modelinfo.py	(revision 6e7ba14b7bd7549d79f388fe396c80945675b80b)
@@ -789,12 +789,13 @@
     info.structure_factor = getattr(kernel_module, 'structure_factor', False)
     # TODO: find Fq by inspection
+    info.effective_radius_type = getattr(kernel_module, 'effective_radius_type', None)
     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', [])
     info.c_code = getattr(kernel_module, 'c_code', None)
+    info.ER = None  # CRUFT
+    info.VR = None  # CRUFT
     # TODO: check the structure of the tests
     info.tests = getattr(kernel_module, 'tests', [])
-    info.ER = getattr(kernel_module, 'ER', None) # type: ignore
-    info.VR = getattr(kernel_module, 'VR', None) # type: ignore
     info.form_volume = getattr(kernel_module, 'form_volume', None) # type: ignore
     info.Iq = getattr(kernel_module, 'Iq', None) # type: ignore
@@ -922,4 +923,7 @@
     #: use those functions.  Form factors are indicated by providing
     #: an :attr:`ER` function.
+    effective_radius_type = None   # type: List[str]
+    #: Returns the occupied volume and the total volume for each parameter set.
+    #: See :attr:`ER` for details on the parameters.
     source = None           # type: List[str]
     #: The set of tests that must pass.  The format of the tests is described
@@ -942,11 +946,4 @@
     #: each parameter set.  Multiplicity parameters will be received as
     #: arrays, with one row per polydispersity condition.
-    ER = None               # type: Optional[Callable[[np.ndarray], np.ndarray]]
-    #: Returns the occupied volume and the total volume for each parameter set.
-    #: See :attr:`ER` for details on the parameters.
-    VR = None               # type: Optional[Callable[[np.ndarray], Tuple[np.ndarray, np.ndarray]]]
-    #: Arbitrary C code containing supporting functions, etc., to be inserted
-    #: after everything in source.  This can include Iq and Iqxy functions with
-    #: the full function signature, including all parameters.
     c_code = None
     #: Returns the form volume for python-based models.  Form volume is needed
Index: sasmodels/models/ellipsoid.c
===================================================================
--- sasmodels/models/ellipsoid.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/ellipsoid.c	(revision 6e7ba14b7bd7549d79f388fe396c80945675b80b)
@@ -4,4 +4,46 @@
     return M_4PI_3*radius_polar*radius_equatorial*radius_equatorial;
 }
+
+
+static double
+radius_from_volume(double radius_polar, double radius_equatorial)
+{
+    return cbrt(radius_polar*radius_equatorial*radius_equatorial);
+}
+
+static double
+radius_from_curvature(double radius_polar, double radius_equatorial)
+{
+    // Trivial cases
+    if (radius_polar == radius_equatorial) return radius_polar;
+    if (radius_polar * radius_equatorial == 0.)  return 0.;
+
+    // see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449
+    const double ratio = (radius_polar < radius_equatorial 
+                          ? radius_polar / radius_equatorial
+                          : radius_equatorial / radius_polar);
+    const double e1 = sqrt(1.0 - ratio*ratio);
+    const double b1 = 1.0 + asin(e1) / (e1 * ratio);
+    const double bL = (1.0 + e1) / (1.0 - e1);
+    const double b2 = 1.0 + 0.5 * ratio * ratio / e1 * log(bL);
+    const double delta = 0.75 * b1 * b2;
+    const double ddd = 2.0 * (delta + 1.0) * radius_polar * radius_equatorial * radius_equatorial;
+    return 0.5 * cbrt(ddd);
+}
+
+static double
+effective_radius(int mode, double radius_polar, double radius_equatorial)
+{
+    if (mode == 1) {
+        return radius_from_curvature(radius_polar, radius_equatorial);
+    } else if (mode == 2) {
+        return radius_from_volume(radius_polar, radius_equatorial);
+    } else if (mode == 3) {
+        return (radius_polar < radius_equatorial ? radius_polar : radius_equatorial);
+    } else {
+        return (radius_polar > radius_equatorial ? radius_polar : radius_equatorial);
+    }
+}
+
 
 static void
Index: sasmodels/models/ellipsoid.py
===================================================================
--- sasmodels/models/ellipsoid.py	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/ellipsoid.py	(revision 6e7ba14b7bd7549d79f388fe396c80945675b80b)
@@ -164,25 +164,5 @@
 source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "ellipsoid.c"]
 have_Fq = True
-
-def ER(radius_polar, radius_equatorial):
-    # see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449
-    ee = np.empty_like(radius_polar)
-    idx = radius_polar > radius_equatorial
-    ee[idx] = (radius_polar[idx] ** 2 - radius_equatorial[idx] ** 2) / radius_polar[idx] ** 2
-    idx = radius_polar < radius_equatorial
-    ee[idx] = (radius_equatorial[idx] ** 2 - radius_polar[idx] ** 2) / radius_equatorial[idx] ** 2
-    idx = radius_polar == radius_equatorial
-    ee[idx] = 2 * radius_polar[idx]
-    valid = (radius_polar * radius_equatorial != 0)
-    bd = 1.0 - ee[valid]
-    e1 = np.sqrt(ee[valid])
-    b1 = 1.0 + np.arcsin(e1) / (e1 * np.sqrt(bd))
-    bL = (1.0 + e1) / (1.0 - e1)
-    b2 = 1.0 + bd / 2 / e1 * np.log(bL)
-    delta = 0.75 * b1 * b2
-
-    ddd = np.zeros_like(radius_polar)
-    ddd[valid] = 2.0 * (delta + 1.0) * radius_polar * radius_equatorial ** 2
-    return 0.5 * ddd ** (1.0 / 3.0)
+effective_radius_type = ["average curvature", "equivalent sphere", "min radius", "max radius"]
 
 def random():
Index: sasmodels/product.py
===================================================================
--- sasmodels/product.py	(revision 6d90684c44557e58fda7371113cb02cefa6e334a)
+++ sasmodels/product.py	(revision 6e7ba14b7bd7549d79f388fe396c80945675b80b)
@@ -35,8 +35,18 @@
 #]
 
-ER_ID = "radius_effective"
-VF_ID = "volfraction"
-BETA_DEFINITION = ("beta_mode", "", 0, [["P*S"],["P*(1+beta*(S-1))"]], "",
-                   "Structure factor dispersion calculation mode")
+RADIUS_ID = "radius_effective"
+VOLFRAC_ID = "volfraction"
+def make_extra_pars(p_info):
+    pars = []
+    if p_info.have_Fq:
+        par = Parameter("structure_factor_mode", "", 0, [["P*S","P*(1+beta*(S-1))"]], "",
+                        "Structure factor calculation")
+        pars.append(par)
+    if p_info.effective_radius_type is not None:
+        par = Parameter("radius_effective_mode", "", 0,
+                        [["unconstrained"] + p_info.effective_radius_type],
+                        "", "Effective radius calculation")
+        pars.append(par)
+    return pars
 
 # TODO: core_shell_sphere model has suppressed the volume ratio calculation
@@ -52,9 +62,9 @@
     # have any magnetic parameters
     if not len(s_info.parameters.kernel_parameters) >= 2:
-        raise TypeError("S needs {} and {} as its first parameters".format(ER_ID, VF_ID))
-    if not s_info.parameters.kernel_parameters[0].id == ER_ID:
-        raise TypeError("S needs {} as first parameter".format(ER_ID))
-    if not s_info.parameters.kernel_parameters[1].id == VF_ID:
-        raise TypeError("S needs {} as second parameter".format(VF_ID))
+        raise TypeError("S needs {} and {} as its first parameters".format(RADIUS_ID, VOLFRAC_ID))
+    if not s_info.parameters.kernel_parameters[0].id == RADIUS_ID:
+        raise TypeError("S needs {} as first parameter".format(RADIUS_ID))
+    if not s_info.parameters.kernel_parameters[1].id == VOLFRAC_ID:
+        raise TypeError("S needs {} as second parameter".format(VOLFRAC_ID))
     if not s_info.parameters.magnetism_index == []:
         raise TypeError("S should not have SLD parameters")
@@ -62,11 +72,10 @@
     s_id, s_name, s_pars = s_info.id, s_info.name, s_info.parameters
 
-    # Create list of parameters for the combined model.  Skip the first
-    # parameter of S, which we verified above is effective radius.  If there
+    # Create list of parameters for the combined model.  If there
     # are any names in P that overlap with those in S, modify the name in S
     # to distinguish it.
     p_set = set(p.id for p in p_pars.kernel_parameters)
     s_list = [(_tag_parameter(par) if par.id in p_set else par)
-              for par in s_pars.kernel_parameters[1:]]
+              for par in s_pars.kernel_parameters]
     # Check if still a collision after renaming.  This could happen if for
     # example S has volfrac and P has both volfrac and volfrac_S.
@@ -74,13 +83,16 @@
         raise TypeError("name collision: P has P.name and P.name_S while S has S.name")
 
+    # make sure effective radius is not a polydisperse parameter in product
+    s_list[0] = copy(s_list[0])
+    s_list[0].polydisperse = False
+
     translate_name = dict((old.id, new.id) for old, new
-                          in zip(s_pars.kernel_parameters[1:], s_list))
-    beta = [Parameter(*BETA_DEFINITION)] if p_info.have_Fq else []
-    combined_pars = p_pars.kernel_parameters + s_list + beta
+                          in zip(s_pars.kernel_parameters, s_list))
+    combined_pars = p_pars.kernel_parameters + s_list + make_extra_pars(p_info)
     parameters = ParameterTable(combined_pars)
     parameters.max_pd = p_pars.max_pd + s_pars.max_pd
     def random():
         combined_pars = p_info.random()
-        s_names = set(par.id for par in s_pars.kernel_parameters[1:])
+        s_names = set(par.id for par in s_pars.kernel_parameters)
         combined_pars.update((translate_name[k], v)
                              for k, v in s_info.random().items()
@@ -225,17 +237,6 @@
         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
-        # the number of 'S' parameters is one less than expected.  The
-        # computed effective radius needs to be added into the weights
-        # vector, especially since it is a polydisperse parameter in the
-        # stand-alone structure factor models.  We will added it at the
-        # end so the remaining offsets don't need to change.
-        s_npars = s_info.parameters.npars-1
+        s_npars = s_info.parameters.npars
         s_length = call_details.length[p_npars:p_npars+s_npars]
         s_offset = call_details.offset[p_npars:p_npars+s_npars]
@@ -243,16 +244,9 @@
         s_offset = np.hstack((nweights, s_offset))
         s_details = make_details(s_info, s_length, s_offset, nweights+1)
-        v, w = weights[:nweights], weights[nweights:]
         s_values = [
-            # scale=1, background=0, radius_effective=p_er, volfraction=s_vr
-            [1., 0., p_er, s_vr],
-            # structure factor parameters start after scale, background and
-            # all the form factor parameters.  Skip the volfraction parameter
-            # as well, since it is computed elsewhere, and go to the end of the
-            # parameter list.
-            values[2+p_npars+1:2+p_npars+s_npars],
-            # no magnetism parameters to include for S
-            # add er into the (value, weights) pairs
-            v, [p_er], w, [1.0]
+            # scale=1, background=0,
+            [1., 0.],
+            values[2+p_npars:2+p_npars+s_npars],
+            weights,
         ]
         spacer = (32 - sum(len(v) for v in s_values)%32)%32
@@ -261,12 +255,24 @@
 
         # beta mode is the first parameter after the structure factor pars
-        beta_index = 2+p_npars+s_npars
-        beta_mode = values[beta_index]
+        extra_offset = 2+p_npars+s_npars
+        if p_info.have_Fq:
+            beta_mode = values[extra_offset]
+            extra_offset += 1
+        else:
+            beta_mode = 0
+        if p_info.effective_radius_type is not None:
+            effective_radius_type = values[extra_offset]
+            extra_offset += 1
+        else:
+            effective_radius_type = 0
 
         # Call the kernels
-        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)
+            F1, F2, volume_avg, effective_radius = self.p_kernel.beta(
+                p_details, p_values, cutoff, magnetic, effective_radius_type)
+            if effective_radius_type > 0:
+                s_values[2] = s_values[2+s_npars+s_offset[0]] = effective_radius
+            s_result = self.s_kernel.Iq(s_details, s_values, cutoff, False)
             combined_scale = scale*volfrac/volume_avg
             # Define lazy results based on intermediate values.
@@ -297,5 +303,9 @@
             final_result = combined_scale*(F2 + (F1**2)*(s_result - 1)) + background
         else:
-            p_result = self.p_kernel.Iq(p_details, p_values, cutoff, magnetic)
+            p_result, effective_radius = self.p_kernel.Pq_Reff(
+                p_details, p_values, cutoff, magnetic, effective_radius_type)
+            if effective_radius_type > 0:
+                s_values[2] = s_values[2+s_npars+s_offset[0]] = effective_radius
+            s_result = self.s_kernel.Iq(s_details, s_values, cutoff, False)
             # remember the parts for plotting later
             self.results = [p_result, s_result]
Index: sasmodels/sasview_model.py
===================================================================
--- sasmodels/sasview_model.py	(revision d533590d9857d7aef962616cfdaa6ecfc2e6751d)
+++ sasmodels/sasview_model.py	(revision 6e7ba14b7bd7549d79f388fe396c80945675b80b)
@@ -272,5 +272,5 @@
     attrs['category'] = model_info.category
     attrs['is_structure_factor'] = model_info.structure_factor
-    attrs['is_form_factor'] = model_info.ER is not None
+    attrs['is_form_factor'] = model_info.effective_radius_type is not None
     attrs['is_multiplicity_model'] = variants[0] > 1
     attrs['multiplicity_info'] = variants
