Index: explore/beta/sasfit_compare.py
===================================================================
--- explore/beta/sasfit_compare.py	(revision 7b0abf8a1cf121fb3e4733dd887a921fb598a7a0)
+++ explore/beta/sasfit_compare.py	(revision 2a12351b8b7f57f6668c3974b81f6be423c68ba3)
@@ -208,4 +208,5 @@
     F1, F2 = np.zeros_like(q), np.zeros_like(q)
     for k, Rpk in enumerate(Rp_val):
+        print("ellipsoid cycle", k, "of", len(Rp_val))
         for i, Rei in enumerate(Re_val):
             theory = ellipsoid_theta(q, Rpk, Rei, sld, sld_solvent)
@@ -343,6 +344,6 @@
     #Sq = Iq/Pq
     #Iq = None#= Sq = None
-    r = I._kernel.results
-    return Theory(Q=q, F1=None, F2=None, P=Pq, S=None, I=None, Seff=r[1], Ibeta=Iq)
+    r = dict(I._kernel.results())
+    return Theory(Q=q, F1=None, F2=None, P=Pq, S=None, I=None, Seff=r["S_eff(Q)"], Ibeta=Iq)
 
 def compare(title, target, actual, fields='F1 F2 P S I Seff Ibeta'):
@@ -401,11 +402,11 @@
         radius_polar=20, radius_equatorial=400, sld=4, sld_solvent=1,
         background=0,
-        radius_polar_pd=.1, radius_polar_pd_type=pd_type,
-        radius_equatorial_pd=.1, radius_equatorial_pd_type=pd_type,
+        radius_polar_pd=0.1, radius_polar_pd_type=pd_type,
+        radius_equatorial_pd=0.1, radius_equatorial_pd_type=pd_type,
         volfraction=0.15,
         radius_effective=270.7543927018,
         #radius_effective=12.59921049894873,
         )
-    target = sasmodels_theory(q, model, beta_mode=1, **pars)
+    target = sasmodels_theory(q, model, effective_radius_mode=0, structure_factor_mode=1, **pars)
     actual = ellipsoid_pe(q, norm='sasview', **pars)
     title = " ".join(("sasmodels", model, pd_type))
Index: sasmodels/core.py
===================================================================
--- sasmodels/core.py	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/core.py	(revision 53998094412e4e529314e9fc45914cc59e148317)
@@ -363,6 +363,7 @@
     model = load_model("cylinder@hardsphere*sphere")
     actual = [p.name for p in model.info.parameters.kernel_parameters]
-    target = ("sld sld_solvent radius length theta phi volfraction"
-              " beta_mode A_sld A_sld_solvent A_radius").split()
+    target = ("sld sld_solvent radius length theta phi"
+              " radius_effective volfraction structure_factor_mode"
+              " A_sld A_sld_solvent A_radius").split()
     assert target == actual, "%s != %s"%(target, actual)
 
Index: sasmodels/generate.py
===================================================================
--- sasmodels/generate.py	(revision 6e45516f19370dc8d3a8ff3fa00959f6e27d9e19)
+++ sasmodels/generate.py	(revision 53998094412e4e529314e9fc45914cc59e148317)
@@ -802,7 +802,10 @@
                               for p in partable.kernel_parameters))
     # Define the function calls
+    call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(mode, v) 0.0"
     if partable.form_volume_parameters:
         refs = _call_pars("_v.", partable.form_volume_parameters)
         call_volume = "#define CALL_VOLUME(_v) form_volume(%s)"%(",".join(refs))
+        if model_info.effective_radius_type:
+            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
@@ -811,4 +814,5 @@
         call_volume = "#define CALL_VOLUME(v) 1.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 53998094412e4e529314e9fc45914cc59e148317)
@@ -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)
+        return Pq
+    __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 and self.dim == '1d' 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
+        F2 = self.result[0:nout*self.q_input.nq:nout]
+        scale = values[0]/(weighted_volume if weighted_volume != 0.0 else 1.0)
+        background = values[1]
+        Pq = scale*F2 + background
+        #print("scale",scale,background)
+        return Pq, 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_header.c
===================================================================
--- sasmodels/kernel_header.c	(revision 108e70e22154db937a6d46f66be6dd70408cce6d)
+++ sasmodels/kernel_header.c	(revision 296c52b88b3fea0e92c54ab768d5f1f91d8598a5)
@@ -68,4 +68,5 @@
          #define NEED_EXPM1
          #define NEED_TGAMMA
+         #define NEED_CBRT
          // expf missing from windows?
          #define expf exp
@@ -85,4 +86,8 @@
 #  define pown(a,b) pow(a,b)
 #endif // !USE_OPENCL
+
+#if defined(NEED_CBRT)
+   #define cbrt(_x) pow(_x, 0.33333333333333333333333)
+#endif
 
 #if defined(NEED_EXPM1)
Index: sasmodels/kernel_iq.c
===================================================================
--- sasmodels/kernel_iq.c	(revision 705307783b90bea1a66cdf1ff38974433eb673cb)
+++ sasmodels/kernel_iq.c	(revision 705307783b90bea1a66cdf1ff38974433eb673cb)
@@ -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
+    int32_t 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 53998094412e4e529314e9fc45914cc59e148317)
@@ -481,13 +481,12 @@
         # at this point, so instead using 32, which is good on the set of
         # architectures tested so far.
+        extra_q = 3  # total weight, weighted volume and weighted radius
         if self.is_2d:
-            # Note: 16 rather than 15 because result is 1 longer than input.
-            width = ((self.nq+16)//16)*16
+            width = ((self.nq+15+extra_q)//16)*16
             self.q = np.empty((width, 2), dtype=dtype)
             self.q[:self.nq, 0] = q_vectors[0]
             self.q[:self.nq, 1] = q_vectors[1]
         else:
-            # Note: 32 rather than 31 because result is 1 longer than input.
-            width = ((self.nq+32)//32)*32
+            width = ((self.nq+31+extra_q)//32)*32
             self.q = np.empty(width, dtype=dtype)
             self.q[:self.nq] = q_vectors[0]
@@ -539,7 +538,7 @@
         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, dtype)
+        nout = 2 if self.info.have_Fq and self.dim == '1d' else 1
+        # plus 3 weight, volume, radius
+        self.result = np.empty(q_input.nq*nout + 3, self.dtype)
 
         # Inputs and outputs for each kernel call
@@ -549,5 +548,5 @@
 
         self.result_b = cl.Buffer(self.queue.context, mf.READ_WRITE,
-                                  q_input.global_size[0] * num_returns * dtype.itemsize)
+                                  q_input.global_size[0] * nout * dtype.itemsize)
         self.q_input = q_input # allocated by GpuInput above
 
@@ -558,30 +557,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 +571,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 53998094412e4e529314e9fc45914cc59e148317)
@@ -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,7 +180,9 @@
         #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)
-        return res
+        radius = ((lambda: 0.0) if effective_radius_type == 0
+                  else (lambda: self._radius(effective_radius_type)))
+        self.result = _loops(
+            self._parameter_vector, self._form, self._volume, radius,
+            self.q_input.nq, call_details, values, cutoff)
 
     def release(self):
@@ -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
@@ -198,57 +213,69 @@
     #                                                              #
     ################################################################
+
+    # WARNING: Trickery ahead
+    # The parameters[] vector is embedded in the closures for form(),
+    # form_volume() and form_radius().  We set the initial vector from
+    # the values for the model parameters. As we loop through the polydispesity
+    # mesh, we update the components with the polydispersity values before
+    # calling the respective functions.
     n_pars = len(parameters)
     parameters[:] = values[2:n_pars+2]
+
     if call_details.num_active == 0:
-        pd_norm = float(form_volume())
-        scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0)
-        background = values[1]
-        return scale*form() + background
-
-    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
-    partial_weight = np.NaN
-    weight = np.NaN
-
-    p0_par = call_details.pd_par[0]
-    p0_length = call_details.pd_length[0]
-    p0_index = p0_length
-    p0_offset = call_details.pd_offset[0]
-
-    pd_par = call_details.pd_par[:call_details.num_active]
-    pd_offset = call_details.pd_offset[:call_details.num_active]
-    pd_stride = call_details.pd_stride[:call_details.num_active]
-    pd_length = call_details.pd_length[:call_details.num_active]
-
-    total = np.zeros(nq, 'd')
-    for loop_index in range(call_details.num_eval):
-        # update polydispersity parameter values
-        if p0_index == p0_length:
-            pd_index = (loop_index//pd_stride)%pd_length
-            parameters[pd_par] = pd_value[pd_offset+pd_index]
-            partial_weight = np.prod(pd_weight[pd_offset+pd_index][1:])
-            p0_index = loop_index%p0_length
-
-        weight = partial_weight * pd_weight[p0_offset + p0_index]
-        parameters[p0_par] = pd_value[p0_offset + p0_index]
-        p0_index += 1
-        if weight > cutoff:
-            # Call the scattering function
-            # Assume that NaNs are only generated if the parameters are bad;
-            # exclude all q for that NaN.  Even better would be to have an
-            # INVALID expression like the C models, but that is too expensive.
-            Iq = np.asarray(form(), 'd')
-            if np.isnan(Iq).any():
-                continue
-
-            # 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
+        total = form()
+        weight_norm = 1.0
+        weighted_volume = form_volume()
+        weighted_radius = form_radius()
+
+    else:
+        pd_value = values[2+n_pars:2+n_pars + call_details.num_weights]
+        pd_weight = values[2+n_pars + call_details.num_weights:]
+
+        weight_norm = 0.0
+        weighted_volume = 0.0
+        weighted_radius = 0.0
+        partial_weight = np.NaN
+        weight = np.NaN
+
+        p0_par = call_details.pd_par[0]
+        p0_length = call_details.pd_length[0]
+        p0_index = p0_length
+        p0_offset = call_details.pd_offset[0]
+
+        pd_par = call_details.pd_par[:call_details.num_active]
+        pd_offset = call_details.pd_offset[:call_details.num_active]
+        pd_stride = call_details.pd_stride[:call_details.num_active]
+        pd_length = call_details.pd_length[:call_details.num_active]
+
+        total = np.zeros(nq, 'd')
+        for loop_index in range(call_details.num_eval):
+            # update polydispersity parameter values
+            if p0_index == p0_length:
+                pd_index = (loop_index//pd_stride)%pd_length
+                parameters[pd_par] = pd_value[pd_offset+pd_index]
+                partial_weight = np.prod(pd_weight[pd_offset+pd_index][1:])
+                p0_index = loop_index%p0_length
+
+            weight = partial_weight * pd_weight[p0_offset + p0_index]
+            parameters[p0_par] = pd_value[p0_offset + p0_index]
+            p0_index += 1
+            if weight > cutoff:
+                # Call the scattering function
+                # Assume that NaNs are only generated if the parameters are bad;
+                # exclude all q for that NaN.  Even better would be to have an
+                # INVALID expression like the C models, but that is expensive.
+                Iq = np.asarray(form(), 'd')
+                if np.isnan(Iq).any():
+                    continue
+
+                # update value and norm
+                total += weight * Iq
+                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 7b9e4ddca7ad3551d8fe3696d2afe75b0018599e)
+++ sasmodels/modelinfo.py	(revision 53998094412e4e529314e9fc45914cc59e148317)
@@ -790,12 +790,14 @@
     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.effective_radius = getattr(kernel_module, 'effective_radius', 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
@@ -923,4 +925,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
@@ -943,11 +948,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/barbell.c
===================================================================
--- sasmodels/models/barbell.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/barbell.c	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -62,4 +62,32 @@
 }
 
+static double
+radius_from_volume(double radius_bell, double radius, double length)
+{
+    const double vol_barbell = form_volume(radius_bell,radius,length);
+    return cbrt(0.75*vol_barbell/M_PI);
+}
+
+static double
+radius_from_totallength(double radius_bell, double radius, double length)
+{
+    const double hdist = sqrt(square(radius_bell) - square(radius));
+    return 0.5*length + hdist + radius_bell;
+}
+
+static double
+effective_radius(int mode, double radius_bell, double radius, double length)
+{
+    if (mode == 1) {
+        return radius_from_volume(radius_bell, radius , length);
+    } else if (mode == 2) {
+        return radius;
+    } else if (mode == 3) {
+        return 0.5*length;
+    } else {
+        return radius_from_totallength(radius_bell,radius,length);
+    }
+}
+
 static void
 Fq(double q,double *F1, double *F2, double sld, double solvent_sld,
Index: sasmodels/models/barbell.py
===================================================================
--- sasmodels/models/barbell.py	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/barbell.py	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -117,4 +117,5 @@
 source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "barbell.c"]
 have_Fq = True
+effective_radius_type = ["equivalent sphere","radius","half length","half total length"]
 
 def random():
Index: sasmodels/models/capped_cylinder.c
===================================================================
--- sasmodels/models/capped_cylinder.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/capped_cylinder.c	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -84,4 +84,32 @@
 }
 
+static double
+radius_from_volume(double radius, double radius_cap, double length)
+{
+    const double vol_cappedcyl = form_volume(radius,radius_cap,length);
+    return cbrt(0.75*vol_cappedcyl/M_PI);
+}
+
+static double
+radius_from_totallength(double radius, double radius_cap, double length)
+{
+    const double hc = radius_cap - sqrt(radius_cap*radius_cap - radius*radius);
+    return 0.5*length + hc;
+}
+
+static double
+effective_radius(int mode, double radius, double radius_cap, double length)
+{
+    if (mode == 1) {
+        return radius_from_volume(radius, radius_cap, length);
+    } else if (mode == 2) {
+        return radius;
+    } else if (mode == 3) {
+        return 0.5*length;
+    } else {
+        return radius_from_totallength(radius, radius_cap,length);
+    }
+}
+
 static void
 Fq(double q,double *F1, double *F2, double sld, double solvent_sld,
Index: sasmodels/models/capped_cylinder.py
===================================================================
--- sasmodels/models/capped_cylinder.py	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/capped_cylinder.py	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -137,4 +137,5 @@
 source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "capped_cylinder.c"]
 have_Fq = True
+effective_radius_type = ["equivalent sphere","radius","half length","half total length"]
 
 def random():
Index: sasmodels/models/core_multi_shell.c
===================================================================
--- sasmodels/models/core_multi_shell.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/core_multi_shell.c	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -19,4 +19,31 @@
 }
 
+static double
+outer_radius(double core_radius, double fp_n, double thickness[])
+{
+  double r = core_radius;
+  int n = (int)(fp_n+0.5);
+  for (int i=0; i < n; i++) {
+    r += thickness[i];
+  }
+  return r;
+}
+
+static double
+effective_radius(int mode, double core_radius, double fp_n, double thickness[])
+{
+    if (mode == 1) {
+        double r = core_radius;
+        int n = (int)(fp_n+0.5);
+        for (int i=0; i < n; i++) {
+            r += thickness[i];
+        }
+        return r;
+        //return outer_radius(core_radius,fp_n,thickness);
+    } else {
+        return core_radius;
+    }
+}
+
 static void
 Fq(double q, double *F1, double *F2, double core_sld, double core_radius,
Index: sasmodels/models/core_multi_shell.py
===================================================================
--- sasmodels/models/core_multi_shell.py	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/core_multi_shell.py	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -100,4 +100,5 @@
 source = ["lib/sas_3j1x_x.c", "core_multi_shell.c"]
 have_Fq = True
+effective_radius_type = ["outer radius", "core radius"]
 
 def random():
@@ -144,8 +145,8 @@
     return np.asarray(z), np.asarray(rho)
 
-def ER(radius, n, thickness):
-    """Effective radius"""
-    n = int(n[0]+0.5)  # n is a control parameter and is not polydisperse
-    return np.sum(thickness[:n], axis=0) + radius
+#def ER(radius, n, thickness):
+#    """Effective radius"""
+#    n = int(n[0]+0.5)  # n is a control parameter and is not polydisperse
+#    return np.sum(thickness[:n], axis=0) + radius
 
 demo = dict(sld_core=6.4,
Index: sasmodels/models/core_shell_bicelle.c
===================================================================
--- sasmodels/models/core_shell_bicelle.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/core_shell_bicelle.c	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -34,4 +34,33 @@
 
     return t;
+}
+
+static double
+radius_from_volume(double radius, double thick_rim, double thick_face, double length)
+{
+    const double volume_bicelle = form_volume(radius,thick_rim,thick_face,length);
+    return cbrt(0.75*volume_bicelle/M_PI);
+}
+
+static double
+radius_from_diagonal(double radius, double thick_rim, double thick_face, double length)
+{
+    const double radius_tot = radius + thick_rim;
+    const double length_tot = length + 2.0*thick_face;
+    return sqrt(radius_tot*radius_tot + 0.25*length_tot*length_tot);
+}
+
+static double
+effective_radius(int mode, double radius, double thick_rim, double thick_face, double length)
+{
+    if (mode == 1) {
+        return radius_from_volume(radius, thick_rim, thick_face, length);
+    } else if (mode == 2) {
+        return radius + thick_rim;
+    } else if (mode == 3) {
+        return 0.5*length + thick_face;
+    } else {
+        return radius_from_diagonal(radius,thick_rim,thick_face,length);
+    }
 }
 
Index: sasmodels/models/core_shell_bicelle.py
===================================================================
--- sasmodels/models/core_shell_bicelle.py	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/core_shell_bicelle.py	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -155,4 +155,5 @@
           "core_shell_bicelle.c"]
 have_Fq = True
+effective_radius_type = ["equivalent sphere","outer rim radius", "half outer thickness","half diagonal"]
 
 def random():
Index: sasmodels/models/core_shell_bicelle_elliptical.c
===================================================================
--- sasmodels/models/core_shell_bicelle_elliptical.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/core_shell_bicelle_elliptical.c	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -8,4 +8,38 @@
 {
     return M_PI*(r_minor+thick_rim)*(r_minor*x_core+thick_rim)*(length+2.0*thick_face);
+}
+
+static double
+radius_from_volume(double r_minor, double x_core, double thick_rim, double thick_face, double length)
+{
+    const double volume_bicelle = form_volume(r_minor, x_core, thick_rim,thick_face,length);
+    return cbrt(0.75*volume_bicelle/M_PI);
+}
+
+static double
+radius_from_diagonal(double r_minor, double x_core, double thick_rim, double thick_face, double length)
+{
+    const double radius_max = (x_core < 1.0 ? r_minor : x_core*r_minor);
+    const double radius_max_tot = radius_max + thick_rim;
+    const double length_tot = length + 2.0*thick_face;
+    return sqrt(radius_max_tot*radius_max_tot + 0.25*length_tot*length_tot);
+}
+
+static double
+effective_radius(int mode, double r_minor, double x_core, double thick_rim, double thick_face, double length)
+{
+    if (mode == 1) {
+        return radius_from_volume(r_minor, x_core, thick_rim, thick_face, length);
+    } else if (mode == 2) {
+        return 0.5*r_minor*(1.0 + x_core) + thick_rim;
+    } else if (mode == 3) {
+        return (x_core < 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim);
+    } else if (mode == 4) {
+        return (x_core > 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim);
+    } else if (mode ==5) {
+        return 0.5*length + thick_face;
+    } else {
+        return radius_from_diagonal(r_minor,x_core,thick_rim,thick_face,length);
+    }
 }
 
Index: sasmodels/models/core_shell_bicelle_elliptical.py
===================================================================
--- sasmodels/models/core_shell_bicelle_elliptical.py	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/core_shell_bicelle_elliptical.py	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -147,4 +147,6 @@
           "core_shell_bicelle_elliptical.c"]
 have_Fq = True
+effective_radius_type = ["equivalent sphere","outer rim average radius","outer rim min radius",
+                         "outer max radius", "half outer thickness","half diagonal"]
 
 def random():
Index: sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c
===================================================================
--- sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -9,4 +9,38 @@
     return M_PI*(  (r_minor + thick_rim)*(r_minor*x_core + thick_rim)* length +
                  square(r_minor)*x_core*2.0*thick_face  );
+}
+
+static double
+radius_from_volume(double r_minor, double x_core, double thick_rim, double thick_face, double length)
+{
+    const double volume_bicelle = form_volume(r_minor, x_core, thick_rim,thick_face,length);
+    return cbrt(0.75*volume_bicelle/M_PI);
+}
+
+static double
+radius_from_diagonal(double r_minor, double x_core, double thick_rim, double thick_face, double length)
+{
+    const double radius_max = (x_core < 1.0 ? r_minor : x_core*r_minor);
+    const double radius_max_tot = radius_max + thick_rim;
+    const double length_tot = length + 2.0*thick_face;
+    return sqrt(radius_max_tot*radius_max_tot + 0.25*length_tot*length_tot);
+}
+
+static double
+effective_radius(int mode, double r_minor, double x_core, double thick_rim, double thick_face, double length)
+{
+    if (mode == 1) {
+        return radius_from_volume(r_minor, x_core, thick_rim, thick_face, length);
+    } else if (mode == 2) {
+        return 0.5*r_minor*(1.0 + x_core) + thick_rim;
+    } else if (mode == 3) {
+        return (x_core < 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim);
+    } else if (mode == 4) {
+        return (x_core > 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim);
+    } else if (mode ==5) {
+        return 0.5*length + thick_face;
+    } else {
+        return radius_from_diagonal(r_minor,x_core,thick_rim,thick_face,length);
+    }
 }
 
Index: sasmodels/models/core_shell_bicelle_elliptical_belt_rough.py
===================================================================
--- sasmodels/models/core_shell_bicelle_elliptical_belt_rough.py	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/core_shell_bicelle_elliptical_belt_rough.py	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -160,4 +160,6 @@
           "core_shell_bicelle_elliptical_belt_rough.c"]
 have_Fq = True
+effective_radius_type = ["equivalent sphere","outer rim average radius","outer rim min radius",
+                         "outer max radius", "half outer thickness","half diagonal"]
 
 demo = dict(scale=1, background=0,
Index: sasmodels/models/core_shell_cylinder.c
===================================================================
--- sasmodels/models/core_shell_cylinder.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/core_shell_cylinder.c	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -11,4 +11,37 @@
 {
     return M_PI*square(radius+thickness)*(length+2.0*thickness);
+}
+
+static double
+radius_from_volume(double radius, double thickness, double length)
+{
+    const double volume_outer_cyl = form_volume(radius,thickness,length);
+    return cbrt(0.75*volume_outer_cyl/M_PI);
+}
+
+static double
+radius_from_diagonal(double radius, double thickness, double length)
+{
+    const double radius_outer = radius + thickness;
+    const double length_outer = length + thickness;
+    return sqrt(radius_outer*radius_outer + 0.25*length_outer*length_outer);
+}
+
+static double
+effective_radius(int mode, double radius, double thickness, double length)
+{
+    if (mode == 1) {
+        return radius_from_volume(radius, thickness, length);
+    } else if (mode == 2) {
+        return radius + thickness;
+    } else if (mode == 3) {
+        return 0.5*length + thickness;
+    } else if (mode == 4) {
+        return (radius < 0.5*length ? radius + thickness : 0.5*length + thickness);
+    } else if (mode == 5) {
+        return (radius > 0.5*length ? radius + thickness : 0.5*length + thickness);
+    } else {
+        return radius_from_diagonal(radius,thickness,length);
+    }
 }
 
Index: sasmodels/models/core_shell_cylinder.py
===================================================================
--- sasmodels/models/core_shell_cylinder.py	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/core_shell_cylinder.py	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -126,13 +126,15 @@
 source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "core_shell_cylinder.c"]
 have_Fq = True
+effective_radius_type = ["equivalent sphere","outer radius","half outer length","half min outer dimension",
+                         "half max outer dimension","half outer diagonal"]
 
-def ER(radius, thickness, length):
-    """
-    Returns the effective radius used in the S*P calculation
-    """
-    radius = radius + thickness
-    length = length + 2 * thickness
-    ddd = 0.75 * radius * (2 * radius * length + (length + radius) * (length + pi * radius))
-    return 0.5 * (ddd) ** (1. / 3.)
+#def ER(radius, thickness, length):
+#    """
+#    Returns the effective radius used in the S*P calculation
+#    """
+#    radius = radius + thickness
+#    length = length + 2 * thickness
+#    ddd = 0.75 * radius * (2 * radius * length + (length + radius) * (length + pi * radius))
+#    return 0.5 * (ddd) ** (1. / 3.)
 
 def VR(radius, thickness, length):
Index: sasmodels/models/core_shell_ellipsoid.c
===================================================================
--- sasmodels/models/core_shell_ellipsoid.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/core_shell_ellipsoid.c	(revision 3c60146f44ec3507ebbe0e1e7c740546e003f5b2)
@@ -36,4 +36,50 @@
     double vol = M_4PI_3*equat_shell*equat_shell*polar_shell;
     return vol;
+}
+
+static double
+radius_from_volume(double radius_equat_core, double x_core, double thick_shell, double x_polar_shell)
+{
+    const double volume_ellipsoid = form_volume(radius_equat_core, x_core, thick_shell, x_polar_shell);
+    return cbrt(0.75*volume_ellipsoid/M_PI);
+}
+
+static double
+radius_from_curvature(double radius_equat_core, double x_core, double thick_shell, double x_polar_shell)
+{
+    // Trivial cases
+    if (1.0 == x_core && 1.0 == x_polar_shell) return radius_equat_core + thick_shell;
+    if ((radius_equat_core + thick_shell)*(radius_equat_core*x_core + thick_shell*x_polar_shell) == 0.)  return 0.;
+
+    // see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449
+    const double radius_equat_tot = radius_equat_core + thick_shell;
+    const double radius_polar_tot = radius_equat_core*x_core + thick_shell*x_polar_shell;
+    const double ratio = (radius_polar_tot < radius_equat_tot
+                          ? radius_polar_tot / radius_equat_tot
+                          : radius_equat_tot / radius_polar_tot);
+    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_tot * radius_equat_tot * radius_equat_tot;
+    return 0.5 * cbrt(ddd);
+}
+
+static double
+effective_radius(int mode, double radius_equat_core, double x_core, double thick_shell, double x_polar_shell)
+{
+    const double radius_equat_tot = radius_equat_core + thick_shell;
+    const double radius_polar_tot = radius_equat_core*x_core + thick_shell*x_polar_shell;
+
+    if (mode == 1) {
+        return radius_from_volume(radius_equat_core, x_core, thick_shell, x_polar_shell);
+    } else if (mode == 2) {
+        return radius_from_curvature(radius_equat_core, x_core, thick_shell, x_polar_shell);
+    } else if (mode == 3) {
+        return (radius_polar_tot < radius_equat_tot ? radius_polar_tot : radius_equat_tot);
+    } else {
+        return (radius_polar_tot > radius_equat_tot ? radius_polar_tot : radius_equat_tot);
+    }
 }
 
Index: sasmodels/models/core_shell_ellipsoid.py
===================================================================
--- sasmodels/models/core_shell_ellipsoid.py	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/core_shell_ellipsoid.py	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -146,13 +146,14 @@
 source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "core_shell_ellipsoid.c"]
 have_Fq = True
-
-def ER(radius_equat_core, x_core, thick_shell, x_polar_shell):
-    """
-        Returns the effective radius used in the S*P calculation
-    """
-    from .ellipsoid import ER as ellipsoid_ER
-    polar_outer = radius_equat_core*x_core + thick_shell*x_polar_shell
-    equat_outer = radius_equat_core + thick_shell
-    return ellipsoid_ER(polar_outer, equat_outer)
+effective_radius_type = ["equivalent sphere","average outer curvature", "min outer radius", "max outer radius"]
+
+#def ER(radius_equat_core, x_core, thick_shell, x_polar_shell):
+#    """
+#        Returns the effective radius used in the S*P calculation
+#    """
+#    from .ellipsoid import ER as ellipsoid_ER
+#    polar_outer = radius_equat_core*x_core + thick_shell*x_polar_shell
+#    equat_outer = radius_equat_core + thick_shell
+#    return ellipsoid_ER(polar_outer, equat_outer)
 
 def random():
Index: sasmodels/models/core_shell_parallelepiped.c
===================================================================
--- sasmodels/models/core_shell_parallelepiped.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/core_shell_parallelepiped.c	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -25,4 +25,40 @@
         2.0 * length_a * length_b * thick_rim_c;
 #endif
+}
+
+static double
+radius_from_volume(double length_a, double length_b, double length_c,
+                   double thick_rim_a, double thick_rim_b, double thick_rim_c)
+{
+    const double volume_paral = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c);
+    return cbrt(0.75*volume_paral/M_PI);
+}
+
+static double
+radius_from_crosssection(double length_a, double length_b, double thick_rim_a, double thick_rim_b)
+{
+    const double area_xsec_paral = length_a*length_b + 2.0*thick_rim_a*length_b + 2.0*thick_rim_b*length_a;
+    return sqrt(area_xsec_paral/M_PI);
+}
+
+static double
+effective_radius(int mode, double length_a, double length_b, double length_c,
+                 double thick_rim_a, double thick_rim_b, double thick_rim_c)
+{
+    if (mode == 1) {
+        return radius_from_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c);
+    } else if (mode == 2) {
+        return 0.5 * (length_a + thick_rim_a);
+    } else if (mode == 3) {
+        return 0.5 * (length_b + thick_rim_b);
+    } else if (mode == 4) {
+        return 0.5 * (length_c + thick_rim_c);
+    } else if (mode == 5) {
+        return radius_from_crosssection(length_a, length_b, thick_rim_a, thick_rim_b);
+    } else if (mode == 6) {
+        return 0.5*sqrt(square(length_a+thick_rim_a) + square(length_b+thick_rim_b));
+    } else {
+        return 0.5*sqrt(square(length_a+thick_rim_a) + square(length_b+thick_rim_b) + square(length_c+thick_rim_c));
+    }
 }
 
Index: sasmodels/models/core_shell_parallelepiped.py
===================================================================
--- sasmodels/models/core_shell_parallelepiped.py	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/core_shell_parallelepiped.py	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -227,16 +227,17 @@
 source = ["lib/gauss76.c", "core_shell_parallelepiped.c"]
 have_Fq = True
-
-
-def ER(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c):
-    """
-        Return equivalent radius (ER)
-    """
-    from .parallelepiped import ER as ER_p
-
-    a = length_a + 2*thick_rim_a
-    b = length_b + 2*thick_rim_b
-    c = length_c + 2*thick_rim_c
-    return ER_p(a, b, c)
+effective_radius_type = ["equivalent sphere","half outer length_a", "half outer length_b", "half outer length_c",
+                         "equivalent circular cross-section","half outer ab diagonal","half outer diagonal"]
+
+#def ER(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c):
+#    """
+#        Return equivalent radius (ER)
+#    """
+#    from .parallelepiped import ER as ER_p
+#
+#    a = length_a + 2*thick_rim_a
+#    b = length_b + 2*thick_rim_b
+#    c = length_c + 2*thick_rim_c
+#    return ER_p(a, b, c)
 
 # VR defaults to 1.0
Index: sasmodels/models/core_shell_sphere.c
===================================================================
--- sasmodels/models/core_shell_sphere.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/core_shell_sphere.c	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -3,4 +3,14 @@
 {
     return M_4PI_3 * cube(radius + thickness);
+}
+
+static double
+effective_radius(int mode, double radius, double thickness)
+{
+    if (mode == 1) {
+        return radius + thickness;
+    } else {
+        return radius;
+    }
 }
 
Index: sasmodels/models/core_shell_sphere.py
===================================================================
--- sasmodels/models/core_shell_sphere.py	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/core_shell_sphere.py	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -78,15 +78,16 @@
 source = ["lib/sas_3j1x_x.c", "lib/core_shell.c", "core_shell_sphere.c"]
 have_Fq = True
+effective_radius_type = ["outer radius", "core radius"]
 
 demo = dict(scale=1, background=0, radius=60, thickness=10,
             sld_core=1.0, sld_shell=2.0, sld_solvent=0.0)
 
-def ER(radius, thickness):
-    """
-        Equivalent radius
-        @param radius: core radius
-        @param thickness: shell thickness
-    """
-    return radius + thickness
+#def ER(radius, thickness):
+#    """
+#        Equivalent radius
+#        @param radius: core radius
+#        @param thickness: shell thickness
+#    """
+#    return radius + thickness
 
 def VR(radius, thickness):
Index: sasmodels/models/cylinder.c
===================================================================
--- sasmodels/models/cylinder.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/cylinder.c	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -11,4 +11,34 @@
 {
     return sas_2J1x_x(qab*radius) * sas_sinx_x(qc*0.5*length);
+}
+
+static double
+radius_from_volume(double radius, double length)
+{
+    return cbrt(0.75*radius*radius*length);
+}
+
+static double
+radius_from_diagonal(double radius, double length)
+{
+    return sqrt(radius*radius + 0.25*length*length);
+}
+
+static double
+effective_radius(int mode, double radius, double length)
+{
+    if (mode == 1) {
+        return radius_from_volume(radius, length);
+    } else if (mode == 2) {
+        return radius;
+    } else if (mode == 3) {
+        return 0.5*length;
+    } else if (mode == 4) {
+        return (radius < 0.5*length ? radius : 0.5*length);
+    } else if (mode == 5) {
+        return (radius > 0.5*length ? radius : 0.5*length);
+    } else {
+        return radius_from_diagonal(radius,length);
+    }
 }
 
Index: sasmodels/models/cylinder.py
===================================================================
--- sasmodels/models/cylinder.py	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/cylinder.py	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -139,11 +139,12 @@
 source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "cylinder.c"]
 have_Fq = True
+effective_radius_type = ["equivalent sphere","radius","half length","half min dimension","half max dimension","half diagonal"]
 
-def ER(radius, length):
-    """
-        Return equivalent radius (ER)
-    """
-    ddd = 0.75 * radius * (2 * radius * length + (length + radius) * (length + pi * radius))
-    return 0.5 * (ddd) ** (1. / 3.)
+#def ER(radius, length):
+#    """
+#        Return equivalent radius (ER)
+#    """
+#    ddd = 0.75 * radius * (2 * radius * length + (length + radius) * (length + pi * radius))
+#    return 0.5 * (ddd) ** (1. / 3.)
 
 def random():
Index: sasmodels/models/ellipsoid.c
===================================================================
--- sasmodels/models/ellipsoid.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/ellipsoid.c	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -4,4 +4,45 @@
     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_volume(radius_polar, radius_equatorial);
+    } else if (mode == 2) {
+        return radius_from_curvature(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 01688448fd87fd51ac83d2d724652266c645773a)
+++ sasmodels/models/ellipsoid.py	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -169,26 +169,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
-    valid = (radius_polar * radius_equatorial != 0) & (radius_polar != radius_equatorial)
-    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 = 2.0 * (delta + 1.0) * (radius_polar * radius_equatorial**2)[valid]
-
-    r = np.zeros_like(radius_polar)
-    r[valid] = 0.5 * cbrt(ddd)
-    idx = radius_polar == radius_equatorial
-    r[idx] = radius_polar[idx]
-    return r
+effective_radius_type = ["equivalent sphere","average curvature", "min radius", "max radius"]
 
 def random():
Index: sasmodels/models/elliptical_cylinder.c
===================================================================
--- sasmodels/models/elliptical_cylinder.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/elliptical_cylinder.c	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -3,4 +3,56 @@
 {
     return M_PI * radius_minor * radius_minor * r_ratio * length;
+}
+
+static double
+radius_from_volume(double radius_minor, double r_ratio, double length)
+{
+    const double volume_ellcyl = form_volume(radius_minor,r_ratio,length);
+    return cbrt(0.75*volume_ellcyl/M_PI);
+}
+
+static double
+radius_from_min_dimension(double radius_minor, double r_ratio, double length)
+{
+    const double rad_min = (r_ratio > 1.0 ? radius_minor : r_ratio*radius_minor);
+    return (rad_min < length ? rad_min : length);
+}
+
+static double
+radius_from_max_dimension(double radius_minor, double r_ratio, double length)
+{
+    const double rad_max = (r_ratio < 1.0 ? radius_minor : r_ratio*radius_minor);
+    return (rad_max > length ? rad_max : length);
+}
+
+static double
+radius_from_diagonal(double radius_minor, double r_ratio, double length)
+{
+    const double radius_max = (r_ratio > 1.0 ? radius_minor*r_ratio : radius_minor);
+    return sqrt(radius_max*radius_max + 0.25*length*length);
+}
+
+static double
+effective_radius(int mode, double radius_minor, double r_ratio, double length)
+{
+    if (mode == 1) {
+        return radius_from_volume(radius_minor, r_ratio, length);
+    } else if (mode == 2) {
+        return 0.5*radius_minor*(1.0 + r_ratio);
+    } else if (mode == 3) {
+        return (r_ratio > 1.0 ? radius_minor : r_ratio*radius_minor);
+    } else if (mode == 4) {
+        return (r_ratio < 1.0 ? radius_minor : r_ratio*radius_minor);
+    } else if (mode == 5) {
+        return sqrt(radius_minor*radius_minor*r_ratio);
+    } else if (mode == 6) {
+        return 0.5*length;
+    } else if (mode == 7) {
+        return radius_from_min_dimension(radius_minor,r_ratio,length);
+    } else if (mode == 8) {
+        return radius_from_max_dimension(radius_minor,r_ratio,length);
+    } else {
+        return radius_from_diagonal(radius_minor,r_ratio,length);
+    }
 }
 
Index: sasmodels/models/elliptical_cylinder.py
===================================================================
--- sasmodels/models/elliptical_cylinder.py	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/elliptical_cylinder.py	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -123,4 +123,6 @@
 source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "elliptical_cylinder.c"]
 have_Fq = True
+effective_radius_type = ["equivalent sphere","average radius","min radius","max radius",
+                         "equivalent circular cross-section","half length","half min dimension","half max dimension","half diagonal"]
 
 demo = dict(scale=1, background=0, radius_minor=100, axis_ratio=1.5, length=400.0,
@@ -128,15 +130,15 @@
             theta_pd=10, phi_pd=2, psi_pd=3)
 
-def ER(radius_minor, axis_ratio, length):
-    """
-        Equivalent radius
-        @param radius_minor: Ellipse minor radius
-        @param axis_ratio: Ratio of major radius over minor radius
-        @param length: Length of the cylinder
-    """
-    radius = sqrt(radius_minor * radius_minor * axis_ratio)
-    ddd = 0.75 * radius * (2 * radius * length
-                           + (length + radius) * (length + pi * radius))
-    return 0.5 * (ddd) ** (1. / 3.)
+#def ER(radius_minor, axis_ratio, length):
+#    """
+#        Equivalent radius
+#        @param radius_minor: Ellipse minor radius
+#        @param axis_ratio: Ratio of major radius over minor radius
+#        @param length: Length of the cylinder
+#    """
+#    radius = sqrt(radius_minor * radius_minor * axis_ratio)
+#    ddd = 0.75 * radius * (2 * radius * length
+#                           + (length + radius) * (length + pi * radius))
+#    return 0.5 * (ddd) ** (1. / 3.)
 
 def random():
Index: sasmodels/models/fuzzy_sphere.c
===================================================================
--- sasmodels/models/fuzzy_sphere.c	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
+++ sasmodels/models/fuzzy_sphere.c	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -0,0 +1,26 @@
+static double form_volume(double radius, double fuzziness)
+{
+    return M_4PI_3*cube(radius);
+}
+
+static double
+effective_radius(int mode, double radius, double fuzziness)
+{
+    if (mode == 1) {
+        return radius;
+    } else {
+        return radius + fuzziness;
+    }
+}
+
+static void Fq(double q, double *F1, double *F2, double sld, double sld_solvent,
+               double radius, double fuzziness)
+{
+    const double qr = q*radius;
+    const double bes = sas_3j1x_x(qr);
+    const double qf = exp(-0.5*square(q*fuzziness));
+    const double contrast = (sld - sld_solvent);
+    const double form = contrast * form_volume(radius,fuzziness) * bes * qf;
+    *F1 = 1.0e-2*form;
+    *F2 = 1.0e-4*form*form;
+}
Index: sasmodels/models/fuzzy_sphere.py
===================================================================
--- sasmodels/models/fuzzy_sphere.py	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/fuzzy_sphere.py	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -78,35 +78,17 @@
               ["sld_solvent", "1e-6/Ang^2",  3, [-inf, inf], "sld",    "Solvent scattering length density"],
               ["radius",      "Ang",        60, [0, inf],    "volume", "Sphere radius"],
-              ["fuzziness",   "Ang",        10, [0, inf],    "",       "std deviation of Gaussian convolution for interface (must be << radius)"],
+              ["fuzziness",   "Ang",        10, [0, inf],    "volume",       "std deviation of Gaussian convolution for interface (must be << radius)"],
              ]
 # pylint: enable=bad-whitespace,line-too-long
 
-source = ["lib/sas_3j1x_x.c"]
+source = ["lib/sas_3j1x_x.c","fuzzy_sphere.c"]
 have_Fq = True
+effective_radius_type = ["radius","radius + fuzziness"]
 
-c_code = """
-static double form_volume(double radius)
-{
-    return M_4PI_3*cube(radius);
-}
-
-static void Fq(double q, double *F1, double *F2, double sld, double sld_solvent,
-               double radius, double fuzziness)
-{
-    const double qr = q*radius;
-    const double bes = sas_3j1x_x(qr);
-    const double qf = exp(-0.5*square(q*fuzziness));
-    const double contrast = (sld - sld_solvent);
-    const double form = contrast * form_volume(radius) * bes * qf;
-    *F1 = 1.0e-2*form;
-    *F2 = 1.0e-4*form*form;
-}
-"""
-
-def ER(radius):
-    """
-    Return radius
-    """
-    return radius
+#def ER(radius):
+#    """
+#    Return radius
+#    """
+#    return radius
 
 # VR defaults to 1.0
Index: sasmodels/models/hollow_cylinder.c
===================================================================
--- sasmodels/models/hollow_cylinder.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/hollow_cylinder.c	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -22,4 +22,34 @@
 }
 
+static double
+radius_from_volume(double radius, double thickness, double length)
+{
+    const double volume_outer_cyl = M_PI*square(radius + thickness)*length;
+    return cbrt(0.75*volume_outer_cyl/M_PI);
+}
+
+static double
+radius_from_diagonal(double radius, double thickness, double length)
+{
+    return sqrt(square(radius + thickness) + 0.25*square(length));
+}
+
+static double
+effective_radius(int mode, double radius, double thickness, double length)
+{
+    if (mode == 1) {
+        return radius_from_volume(radius, thickness, length);
+    } else if (mode == 2) {
+        return radius + thickness;
+    } else if (mode == 3) {
+        return 0.5*length;
+    } else if (mode == 4) {
+        return (radius + thickness < 0.5*length ? radius + thickness : 0.5*length);
+    } else if (mode == 5) {
+        return (radius + thickness > 0.5*length ? radius + thickness : 0.5*length);
+    } else {
+        return radius_from_diagonal(radius,thickness,length);
+    }
+}
 
 static void
Index: sasmodels/models/hollow_cylinder.py
===================================================================
--- sasmodels/models/hollow_cylinder.py	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/hollow_cylinder.py	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -90,23 +90,25 @@
 source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "hollow_cylinder.c"]
 have_Fq = True
+effective_radius_type = ["equivalent sphere","outer radius","half length",
+                         "half outer min dimension","half outer max dimension","half outer diagonal"]
 
 # pylint: disable=W0613
-def ER(radius, thickness, length):
-    """
-    :param radius:      Cylinder core radius
-    :param thickness:   Cylinder wall thickness
-    :param length:      Cylinder length
-    :return:            Effective radius
-    """
-    router = radius + thickness
-    if router == 0 or length == 0:
-        return 0.0
-    len1 = router
-    len2 = length/2.0
-    term1 = len1*len1*2.0*len2/2.0
-    term2 = 1.0 + (len2/len1)*(1.0 + 1/len2/2.0)*(1.0 + pi*len1/len2/2.0)
-    ddd = 3.0*term1*term2
-    diam = pow(ddd, (1.0/3.0))
-    return diam
+#def ER(radius, thickness, length):
+#    """
+#    :param radius:      Cylinder core radius
+#    :param thickness:   Cylinder wall thickness
+#    :param length:      Cylinder length
+#    :return:            Effective radius
+#    """
+#    router = radius + thickness
+#    if router == 0 or length == 0:
+#        return 0.0
+#    len1 = router
+#    len2 = length/2.0
+#    term1 = len1*len1*2.0*len2/2.0
+#    term2 = 1.0 + (len2/len1)*(1.0 + 1/len2/2.0)*(1.0 + pi*len1/len2/2.0)
+#    ddd = 3.0*term1*term2
+#    diam = pow(ddd, (1.0/3.0))
+#    return diam
 
 def VR(radius, thickness, length):
Index: sasmodels/models/hollow_rectangular_prism.c
===================================================================
--- sasmodels/models/hollow_rectangular_prism.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/hollow_rectangular_prism.c	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -11,4 +11,24 @@
     double vol_shell = vol_total - vol_core;
     return vol_shell;
+}
+
+static double
+effective_radius(int mode, double length_a, double b2a_ratio, double c2a_ratio, double thickness)
+{
+    if (mode == 1) {
+        return cbrt(0.75*cube(length_a)*b2a_ratio*c2a_ratio/M_PI);
+    } else if (mode == 2) {
+        return 0.5 * length_a;
+    } else if (mode == 3) {
+        return 0.5 * length_a*b2a_ratio;
+    } else if (mode == 4) {
+        return 0.5 * length_a*c2a_ratio;
+    } else if (mode == 5) {
+        return length_a*sqrt(b2a_ratio/M_PI);
+    } else if (mode == 6) {
+        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio)));
+    } else {
+        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio) + square(c2a_ratio)));
+    }
 }
 
Index: sasmodels/models/hollow_rectangular_prism.py
===================================================================
--- sasmodels/models/hollow_rectangular_prism.py	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/hollow_rectangular_prism.py	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -143,18 +143,20 @@
 source = ["lib/gauss76.c", "hollow_rectangular_prism.c"]
 have_Fq = True
-
-def ER(length_a, b2a_ratio, c2a_ratio, thickness):
-    """
-    Return equivalent radius (ER)
-    thickness parameter not used
-    """
-    b_side = length_a * b2a_ratio
-    c_side = length_a * c2a_ratio
-
-    # surface average radius (rough approximation)
-    surf_rad = sqrt(length_a * b_side / pi)
-
-    ddd = 0.75 * surf_rad * (2 * surf_rad * c_side + (c_side + surf_rad) * (c_side + pi * surf_rad))
-    return 0.5 * (ddd) ** (1. / 3.)
+effective_radius_type = ["equivalent sphere","half length_a", "half length_b", "half length_c",
+                         "equivalent outer circular cross-section","half ab diagonal","half diagonal"]
+
+#def ER(length_a, b2a_ratio, c2a_ratio, thickness):
+#    """
+#    Return equivalent radius (ER)
+#    thickness parameter not used
+#    """
+#    b_side = length_a * b2a_ratio
+#    c_side = length_a * c2a_ratio
+#
+#    # surface average radius (rough approximation)
+#    surf_rad = sqrt(length_a * b_side / pi)
+#
+#    ddd = 0.75 * surf_rad * (2 * surf_rad * c_side + (c_side + surf_rad) * (c_side + pi * surf_rad))
+#    return 0.5 * (ddd) ** (1. / 3.)
 
 def VR(length_a, b2a_ratio, c2a_ratio, thickness):
Index: sasmodels/models/hollow_rectangular_prism_thin_walls.c
===================================================================
--- sasmodels/models/hollow_rectangular_prism_thin_walls.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/hollow_rectangular_prism_thin_walls.c	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -6,4 +6,24 @@
     double vol_shell = 2.0 * (length_a*length_b + length_a*length_c + length_b*length_c);
     return vol_shell;
+}
+
+static double
+effective_radius(int mode, double length_a, double b2a_ratio, double c2a_ratio)
+{
+    if (mode == 1) {
+        return cbrt(0.75*cube(length_a)*b2a_ratio*c2a_ratio/M_PI);
+    } else if (mode == 2) {
+        return 0.5 * length_a;
+    } else if (mode == 3) {
+        return 0.5 * length_a*b2a_ratio;
+    } else if (mode == 4) {
+        return 0.5 * length_a*c2a_ratio;
+    } else if (mode == 5) {
+        return length_a*sqrt(b2a_ratio/M_PI);
+    } else if (mode == 6) {
+        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio)));
+    } else {
+        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio) + square(c2a_ratio)));
+    }
 }
 
Index: sasmodels/models/hollow_rectangular_prism_thin_walls.py
===================================================================
--- sasmodels/models/hollow_rectangular_prism_thin_walls.py	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/hollow_rectangular_prism_thin_walls.py	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -103,17 +103,19 @@
 source = ["lib/gauss76.c", "hollow_rectangular_prism_thin_walls.c"]
 have_Fq = True
+effective_radius_type = ["equivalent sphere","half length_a", "half length_b", "half length_c",
+                         "equivalent outer circular cross-section","half ab diagonal","half diagonal"]
 
-def ER(length_a, b2a_ratio, c2a_ratio):
-    """
-        Return equivalent radius (ER)
-    """
-    b_side = length_a * b2a_ratio
-    c_side = length_a * c2a_ratio
-
-    # surface average radius (rough approximation)
-    surf_rad = sqrt(length_a * b_side / pi)
-
-    ddd = 0.75 * surf_rad * (2 * surf_rad * c_side + (c_side + surf_rad) * (c_side + pi * surf_rad))
-    return 0.5 * (ddd) ** (1. / 3.)
+#def ER(length_a, b2a_ratio, c2a_ratio):
+#    """
+#        Return equivalent radius (ER)
+#    """
+#    b_side = length_a * b2a_ratio
+#    c_side = length_a * c2a_ratio
+#
+#    # surface average radius (rough approximation)
+#    surf_rad = sqrt(length_a * b_side / pi)
+#
+#    ddd = 0.75 * surf_rad * (2 * surf_rad * c_side + (c_side + surf_rad) * (c_side + pi * surf_rad))
+#    return 0.5 * (ddd) ** (1. / 3.)
 
 def VR(length_a, b2a_ratio, c2a_ratio):
Index: sasmodels/models/mono_gauss_coil.c
===================================================================
--- sasmodels/models/mono_gauss_coil.c	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
+++ sasmodels/models/mono_gauss_coil.c	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -0,0 +1,31 @@
+static double form_volume(double rg)
+{
+    return 1.0;
+}
+
+static double
+effective_radius(int mode, double rg)
+{
+    if (mode == 1) {
+        return rg;
+    } else if (mode == 2) {
+        return 2.0*rg;
+    } else if (mode == 3) {
+        return 3.0*rg;
+    } else {
+        return sqrt(5.0/3.0)*rg;
+    }
+}
+
+double Iq(double q, double i_zero, double rg)
+{
+    const double uarg = square(q*rg);
+    const double inten;
+    if (q == 0) {
+        inten = i_zero;
+    } else {
+        inten = 2.0*i_zero * (exp(-uarg) + uarg - 1.0)/square(uarg);
+    }
+
+    return inten;
+}
Index: sasmodels/models/mono_gauss_coil.py
===================================================================
--- sasmodels/models/mono_gauss_coil.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
+++ sasmodels/models/mono_gauss_coil.py	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -54,5 +54,5 @@
 
 import numpy as np
-from numpy import inf, exp, errstate
+from numpy import inf
 
 name = "mono_gauss_coil"
@@ -69,18 +69,24 @@
 parameters = [
     ["i_zero", "1/cm", 70.0, [0.0, inf], "", "Intensity at q=0"],
-    ["rg", "Ang", 75.0, [0.0, inf], "", "Radius of gyration"],
+    ["rg", "Ang", 75.0, [0.0, inf], "volume", "Radius of gyration"],
     ]
+
+source = ["mono_gauss_coil.c"]
+have_Fq = False
+effective_radius_type = ["R_g","2R_g","3R_g","(5/3)^0.5*R_g"]
+
+
 # pylint: enable=bad-whitespace, line-too-long
 
-# NB: Scale and Background are implicit parameters on every model
-def Iq(q, i_zero, rg):
-    # pylint: disable = missing-docstring
-    z = (q * rg)**2
-
-    with errstate(invalid='ignore'):
-        inten = (i_zero * 2.0) * (exp(-z) + z - 1.0)/z**2
-        inten[q == 0] = i_zero
-    return inten
-Iq.vectorized = True # Iq accepts an array of q values
+## NB: Scale and Background are implicit parameters on every model
+#def Iq(q, i_zero, rg):
+#    # pylint: disable = missing-docstring
+#    z = (q * rg)**2
+#
+#    with errstate(invalid='ignore'):
+#        inten = (i_zero * 2.0) * (exp(-z) + z - 1.0)/z**2
+#        inten[q == 0] = i_zero
+#    return inten
+#Iq.vectorized = True # Iq accepts an array of q values
 
 def random():
Index: sasmodels/models/multilayer_vesicle.c
===================================================================
--- sasmodels/models/multilayer_vesicle.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/multilayer_vesicle.c	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -47,4 +47,10 @@
 }
 
+static double
+effective_radius(int mode, double radius, double thick_shell, double thick_solvent, double fp_n_shells)
+{
+    return radius + fp_n_shells*thick_shell + (fp_n_shells - 1.0)*thick_solvent;
+}
+
 static void
 Fq(double q,
Index: sasmodels/models/multilayer_vesicle.py
===================================================================
--- sasmodels/models/multilayer_vesicle.py	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/multilayer_vesicle.py	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -146,8 +146,9 @@
 source = ["lib/sas_3j1x_x.c", "multilayer_vesicle.c"]
 have_Fq = True
+effective_radius_type = ["outer radius"]
 
-def ER(radius, thick_shell, thick_solvent, n_shells):
-    n_shells = int(n_shells+0.5)
-    return radius + n_shells * (thick_shell + thick_solvent) - thick_solvent
+#def ER(radius, thick_shell, thick_solvent, n_shells):
+#    n_shells = int(n_shells+0.5)
+#    return radius + n_shells * (thick_shell + thick_solvent) - thick_solvent
 
 def random():
Index: sasmodels/models/onion.c
===================================================================
--- sasmodels/models/onion.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/onion.c	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -40,4 +40,15 @@
 }
 
+static double
+effective_radius(int mode, double radius_core, double n_shells, double thickness[])
+{
+  int n = (int)(n_shells+0.5);
+  double r = radius_core;
+  for (int i=0; i < n; i++) {
+    r += thickness[i];
+  }
+  return r;
+}
+
 static void
 Fq(double q, double *F1, double *F2, double sld_core, double radius_core, double sld_solvent,
Index: sasmodels/models/onion.py
===================================================================
--- sasmodels/models/onion.py	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/onion.py	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -316,4 +316,5 @@
 single = False
 have_Fq = True
+effective_radius_type = ["outer radius"]
 
 profile_axes = ['Radius (A)', 'SLD (1e-6/A^2)']
@@ -366,8 +367,8 @@
     return np.asarray(z), np.asarray(rho)
 
-def ER(radius_core, n_shells, thickness):
-    """Effective radius"""
-    n = int(n_shells[0]+0.5)
-    return np.sum(thickness[:n], axis=0) + radius_core
+#def ER(radius_core, n_shells, thickness):
+#    """Effective radius"""
+#    n = int(n_shells[0]+0.5)
+#    return np.sum(thickness[:n], axis=0) + radius_core
 
 demo = {
Index: sasmodels/models/parallelepiped.c
===================================================================
--- sasmodels/models/parallelepiped.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/parallelepiped.c	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -3,4 +3,24 @@
 {
     return length_a * length_b * length_c;
+}
+
+static double
+effective_radius(int mode, double length_a, double length_b, double length_c)
+{
+    if (mode == 1) {
+        return cbrt(0.75*length_a*length_b*length_c/M_PI);
+    } else if (mode == 2) {
+        return 0.5 * length_a;
+    } else if (mode == 3) {
+        return 0.5 * length_b;
+    } else if (mode == 4) {
+        return 0.5 * length_c;
+    } else if (mode == 5) {
+        return sqrt(length_a*length_b/M_PI);
+    } else if (mode == 6) {
+        return 0.5*sqrt(length_a*length_a + length_b*length_b);
+    } else {
+        return 0.5*sqrt(length_a*length_a + length_b*length_b + length_c*length_c);
+    }
 }
 
Index: sasmodels/models/parallelepiped.py
===================================================================
--- sasmodels/models/parallelepiped.py	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/parallelepiped.py	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -231,20 +231,22 @@
 source = ["lib/gauss76.c", "parallelepiped.c"]
 have_Fq = True
-
-def ER(length_a, length_b, length_c):
-    """
-    Return effective radius (ER) for P(q)*S(q)
-    """
-    # now that axes can be in any size order, need to sort a,b,c
-    # where a~b and c is either much smaller or much larger
-    abc = np.vstack((length_a, length_b, length_c))
-    abc = np.sort(abc, axis=0)
-    selector = (abc[1] - abc[0]) > (abc[2] - abc[1])
-    length = np.where(selector, abc[0], abc[2])
-    # surface average radius (rough approximation)
-    radius = sqrt(np.where(~selector, abc[0]*abc[1], abc[1]*abc[2]) / pi)
-
-    ddd = 0.75 * radius * (2*radius*length + (length + radius)*(length + pi*radius))
-    return 0.5 * (ddd) ** (1. / 3.)
+effective_radius_type = ["equivalent sphere","half length_a", "half length_b", "half length_c",
+                         "equivalent circular cross-section","half ab diagonal","half diagonal"]
+
+#def ER(length_a, length_b, length_c):
+#    """
+#    Return effective radius (ER) for P(q)*S(q)
+#    """
+#    # now that axes can be in any size order, need to sort a,b,c
+#    # where a~b and c is either much smaller or much larger
+#    abc = np.vstack((length_a, length_b, length_c))
+#    abc = np.sort(abc, axis=0)
+#    selector = (abc[1] - abc[0]) > (abc[2] - abc[1])
+#    length = np.where(selector, abc[0], abc[2])
+#    # surface average radius (rough approximation)
+#    radius = sqrt(np.where(~selector, abc[0]*abc[1], abc[1]*abc[2]) / pi)
+#
+#    ddd = 0.75 * radius * (2*radius*length + (length + radius)*(length + pi*radius))
+#    return 0.5 * (ddd) ** (1. / 3.)
 
 # VR defaults to 1.0
Index: sasmodels/models/pringle.c
===================================================================
--- sasmodels/models/pringle.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
+++ sasmodels/models/pringle.c	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -104,4 +104,14 @@
 }
 
+static double
+effective_radius(int mode, double radius, double thickness, double alpha, double beta)
+{
+    if (mode == 1) {
+        return cbrt(0.75*radius*radius*thickness);
+    } else {
+        return radius;
+    }
+}
+
 double Iq(
     double q,
Index: sasmodels/models/pringle.py
===================================================================
--- sasmodels/models/pringle.py	(revision ef07e950965796e2a3a8156fcade98ce5e403001)
+++ sasmodels/models/pringle.py	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -74,12 +74,13 @@
 source = ["lib/polevl.c", "lib/sas_J0.c", "lib/sas_J1.c", \
           "lib/sas_JN.c", "lib/gauss76.c", "pringle.c"]
+effective_radius_type = ["equivalent sphere","radius"]
 
-def ER(radius, thickness, alpha, beta):
-    """
-    Return equivalent radius (ER)
-    """
-    ddd = 0.75 * radius * (2 * radius * thickness + (thickness + radius) \
-                           * (thickness + pi * radius))
-    return 0.5 * (ddd) ** (1. / 3.)
+#def ER(radius, thickness, alpha, beta):
+#    """
+#    Return equivalent radius (ER)
+#    """
+#    ddd = 0.75 * radius * (2 * radius * thickness + (thickness + radius) \
+#                           * (thickness + pi * radius))
+#    return 0.5 * (ddd) ** (1. / 3.)
 
 def random():
Index: sasmodels/models/raspberry.c
===================================================================
--- sasmodels/models/raspberry.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/raspberry.c	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -1,3 +1,3 @@
-double form_volume(double radius_lg);
+double form_volume(double radius_lg, double radius_sm, double penetration);
 
 double Iq(double q,
@@ -6,5 +6,5 @@
           double radius_lg, double radius_sm, double penetration);
 
-double form_volume(double radius_lg)
+double form_volume(double radius_lg, double radius_sm, double penetration)
 {
     //Because of the complex structure, volume normalization must
@@ -12,4 +12,14 @@
     double volume=1.0;
     return volume;
+}
+
+static double
+effective_radius(int mode, double radius_lg, double radius_sm, double penetration)
+{
+    if (mode == 1) {
+        return radius_lg;
+    } else {
+        return radius_lg + 2.0*radius_sm - penetration;
+    }
 }
 
Index: sasmodels/models/raspberry.py
===================================================================
--- sasmodels/models/raspberry.py	(revision ef07e950965796e2a3a8156fcade98ce5e403001)
+++ sasmodels/models/raspberry.py	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -145,11 +145,12 @@
               ["radius_lg", "Ang", 5000, [0, inf], "volume",
                "radius of large spheres"],
-              ["radius_sm", "Ang", 100, [0, inf], "",
+              ["radius_sm", "Ang", 100, [0, inf], "volume",
                "radius of small spheres"],
-              ["penetration", "Ang", 0, [-1, 1], "",
+              ["penetration", "Ang", 0, [-1, 1], "volume",
                "fractional penetration depth of small spheres into large sphere"],
              ]
 
 source = ["lib/sas_3j1x_x.c", "raspberry.c"]
+effective_radius_type = ["radius_large","radius_outer"]
 
 def random():
Index: sasmodels/models/rectangular_prism.c
===================================================================
--- sasmodels/models/rectangular_prism.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/rectangular_prism.c	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -3,4 +3,24 @@
 {
     return length_a * (length_a*b2a_ratio) * (length_a*c2a_ratio);
+}
+
+static double
+effective_radius(int mode, double length_a, double b2a_ratio, double c2a_ratio)
+{
+    if (mode == 1) {
+        return cbrt(0.75*cube(length_a)*b2a_ratio*c2a_ratio/M_PI);
+    } else if (mode == 2) {
+        return 0.5 * length_a;
+    } else if (mode == 3) {
+        return 0.5 * length_a*b2a_ratio;
+    } else if (mode == 4) {
+        return 0.5 * length_a*c2a_ratio;
+    } else if (mode == 5) {
+        return length_a*sqrt(b2a_ratio/M_PI);
+    } else if (mode == 6) {
+        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio)));
+    } else {
+        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio) + square(c2a_ratio)));
+    }
 }
 
Index: sasmodels/models/rectangular_prism.py
===================================================================
--- sasmodels/models/rectangular_prism.py	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/rectangular_prism.py	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -136,17 +136,19 @@
 source = ["lib/gauss76.c", "rectangular_prism.c"]
 have_Fq = True
+effective_radius_type = ["equivalent sphere","half length_a", "half length_b", "half length_c",
+                         "equivalent circular cross-section","half ab diagonal","half diagonal"]
 
-def ER(length_a, b2a_ratio, c2a_ratio):
-    """
-        Return equivalent radius (ER)
-    """
-    b_side = length_a * b2a_ratio
-    c_side = length_a * c2a_ratio
-
-    # surface average radius (rough approximation)
-    surf_rad = sqrt(length_a * b_side / pi)
-
-    ddd = 0.75 * surf_rad * (2 * surf_rad * c_side + (c_side + surf_rad) * (c_side + pi * surf_rad))
-    return 0.5 * (ddd) ** (1. / 3.)
+#def ER(length_a, b2a_ratio, c2a_ratio):
+#    """
+#        Return equivalent radius (ER)
+#    """
+#    b_side = length_a * b2a_ratio
+#    c_side = length_a * c2a_ratio
+#
+#    # surface average radius (rough approximation)
+#    surf_rad = sqrt(length_a * b_side / pi)
+#
+#    ddd = 0.75 * surf_rad * (2 * surf_rad * c_side + (c_side + surf_rad) * (c_side + pi * surf_rad))
+#    return 0.5 * (ddd) ** (1. / 3.)
 
 def random():
Index: sasmodels/models/sphere.c
===================================================================
--- sasmodels/models/sphere.c	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
+++ sasmodels/models/sphere.c	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -0,0 +1,19 @@
+static double form_volume(double radius)
+{
+    return M_4PI_3*cube(radius);
+}
+
+static double
+effective_radius(int mode, double radius)
+{
+    return radius;
+}
+
+static void Fq(double q, double *f1, double *f2, double sld, double sld_solvent, double radius)
+{
+    const double bes = sas_3j1x_x(q*radius);
+    const double contrast = (sld - sld_solvent);
+    const double form = contrast * form_volume(radius) * bes;
+    *f1 = 1.0e-2*form;
+    *f2 = 1.0e-4*form*form;
+}
Index: sasmodels/models/sphere.py
===================================================================
--- sasmodels/models/sphere.py	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/sphere.py	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -67,28 +67,13 @@
              ]
 
-source = ["lib/sas_3j1x_x.c"]
+source = ["lib/sas_3j1x_x.c","sphere.c"]
 have_Fq = True
+effective_radius_type = ["radius"]
 
-c_code = """
-static double form_volume(double radius)
-{
-    return M_4PI_3*cube(radius);
-}
-
-static void Fq(double q, double *f1, double *f2, double sld, double sld_solvent, double radius)
-{
-    const double bes = sas_3j1x_x(q*radius);
-    const double contrast = (sld - sld_solvent);
-    const double form = contrast * form_volume(radius) * bes;
-    *f1 = 1.0e-2*form;
-    *f2 = 1.0e-4*form*form;
-}
-"""
-
-def ER(radius):
-    """
-    Return equivalent radius (ER)
-    """
-    return radius
+#def ER(radius):
+#    """
+#    Return equivalent radius (ER)
+#    """
+#    return radius
 
 # VR defaults to 1.0
Index: sasmodels/models/spherical_sld.c
===================================================================
--- sasmodels/models/spherical_sld.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/spherical_sld.c	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -10,4 +10,15 @@
     }
     return M_4PI_3*cube(r);
+}
+
+static double
+effective_radius(int mode, double fp_n_shells, double thickness[], double interface[])
+{
+    int n_shells= (int)(fp_n_shells + 0.5);
+    double r = 0.0;
+    for (int i=0; i < n_shells; i++) {
+        r += thickness[i] + interface[i];
+    }
+    return r;
 }
 
Index: sasmodels/models/spherical_sld.py
===================================================================
--- sasmodels/models/spherical_sld.py	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/spherical_sld.py	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -210,4 +210,5 @@
 single = False  # TODO: fix low q behaviour
 have_Fq = True
+effective_radius_type = ["outer radius"]
 
 profile_axes = ['Radius (A)', 'SLD (1e-6/A^2)']
@@ -255,10 +256,10 @@
 
 
-def ER(n_shells, thickness, interface):
-    """Effective radius"""
-    n_shells = int(n_shells + 0.5)
-    total = (np.sum(thickness[:n_shells], axis=1)
-             + np.sum(interface[:n_shells], axis=1))
-    return total
+#def ER(n_shells, thickness, interface):
+#    """Effective radius"""
+#    n_shells = int(n_shells + 0.5)
+#    total = (np.sum(thickness[:n_shells], axis=1)
+#             + np.sum(interface[:n_shells], axis=1))
+#    return total
 
 
Index: sasmodels/models/triaxial_ellipsoid.c
===================================================================
--- sasmodels/models/triaxial_ellipsoid.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/triaxial_ellipsoid.c	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -7,4 +7,35 @@
 }
 
+static double
+radius_from_volume(double radius_equat_minor, double radius_equat_major, double radius_polar)
+{
+    return cbrt(radius_equat_minor*radius_equat_major*radius_polar);
+}
+
+static double
+radius_from_min_dimension(double radius_equat_minor, double radius_equat_major, double radius_polar)
+{
+    const double rad_equat_min = (radius_equat_minor < radius_equat_major ? radius_equat_minor : radius_equat_major);
+    return (rad_equat_min < radius_polar ? rad_equat_min : radius_polar);
+}
+
+static double
+radius_from_max_dimension(double radius_equat_minor, double radius_equat_major, double radius_polar)
+{
+    const double rad_equat_max = (radius_equat_minor < radius_equat_major ? radius_equat_major : radius_equat_minor);
+    return (rad_equat_max > radius_polar ? rad_equat_max : radius_polar);
+}
+
+static double
+effective_radius(int mode, double radius_equat_minor, double radius_equat_major, double radius_polar)
+{
+    if (mode == 1) {
+        return radius_from_volume(radius_equat_minor,radius_equat_major, radius_polar);
+    } else if (mode == 2) {
+        return radius_from_min_dimension(radius_equat_minor,radius_equat_major, radius_polar);
+    } else {
+        return radius_from_max_dimension(radius_equat_minor,radius_equat_major, radius_polar);
+    }
+}
 
 static void
Index: sasmodels/models/triaxial_ellipsoid.py
===================================================================
--- sasmodels/models/triaxial_ellipsoid.py	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/triaxial_ellipsoid.py	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -158,4 +158,5 @@
 source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "triaxial_ellipsoid.c"]
 have_Fq = True
+effective_radius_type = ["equivalent sphere","min radius", "max radius"]
 
 def ER(radius_equat_minor, radius_equat_major, radius_polar):
Index: sasmodels/models/vesicle.c
===================================================================
--- sasmodels/models/vesicle.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/vesicle.c	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -6,4 +6,9 @@
 }
 
+static double
+effective_radius(int mode, double radius, double thickness)
+{
+    return radius + thickness;
+}
 
 static void
Index: sasmodels/models/vesicle.py
===================================================================
--- sasmodels/models/vesicle.py	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/vesicle.py	(revision d277229f02cb7770fbb4ea2ad506d96bed705d75)
@@ -95,13 +95,14 @@
 source = ["lib/sas_3j1x_x.c", "vesicle.c"]
 have_Fq = True
+effective_radius_type = ["outer radius"]
 
-def ER(radius, thickness):
-    '''
-    returns the effective radius used in the S*P calculation
-
-    :param radius: core radius
-    :param thickness: shell thickness
-    '''
-    return radius + thickness
+#def ER(radius, thickness):
+#    '''
+#    returns the effective radius used in the S*P calculation
+#
+#    :param radius: core radius
+#    :param thickness: shell thickness
+#    '''
+#    return radius + thickness
 
 def VR(radius, thickness):
Index: sasmodels/product.py
===================================================================
--- sasmodels/product.py	(revision c0131d44f4e88fca400e52fffee1c461d7b4662d)
+++ sasmodels/product.py	(revision 33d7be33b92c42bb96a5c8ef527d96896321b221)
@@ -24,5 +24,5 @@
 # pylint: disable=unused-import
 try:
-    from typing import Tuple, Callable
+    from typing import Tuple, Callable, Union
 except ImportError:
     pass
@@ -37,7 +37,6 @@
 #]
 
-ER_ID = "radius_effective"
-VF_ID = "volfraction"
-
+RADIUS_ID = "radius_effective"
+VOLFRAC_ID = "volfraction"
 def make_extra_pars(p_info):
     pars = []
@@ -51,4 +50,13 @@
                 "Structure factor calculation")
         pars.append(par)
+    if p_info.effective_radius_type is not None:
+        par = parse_parameter(
+                "radius_effective_mode",
+                "",
+                0,
+                [["unconstrained"] + p_info.effective_radius_type],
+                "",
+                "Effective radius calculation")
+        pars.append(par)
     return pars
 
@@ -65,9 +73,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")
@@ -75,11 +83,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.
@@ -87,6 +94,10 @@
         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))
+                          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)
@@ -94,5 +105,5 @@
     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()
@@ -157,6 +168,6 @@
     return par
 
-def _intermediates(P, S):
-    # type: (np.ndarray, np.ndarray) -> OrderedDict[str, np.ndarray]
+def _intermediates(P, S, effective_radius):
+    # type: (np.ndarray, np.ndarray, float) -> OrderedDict[str, np.ndarray]
     """
     Returns intermediate results for standard product (P(Q)*S(Q))
@@ -165,10 +176,17 @@
         ("P(Q)", P),
         ("S(Q)", S),
+        ("effective_radius", effective_radius),
     ))
 
-def _intermediates_beta(F1, F2, S, scale, bg):
-    # type: (np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray) -> OrderedDict[str, np.ndarray]
-    """
-    Returns intermediate results for beta approximation-enabled product
+def _intermediates_beta(F1,              # type: np.ndarray
+                        F2,              # type: np.ndarray
+                        S,               # type: np.ndarray
+                        scale,           # type: np.ndarray
+                        bg,              # type: np.ndarray
+                        effective_radius # type: float
+                        ):
+    # type: (...) -> OrderedDict[str, Union[np.ndarray, float]]
+    """
+    Returns intermediate results for beta approximation-enabled product. The result may be an array or a float.
     """
     # TODO: 1. include calculated Q vector
@@ -179,4 +197,5 @@
         ("beta(Q)", F1**2 / F2),
         ("S_eff(Q)", 1 + (F1**2 / F2)*(S-1)),
+        ("effective_radius", effective_radius),
         # ("I(Q)", scale*(F2 + (F1**2)*(S-1)) + bg),
     ))
@@ -262,17 +281,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]
@@ -280,16 +288,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
@@ -298,21 +299,39 @@
 
         # 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 = int(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:
+                # Plug R_eff from p into S model (initial value and pd value)
+                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
 
-            self.results = lambda: _intermediates_beta(F1, F2, s_result, volfrac/volume_avg, background)
+            self.results = lambda: _intermediates_beta(F1, F2, s_result, volfrac/volume_avg, background, effective_radius)
             final_result = combined_scale*(F2 + (F1**2)*(s_result - 1)) + background
 
         else:
-            p_result = self.p_kernel.Iq(p_details, p_values, cutoff, magnetic)
-
-            self.results = lambda: _intermediates(p_result, s_result)
+            p_result, effective_radius = self.p_kernel.Pq_Reff(
+                p_details, p_values, cutoff, magnetic, effective_radius_type)
+            if effective_radius_type > 0:
+                # Plug R_eff from p into S model (initial value and pd value)
+                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 = lambda: _intermediates(p_result, s_result, effective_radius)
             final_result = scale*(p_result*s_result) + background
 
Index: sasmodels/sasview_model.py
===================================================================
--- sasmodels/sasview_model.py	(revision 84f29620067dee7dc1c9d0c004c375ea633f98a9)
+++ sasmodels/sasview_model.py	(revision 84f29620067dee7dc1c9d0c004c375ea633f98a9)
@@ -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
@@ -381,5 +381,8 @@
                 continue
             self.params[p.name] = p.default
-            self.details[p.id] = [p.units, p.limits[0], p.limits[1]]
+            if p.limits and type(p.limits) is list and len(p.limits) > 1:
+                self.details[p.id] = [p.units if not p.choices else p.choices, p.limits[0], p.limits[1]]
+            else:
+                self.details[p.id] = [p.units if not p.choices else p.choices, None, None]
             if p.polydisperse:
                 self.details[p.id+".width"] = [
