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: explore/beta/sasfit_compare_new.py
===================================================================
--- explore/beta/sasfit_compare_new.py	(revision 3f818b2beab7c1e1a82656f61f49a19156bec0f7)
+++ explore/beta/sasfit_compare_new.py	(revision 3f818b2beab7c1e1a82656f61f49a19156bec0f7)
@@ -0,0 +1,913 @@
+from __future__ import division, print_function
+# Make sasmodels available on the path
+import sys, os
+BETA_DIR = os.path.dirname(os.path.realpath(__file__))
+SASMODELS_DIR = os.path.dirname(os.path.dirname(BETA_DIR))
+sys.path.insert(0, SASMODELS_DIR)
+
+from collections import namedtuple
+
+from matplotlib import pyplot as plt
+import numpy as np
+from numpy import pi, sin, cos, sqrt, fabs
+from numpy.polynomial.legendre import leggauss
+from scipy.special import j1 as J1
+from numpy import inf
+from scipy.special import gammaln  # type: ignore
+
+Theory = namedtuple('Theory', 'Q F1 F2 P S I Seff Ibeta')
+Theory.__new__.__defaults__ = (None,) * len(Theory._fields)
+
+#Used to calculate F(q) for the cylinder, sphere, ellipsoid models
+# RKH adding vesicle and hollow_cylinder to test sasview special cases of ER and VR
+def sas_sinx_x(x):
+    with np.errstate(all='ignore'):
+        retvalue = sin(x)/x
+    retvalue[x == 0.] = 1.
+    return retvalue
+
+def sas_2J1x_x(x):
+    with np.errstate(all='ignore'):
+        retvalue = 2*J1(x)/x
+    retvalue[x == 0] = 1.
+    return retvalue
+
+def sas_3j1x_x(x):
+    """return 3*j1(x)/x"""
+    retvalue = np.empty_like(x)
+    with np.errstate(all='ignore'):
+        # GSL bessel_j1 taylor expansion
+        index = (x < 0.25)
+        y = x[index]**2
+        c1 = -1.0/10.0
+        c2 = +1.0/280.0
+        c3 = -1.0/15120.0
+        c4 = +1.0/1330560.0
+        c5 = -1.0/172972800.0
+        retvalue[index] = 1.0 + y*(c1 + y*(c2 + y*(c3 + y*(c4 + y*c5))))
+        index = ~index
+        y = x[index]
+        retvalue[index] = 3*(sin(y) - y*cos(y))/y**3
+    retvalue[x == 0.] = 1.
+    return retvalue
+
+#Used to cross check my models with sasview models
+def build_model(model_name, q, **pars):
+    from sasmodels.core import load_model_info, build_model as build_sasmodel
+    from sasmodels.data import empty_data1D
+    from sasmodels.direct_model import DirectModel
+    model_info = load_model_info(model_name)
+    model = build_sasmodel(model_info, dtype='double!')
+    data = empty_data1D(q)
+    calculator = DirectModel(data, model,cutoff=0)
+    calculator.pars = pars.copy()
+    calculator.pars.setdefault('background', 0)
+    return calculator
+
+#gives the hardsphere structure factor that sasview uses
+def _hardsphere_simple(q, radius_effective, volfraction):
+    CUTOFFHS = 0.05
+    if fabs(radius_effective) < 1.E-12:
+        HARDSPH = 1.0
+        return HARDSPH
+    X = 1.0/(1.0 -volfraction)
+    D = X*X
+    A = (1.+2.*volfraction)*D
+    A *= A
+    X = fabs(q*radius_effective*2.0)
+    if X < 5.E-06:
+        HARDSPH = 1./A
+        return HARDSPH
+    X2 = X*X
+    B = (1.0 +0.5*volfraction)*D
+    B *= B
+    B *= -6.*volfraction
+    G = 0.5*volfraction*A
+    if X < CUTOFFHS:
+        FF = 8.0*A +6.0*B + 4.0*G + ( -0.8*A -B/1.5 -0.5*G +(A/35. +0.0125*B +0.02*G)*X2)*X2
+        HARDSPH = 1./(1. + volfraction*FF )
+        return HARDSPH
+    X4 = X2*X2
+    S, C = sin(X), cos(X)
+    FF =  ((G*( (4.*X2 -24.)*X*S -(X4 -12.*X2 +24.)*C +24. )/X2 + B*(2.*X*S -(X2-2.)*C -2.) )/X + A*(S-X*C))/X
+    HARDSPH = 1./(1. + 24.*volfraction*FF/X2)
+    return HARDSPH
+
+def hardsphere_simple(q, radius_effective, volfraction):
+    SQ = [_hardsphere_simple(qk, radius_effective, volfraction) for qk in q]
+    return np.array(SQ)
+
+#Used in gaussian quadrature for polydispersity
+#returns values and the probability of those values based on gaussian distribution
+N_GAUSS = 35
+NSIGMA_GAUSS = 3
+def gaussian_distribution(center, sigma, lb, ub):
+    #3 standard deviations covers approx. 99.7%
+    if sigma != 0:
+        nsigmas = NSIGMA_GAUSS
+        x = np.linspace(center-sigma*nsigmas, center+sigma*nsigmas, num=N_GAUSS)
+        x = x[(x >= lb) & (x <= ub)]
+        px = np.exp((x-center)**2 / (-2.0 * sigma * sigma))
+        return x, px
+    else:
+        return np.array([center]), np.array([1])
+
+N_SCHULZ = 80
+NSIGMA_SCHULZ = 8
+def schulz_distribution(center, sigma, lb, ub):
+    if sigma != 0:
+        nsigmas = NSIGMA_SCHULZ
+        x = np.linspace(center-sigma*nsigmas, center+sigma*nsigmas, num=N_SCHULZ)
+        x = x[(x >= lb) & (x <= ub)]
+        R = x/center
+        z = (center/sigma)**2
+        arg = z*np.log(z) + (z-1)*np.log(R) - R*z - np.log(center) - gammaln(z)
+        px = np.exp(arg)
+        return x, px
+    else:
+        return np.array([center]), np.array([1])
+
+#returns the effective radius used in sasview
+def ER_ellipsoid(radius_polar, radius_equatorial):
+    ee = np.empty_like(radius_polar)
+    if radius_polar > radius_equatorial:
+        ee = (radius_polar**2 - radius_equatorial**2)/radius_polar**2
+    elif radius_polar < radius_equatorial:
+        ee = (radius_equatorial**2 - radius_polar**2) / radius_equatorial**2
+    else:
+        ee = 2*radius_polar
+    if radius_polar * radius_equatorial != 0:
+        bd = 1.0 - ee
+        e1 = np.sqrt(ee)
+        b1 = 1.0 + np.arcsin(e1) / (e1*np.sqrt(bd))
+        bL = (1.0 + e1) / (1.0 - e1)
+        b2 = 1.0 + bd / 2 / e1 * np.log(bL)
+        delta = 0.75 * b1 * b2
+    ddd = np.zeros_like(radius_polar)
+    ddd = 2.0*(delta + 1.0)*radius_polar*radius_equatorial**2
+    return 0.5*ddd**(1.0 / 3.0)
+
+def ellipsoid_volume(radius_polar, radius_equatorial):
+    volume = (4./3.)*pi*radius_polar*radius_equatorial**2
+    return volume
+
+# F1 is F(q)
+# F2 is F(g)^2
+#IQM is I(q) with monodispersity
+#IQSM is I(q) with structure factor S(q) and monodispersity
+#IQBM is I(q) with Beta Approximation and monodispersity
+#SQ is monodisperse approach for structure factor
+#SQ_EFF is the effective structure factor from beta approx
+def ellipsoid_theta(q, radius_polar, radius_equatorial, sld, sld_solvent,
+                    volfraction=0, radius_effective=None):
+    #creates values z and corresponding probabilities w from legendre-gauss quadrature
+    volume = ellipsoid_volume(radius_polar, radius_equatorial)
+    z, w = leggauss(76)
+    F1 = np.zeros_like(q)
+    F2 = np.zeros_like(q)
+    #use a u subsition(u=cos) and then u=(z+1)/2 to change integration from
+    #0->2pi with respect to alpha to -1->1 with respect to z, allowing us to use
+    #legendre-gauss quadrature
+    for k, qk in enumerate(q):
+        r = sqrt(radius_equatorial**2*(1-((z+1)/2)**2)+radius_polar**2*((z+1)/2)**2)
+        form = (sld-sld_solvent)*volume*sas_3j1x_x(qk*r)
+        F2[k] = np.sum(w*form**2)
+        F1[k] = np.sum(w*form)
+    #the 1/2 comes from the change of variables mentioned above
+    F2 = F2/2.0
+    F1 = F1/2.0
+    if radius_effective is None:
+        radius_effective = ER_ellipsoid(radius_polar,radius_equatorial)
+    SQ = hardsphere_simple(q, radius_effective, volfraction)
+    SQ_EFF = 1 + F1**2/F2*(SQ - 1)
+    IQM = 1e-4*F2/volume
+    IQSM = volfraction*IQM*SQ
+    IQBM = volfraction*IQM*SQ_EFF
+    return Theory(Q=q, F1=F1, F2=F2, P=IQM, S=SQ, I=IQSM, Seff=SQ_EFF, Ibeta=IQBM)
+
+#IQD is I(q) polydispursed, IQSD is I(q)S(q) polydispursed, etc.
+#IQBD HAS NOT BEEN CROSS CHECKED AT ALL
+def ellipsoid_pe(q, radius_polar, radius_equatorial, sld, sld_solvent,
+                 radius_polar_pd=0.1, radius_equatorial_pd=0.1,
+                 radius_polar_pd_type='gaussian',
+                 radius_equatorial_pd_type='gaussian',
+                 volfraction=0, radius_effective=None,
+                 background=0, scale=1,
+                 norm='sasview'):
+    if norm not in ['sasview', 'sasfit', 'yun']:
+        raise TypeError("unknown norm "+norm)
+    if radius_polar_pd_type == 'gaussian':
+        Rp_val, Rp_prob = gaussian_distribution(radius_polar, radius_polar_pd*radius_polar, 0, inf)
+    elif radius_polar_pd_type == 'schulz':
+        Rp_val, Rp_prob = schulz_distribution(radius_polar, radius_polar_pd*radius_polar, 0, inf)
+    if radius_equatorial_pd_type == 'gaussian':
+        Re_val, Re_prob = gaussian_distribution(radius_equatorial, radius_equatorial_pd*radius_equatorial, 0, inf)
+    elif radius_equatorial_pd_type == 'schulz':
+        Re_val, Re_prob = schulz_distribution(radius_equatorial, radius_equatorial_pd*radius_equatorial, 0, inf)
+    total_weight = total_volume = 0
+    radius_eff = 0
+    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)
+            volume = ellipsoid_volume(Rpk, Rei)
+            weight = Rp_prob[k]*Re_prob[i]
+            total_weight += weight
+            total_volume += weight*volume
+            F1 += theory.F1*weight
+            F2 += theory.F2*weight
+            radius_eff += weight*ER_ellipsoid(Rpk, Rei)
+    F1 /= total_weight
+    F2 /= total_weight
+    average_volume = total_volume/total_weight
+    if radius_effective is None:
+        radius_effective = radius_eff/total_weight
+    if norm == 'sasfit':
+        IQD = F2
+    elif norm == 'sasview':
+        # Note: internally, sasview uses F2/total_volume because:
+        #   average_volume = total_volume/total_weight
+        #   F2/total_weight / average_volume
+        #     = F2/total_weight / total_volume/total_weight
+        #     = F2/total_volume
+        IQD = F2/average_volume*1e-4*volfraction
+        F1 *= 1e-2  # Yun is using sld in 1/A^2, not 1e-6/A^2
+        F2 *= 1e-4
+    elif norm == 'yun':
+        F1 *= 1e-6  # Yun is using sld in 1/A^2, not 1e-6/A^2
+        F2 *= 1e-12
+        IQD = F2/average_volume*1e8*volfraction
+    SQ = hardsphere_simple(q, radius_effective, volfraction)
+    beta = F1**2/F2
+    SQ_EFF = 1 + beta*(SQ - 1)
+    IQSD = IQD*SQ
+    IQBD = IQD*SQ_EFF
+    return Theory(Q=q, F1=F1, F2=F2, P=IQD, S=SQ, I=IQSD, Seff=SQ_EFF, Ibeta=IQBD)
+
+#polydispersity for sphere
+def sphere_r(q,radius,sld,sld_solvent,
+             radius_pd=0.1, radius_pd_type='gaussian',
+             volfraction=0, radius_effective=None,
+             background=0, scale=1,
+             norm='sasview'):
+    if norm not in ['sasview', 'sasfit', 'yun']:
+        raise TypeError("unknown norm "+norm)
+    if radius_pd_type == 'gaussian':
+        radius_val, radius_prob = gaussian_distribution(radius, radius_pd*radius, 0, inf)
+    elif radius_pd_type == 'schulz':
+        radius_val, radius_prob = schulz_distribution(radius, radius_pd*radius, 0, inf)
+    total_weight = total_volume = 0
+    F1 = np.zeros_like(q)
+    F2 = np.zeros_like(q)
+    for k, rk in enumerate(radius_val):
+        volume = 4./3.*pi*rk**3
+        total_weight += radius_prob[k]
+        total_volume += radius_prob[k]*volume
+        form = (sld-sld_solvent)*volume*sas_3j1x_x(q*rk)
+        F2 += radius_prob[k]*form**2
+        F1 += radius_prob[k]*form
+    F1 /= total_weight
+    F2 /= total_weight
+    average_volume = total_volume/total_weight
+
+    if radius_effective is None:
+        radius_effective = radius
+        print("radius_effective  ",radius_effective)
+        average_volume = total_volume/total_weight
+    if norm == 'sasfit':
+        IQD = F2
+    elif norm == 'sasview':
+        IQD = F2/average_volume*1e-4*volfraction
+    elif norm == 'yun':
+        F1 *= 1e-6  # Yun is using sld in 1/A^2, not 1e-6/A^2
+        F2 *= 1e-12
+        IQD = F2/average_volume*1e8*volfraction
+    SQ = hardsphere_simple(q, radius_effective, volfraction)
+    beta = F1**2/F2
+    SQ_EFF = 1 + beta*(SQ - 1)
+    IQSD = IQD*SQ
+    IQBD = IQD*SQ_EFF
+    return Theory(Q=q, F1=F1, F2=F2, P=IQD, S=SQ, I=IQSD, Seff=SQ_EFF, Ibeta=IQBD)
+
+#polydispersity for vesicle
+def vesicle_pe(q,radius=20, thickness=10, sld=4, sld_solvent=1, volfraction=0.03,
+        radius_pd=0.1, thickness_pd=0.2, radius_pd_type="gaussian", thickness_pd_type="gaussian",
+        radius_effective=None, background=0, scale=1,
+                 norm='sasview'):
+    if norm not in ['sasview', 'sasfit', 'yun']:
+        raise TypeError("unknown norm "+norm)
+    if radius_pd_type == 'gaussian':
+        R_val, R_prob = gaussian_distribution(radius, radius_pd*radius, 0, inf)
+    elif radius_pd_type == 'schulz':
+        R_val, R_prob = schulz_distribution(radius, radius_pd*radius, 0, inf)
+    if thickness_pd_type == 'gaussian':
+        T_val, T_prob = gaussian_distribution(thickness, thickness_pd*thickness, 0, inf)
+    elif thickness_pd_type == 'schulz':
+        T_val, T_prob = schulz_distribution(thickness, thickness_pd*thickness, 0, inf)
+    total_weight = total_volume = total_shell = 0
+    radius_eff = 0
+    F1, F2 = np.zeros_like(q), np.zeros_like(q)
+    for k, Rk in enumerate(R_val):
+        print("vesicle cycle", k, "of", len(R_val)," Rk ",Rk)
+        volume_in = 4./3.*pi*Rk**3
+        form_in = (sld-sld_solvent)*volume_in*sas_3j1x_x(q*Rk)
+        for i, Ti in enumerate(T_val):
+            Rout = Rk + Ti
+            volume_out = 4./3.*pi*Rout**3
+            form_out = (sld-sld_solvent)*volume_out*sas_3j1x_x(q*Rout)
+            form = form_out - form_in
+            volume = volume_out - volume_in
+            weight = R_prob[k]*T_prob[i]
+            total_weight += weight
+            total_shell += weight*volume
+            total_volume += weight*volume_out
+            F1 += form*weight
+            F2 += form**2*weight
+            radius_eff += weight*Rout
+    F1 /= total_weight
+    F2 /= total_weight
+    average_volume = total_volume/total_weight
+    if radius_effective is None:
+        radius_effective = radius_eff/total_weight
+        print("radius_effective  ",radius_effective)
+    if norm == 'sasfit':
+        IQD = F2
+    elif norm == 'sasview':
+        # Note: internally, sasview uses F2/total_volume because:
+        #   average_volume = total_volume/total_weight
+        #   F2/total_weight / average_volume
+        #     = F2/total_weight / total_volume/total_weight
+        #     = F2/total_volume
+        IQD = F2/average_volume*1e-4*volfraction
+        F1 *= 1e-2  # Yun is using sld in 1/A^2, not 1e-6/A^2
+        F2 *= 1e-4
+    elif norm == 'yun':
+        F1 *= 1e-6  # Yun is using sld in 1/A^2, not 1e-6/A^2
+        F2 *= 1e-12
+        IQD = F2/average_volume*1e8*volfraction
+    # RKH SORT THIS OUT WHEN HAVE NEW MODEL
+    vfff=volfraction*(30./20.)**3
+    print("radius_effective, fudged volfaction for S(Q) ",radius_effective,vfff)
+    SQ = hardsphere_simple(q, radius_effective, vfff)
+    beta = F1**2/F2
+    SQ_EFF = 1 + beta*(SQ - 1)
+    IQSD = IQD*SQ
+    IQBD = IQD*SQ_EFF
+    return Theory(Q=q, F1=F1, F2=F2, P=IQD, S=SQ, I=IQSD, Seff=SQ_EFF, Ibeta=IQBD)
+#
+#polydispersity for hollow_cylinder
+def hollow_cylinder_pe(q,radius=20, thickness=10, sld=4, sld_solvent=1, volfraction=0.15,
+        radius_pd=0.1, thickness_pd=0.2, length_pd=0.05, radius_pd_type="gaussian", thickness_pd_type="gaussian",
+        radius_effective=None, background=0, scale=1,
+                 norm='sasview'):
+    if norm not in ['sasview', 'sasfit', 'yun']:
+        raise TypeError("unknown norm "+norm)
+    if radius_pd_type == 'gaussian':
+        R_val, R_prob = gaussian_distribution(radius, radius_pd*radius, 0, inf)
+    elif radius_pd_type == 'schulz':
+        R_val, R_prob = schulz_distribution(radius, radius_pd*radius, 0, inf)
+    if thickness_pd_type == 'gaussian':
+        T_val, T_prob = gaussian_distribution(thickness, thickness_pd*thickness, 0, inf)
+    elif thickness_pd_type == 'schulz':
+        T_val, T_prob = schulz_distribution(thickness, thickness_pd*thickness, 0, inf)
+    if length_pd_type == 'gaussian':
+        T_val, T_prob = gaussian_distribution(length, length_pd*length, 0, inf)
+    elif length_pd_type == 'schulz':
+        L_val, L_prob = schulz_distribution(length, length_pd*length, 0, inf)
+    total_weight = total_volume = total_shell = 0
+    radius_eff = 0
+    F1, F2 = np.zeros_like(q), np.zeros_like(q)
+    for k, Rk in enumerate(R_val):
+        print("hollow_cylinder cycle", k, "of", len(R_val)," Rk ",Rk)
+        for i, Ti in enumerate(T_val):
+            for j, Lj in enumerate(L_val):
+                Rout = Rk + Ti
+                volume_out = pi*Rout**2*Lj
+                volume_in = pi*Rk**2*Lj
+                volume = volume_out - volume_in
+                theory = hollow_cylinder_theta(q, Rk, Ti, Lj, sld, sld_solvent,
+                    volfraction=0, radius_effective=None)
+                weight = R_prob[k]*T_prob[i]*L_prob[j]
+                total_weight += weight
+                total_shell += weight*volume
+                total_volume += weight*volume_out
+                F1 += theory.F1*weight
+                F2 += theory.F2*weight
+                radius_eff += weight*ER_hollow_cylinder(radius, thickness, length)
+    F1 /= total_weight
+    F2 /= total_weight
+    average_volume = total_volume/total_weight
+    if radius_effective is None:
+        radius_effective = radius_eff/total_weight
+        print("radius_effective  ",radius_effective)
+    if norm == 'sasfit':
+        IQD = F2
+    elif norm == 'sasview':
+        # Note: internally, sasview uses F2/total_volume because:
+        #   average_volume = total_volume/total_weight
+        #   F2/total_weight / average_volume
+        #     = F2/total_weight / total_volume/total_weight
+        #     = F2/total_volume
+        IQD = F2/average_volume*1e-4*volfraction
+        F1 *= 1e-2  # Yun is using sld in 1/A^2, not 1e-6/A^2
+        F2 *= 1e-4
+    elif norm == 'yun':
+        F1 *= 1e-6  # Yun is using sld in 1/A^2, not 1e-6/A^2
+        F2 *= 1e-12
+        IQD = F2/average_volume*1e8*volfraction
+# RKH SORT THIS OUT WHEN HAVE NEW MODEL
+    vfff=volfraction
+    print("radius_effective, volfaction for S(Q) ",radius_effective,vfff)
+    SQ = hardsphere_simple(q, radius_effective, vfff)
+    beta = F1**2/F2
+    SQ_EFF = 1 + beta*(SQ - 1)
+    IQSD = IQD*SQ
+    IQBD = IQD*SQ_EFF
+    return Theory(Q=q, F1=F1, F2=F2, P=IQD, S=SQ, I=IQSD, Seff=SQ_EFF, Ibeta=IQBD)
+
+# RKH copying Paul K's ellipsoid here, this computes everything for monodisperse case
+# and supplies partial results for polydisperse case.
+def hollow_cylinder_theta(q, radius, thickness, length, sld, sld_solvent,
+                    volfraction=0, radius_effective=None):
+    #creates values z and corresponding probabilities w from legendre-gauss quadrature
+    z, w = leggauss(76)
+    F1 = np.zeros_like(q)
+    F2 = np.zeros_like(q)
+    #use a u subsition(u=cos) and then u=(z+1)/2 to change integration from
+    #0->2pi with respect to alpha to -1->1 with respect to z, allowing us to use
+    #legendre-gauss quadrature
+    lower = 0.0
+    upper = 1.0
+    gamma_sq = (radius/(radius+thickness))**2
+    for k, qk in enumerate(q):
+        for i, zt in enumerate(z):
+            # quadrature loop
+            cos_theta = 0.5*( zt* (upper-lower) + lower + upper  )
+            sin_theta = sqrt(1.0 - cos_theta*cos_theta)
+            aaa=(radius+thickness)*qk*sin_theta
+            ## sas_J1 expects an array, so try direct use of J1
+            lam1 = J1(aaa)/aaa
+            aaa=radius*qk*sin_theta
+            lam2 = J1(aaa)/aaa
+            psi = (lam1 - gamma_sq*lam2)/(1.0 - gamma_sq)
+            #Note: lim_{thickness -> 0} psi = sas_J0(radius*qab)
+            #Note: lim_{radius -> 0} psi = sas_2J1x_x(thickness*qab)
+            aaa=0.5*length*qk*cos_theta
+            t2 = sin(aaa)/aaa
+#            t2 = sas_sinx_x(0.5*length*qk*cos_theta)
+            form = psi*t2
+            F1[i] += w[i] * form
+            F2[i] += w[i] * form * form
+    volume = pi*((radius + thickness)**2 - radius**2)*length
+    s = (sld - sld_solvent) * volume
+    F2 = F2*s*(upper - lower)/2.0
+    F1 = F1*s*s*(upper - lower)/2.0
+    if radius_effective is None:
+        radius_effective = ER_hollow_cylinder(radius, thickness, length)
+    SQ = hardsphere_simple(q, radius_effective, volfraction)
+    SQ_EFF = 1 + F1**2/F2*(SQ - 1)
+    IQM = 1e-4*F2/volume
+    IQSM = volfraction*IQM*SQ
+    IQBM = volfraction*IQM*SQ_EFF
+    return Theory(Q=q, F1=F1, F2=F2, P=IQM, S=SQ, I=IQSM, Seff=SQ_EFF, Ibeta=IQBM)
+#
+def ER_hollow_cylinder(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
+
+###############################################################################
+###############################################################################
+###############################################################################
+##################                                           ##################
+##################                   TESTS                   ##################
+##################                                           ##################
+###############################################################################
+###############################################################################
+###############################################################################
+
+def popn(d, keys):
+    """
+    Splits a dict into two, with any key of *d* which is in *keys* removed
+    from *d* and added to *b*. Returns *b*.
+    """
+    b = {}
+    for k in keys:
+        try:
+            b[k] = d.pop(k)
+        except KeyError:
+            pass
+    return b
+
+def sasmodels_theory(q, Pname, **pars):
+    """
+    Call sasmodels to compute the model with and without a hard sphere
+    structure factor.
+    """
+    #mono_pars = {k: (0 if k.endswith('_pd') else v) for k, v in pars.items()}
+    Ppars = pars.copy()
+    Spars = popn(Ppars, ['radius_effective', 'volfraction'])
+    Ipars = pars.copy()
+
+    # Autofill npts and nsigmas for the given pd type
+    for k, v in pars.items():
+        if k.endswith("_pd_type"):
+            if v == "gaussian":
+                n, nsigmas = N_GAUSS, NSIGMA_GAUSS
+            elif v == "schulz":
+                n, nsigmas = N_SCHULZ, NSIGMA_SCHULZ
+            Ppars.setdefault(k.replace("_pd_type", "_pd_n"), n)
+            Ppars.setdefault(k.replace("_pd_type", "_pd_nsigma"), nsigmas)
+            Ipars.setdefault(k.replace("_pd_type", "_pd_n"), n)
+            Ipars.setdefault(k.replace("_pd_type", "_pd_nsigma"), nsigmas)
+
+    #Ppars['scale'] = Spars.get('volfraction', 1)
+    P = build_model(Pname, q)
+    S = build_model("hardsphere", q)
+    I = build_model(Pname+"@hardsphere", q)
+    Pq = P(**Ppars)*pars.get('volfraction', 1)
+    Sq = S(**Spars)
+    Iq = I(**Ipars)
+    #Iq = Pq*Sq*pars.get('volfraction', 1)
+    #Sq = Iq/Pq
+    #Iq = None#= Sq = None
+    r = dict(I._kernel.results())
+    #print(r)
+    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'):
+    """
+    Plot fields in common between target and actual, along with relative error.
+    """
+    available = [s for s in fields.split()
+                 if getattr(target, s) is not None and getattr(actual, s) is not None]
+    rows = len(available)
+    for row, field in enumerate(available):
+        Q = target.Q
+        I1, I2 = getattr(target, field), getattr(actual, field)
+        plt.subplot(rows, 2, 2*row+1)
+        plt.loglog(Q, abs(I1), label="target "+field)
+        plt.loglog(Q, abs(I2), label="value "+field)
+        #plt.legend(loc="upper left", bbox_to_anchor=(1,1))
+        plt.legend(loc='lower left')
+        plt.subplot(rows, 2, 2*row+2)
+        plt.semilogx(Q, I2/I1 - 1, label="relative error")
+        #plt.semilogx(Q, I1/I2 - 1, label="relative error")
+    plt.tight_layout()
+    plt.suptitle(title)
+    plt.show()
+
+def data_file(name):
+    return os.path.join(BETA_DIR, 'data_files', name)
+
+def load_sasfit(path):
+    data = np.loadtxt(path, dtype=str, delimiter=';').T
+    data = np.vstack((map(float, v) for v in data[0:2]))
+    return data
+
+COMPARISON = {}  # Type: Dict[(str,str,str)] -> Callable[(), None]
+
+def compare_sasview_sphere(pd_type='schulz'):
+    q = np.logspace(-5, 0, 250)
+    model = 'sphere'
+    pars = dict(
+        radius=20, sld=4, sld_solvent=1,
+        background=0,
+        radius_pd=.1, radius_pd_type=pd_type,
+        volfraction=0.15,
+        radius_effective=12.59921049894873,  # equivalent average sphere radius
+        )
+    target = sasmodels_theory(q, model, effective_radius_mode=0, structure_factor_mode=1, **pars)
+    actual = sphere_r(q, norm='sasview', **pars)
+    title = " ".join(("sasmodels", model, pd_type))
+    compare(title, target, actual)
+COMPARISON[('sasview', 'sphere', 'gaussian')] = lambda: compare_sasview_sphere(pd_type='gaussian')
+COMPARISON[('sasview', 'sphere', 'schulz')] = lambda: compare_sasview_sphere(pd_type='schulz')
+
+
+def compare_sasview_vesicle(pd_type='gaussian'):
+    q = np.logspace(-5, 0, 250)
+    model = 'vesicle'
+    print("F2 seems OK, but big issues with S(Q), wo work in progress")
+# NOTE parameers for vesicle model are soon to be changed to remove "volfraction"
+    pars = dict(
+        radius=20, thickness=10, sld=4, sld_solvent=1, volfraction=0.03,
+        background=0,
+        radius_pd=0.1, thickness_pd=0.0, radius_pd_type=pd_type, thickness_pd_type=pd_type
+        #,radius_effective=12.59921049894873,  # equivalent average sphere radius
+        )
+    target = sasmodels_theory(q, model, effective_radius_mode=0, structure_factor_mode=1, **pars)
+    actual = vesicle_pe(q, norm='sasview', **pars)
+    title = " ".join(("sasmodels", model, pd_type))
+    compare(title, target, actual)
+COMPARISON[('sasview', 'vesicle', 'gaussian')] = lambda: compare_sasview_vesicle(pd_type='gaussian')
+COMPARISON[('sasview', 'vesicle', 'schulz')] = lambda: compare_sasview_vesicle(pd_type='schulz')
+
+def compare_sasview_hollow_cylinder(pd_type='gaussian'):
+    q = np.logspace(-5, 0, 250)
+    model = 'hollow_cylinder'
+    print(model)
+    print("just about works for monodisperse, polydisperse needs work")
+# NOTE parameters for hollow_cylinder model are soon to be changed to remove "volfraction"
+# setting all three pd to zero is OK
+    pars = dict(
+        radius=20, thickness=10, length=80, sld=4, sld_solvent=1,
+        background=0,
+        radius_pd=0.1, thickness_pd=0.0, length_pd=0.0, radius_pd_type=pd_type, thickness_pd_type=pd_type, length_pd_type=pd_type
+        #,radius_effective=12.59921049894873,  # equivalent average sphere radius
+        )
+    target = sasmodels_theory(q, model, effective_radius_mode=0, structure_factor_mode=1, **pars)
+    actual = hollow_cylinder_pe(q, norm='sasview', **pars)
+# RKH monodisp was OK,    actual = hollow_cylinder_theta(q,radius=20, thickness=10, length=80, sld=4, sld_solvent=1 )
+    title = " ".join(("sasmodels", model, pd_type))
+    compare(title, target, actual)
+COMPARISON[('sasview', 'hollow_cylinder', 'gaussian')] = lambda: compare_sasview_hollow_cylinder(pd_type='gaussian')
+COMPARISON[('sasview', 'hollow_cylinder', 'schulz')] = lambda: compare_sasview_hollow_cylinder(pd_type='schulz')
+
+
+def compare_sasview_ellipsoid(pd_type='gaussian'):
+    q = np.logspace(-5, 0, 50)
+    model = 'ellipsoid'
+    pars = dict(
+        radius_polar=20, radius_equatorial=400, sld=4, sld_solvent=1,
+        background=0,
+        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, effective_radius_mode=0, structure_factor_mode=1, **pars)
+    actual = ellipsoid_pe(q, norm='sasview', **pars)
+# RKH test       actual = ellipsoid_theta(q, radius_polar=20, radius_equatorial=400, sld=4, sld_solvent=1, volfraction=0.15, radius_effective=270.)
+    title = " ".join(("sasmodels", model, pd_type))
+    compare(title, target, actual)
+COMPARISON[('sasview', 'ellipsoid', 'gaussian')] = lambda: compare_sasview_ellipsoid(pd_type='gaussian')
+COMPARISON[('sasview', 'ellipsoid', 'schulz')] = lambda: compare_sasview_ellipsoid(pd_type='schulz')
+
+def compare_yun_ellipsoid_mono():
+    pars = {
+        'radius_polar': 20, 'radius_polar_pd': 0, 'radius_polar_pd_type': 'gaussian',
+        'radius_equatorial': 10, 'radius_equatorial_pd': 0, 'radius_equatorial_pd_type': 'gaussian',
+        'sld': 2, 'sld_solvent': 1,
+        'volfraction': 0.15,
+        # Yun uses radius for same volume sphere for effective radius
+        # whereas sasview uses the average curvature.
+        'radius_effective': 12.59921049894873,
+    }
+
+    data = np.loadtxt(data_file('yun_ellipsoid.dat'),skiprows=2).T
+    Q = data[0]
+    F1 = data[1]
+    P = data[3]*pars['volfraction']
+    S = data[5]
+    Seff = data[6]
+    target = Theory(Q=Q, F1=F1, P=P, S=S, Seff=Seff)
+    actual = ellipsoid_pe(Q, norm='yun', **pars)
+    title = " ".join(("yun", "ellipsoid", "no pd"))
+    #compare(title, target, actual, fields="P S I Seff Ibeta")
+    compare(title, target, actual)
+COMPARISON[('yun', 'ellipsoid', 'gaussian')] = compare_yun_ellipsoid_mono
+COMPARISON[('yun', 'ellipsoid', 'schulz')] = compare_yun_ellipsoid_mono
+
+def compare_yun_sphere_gauss():
+    # Note: yun uses gauss limits from R0/10 to R0 + 5 sigma steps sigma/100
+    # With pd = 0.1, that's 14 sigma and 1400 points.
+    pars = {
+        'radius': 20, 'radius_pd': 0.1, 'radius_pd_type': 'gaussian',
+        'sld': 6, 'sld_solvent': 0,
+        'volfraction': 0.1,
+    }
+
+    data = np.loadtxt(data_file('testPolydisperseGaussianSphere.dat'), skiprows=2).T
+    Q = data[0]
+    F1 = data[1]
+    F2 = data[2]
+    P = data[3]
+    S = data[5]
+    Seff = data[6]
+    target = Theory(Q=Q, F1=F1, P=P, S=S, Seff=Seff)
+    actual = sphere_r(Q, norm='yun', **pars)
+    title = " ".join(("yun", "sphere", "10% dispersion 10% Vf"))
+    compare(title, target, actual)
+    data = np.loadtxt(data_file('testPolydisperseGaussianSphere2.dat'), skiprows=2).T
+    pars.update(radius_pd=0.15)
+    Q = data[0]
+    F1 = data[1]
+    F2 = data[2]
+    P = data[3]
+    S = data[5]
+    Seff = data[6]
+    target = Theory(Q=Q, F1=F1, P=P, S=S, Seff=Seff)
+    actual = sphere_r(Q, norm='yun', **pars)
+    title = " ".join(("yun", "sphere", "15% dispersion 10% Vf"))
+    compare(title, target, actual)
+COMPARISON[('yun', 'sphere', 'gaussian')] = compare_yun_sphere_gauss
+
+
+def compare_sasfit_sphere_gauss():
+    #N=1,s=2,X0=20,distr radius R=20,eta_core=4,eta_solv=1,.3
+    pars = {
+        'radius': 20, 'radius_pd': 0.1, 'radius_pd_type': 'gaussian',
+        'sld': 4, 'sld_solvent': 1,
+        'volfraction': 0.3,
+    }
+
+    Q, IQ = load_sasfit(data_file('sasfit_sphere_IQD.txt'))
+    Q, IQSD = load_sasfit(data_file('sasfit_sphere_IQSD.txt'))
+    Q, IQBD = load_sasfit(data_file('sasfit_sphere_IQBD.txt'))
+    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_sphere_sq.txt'))
+    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_sphere_sqeff.txt'))
+    target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=SQ, I=IQSD, Seff=SQ_EFF, Ibeta=IQBD)
+    actual = sphere_r(Q, norm="sasfit", **pars)
+    title = " ".join(("sasfit", "sphere", "pd=10% gaussian"))
+    compare(title, target, actual)
+    #compare(title, target, actual, fields="P")
+COMPARISON[('sasfit', 'sphere', 'gaussian')] = compare_sasfit_sphere_gauss
+
+def compare_sasfit_sphere_schulz():
+    #radius=20,sld=4,sld_solvent=1,volfraction=0.3,radius_pd=0.1
+    #We have scaled the output from sasfit by 1e-4*volume*volfraction
+    #0.10050378152592121
+    pars = {
+        'radius': 20, 'radius_pd': 0.1, 'radius_pd_type': 'schulz',
+        'sld': 4, 'sld_solvent': 1,
+        'volfraction': 0.3,
+    }
+
+    Q, IQ = load_sasfit(data_file('richard_test.txt'))
+    Q, IQSD = load_sasfit(data_file('richard_test2.txt'))
+    Q, IQBD = load_sasfit(data_file('richard_test3.txt'))
+    target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=None, I=IQSD, Seff=None, Ibeta=IQBD)
+    actual = sphere_r(Q, norm="sasfit", **pars)
+    title = " ".join(("sasfit", "sphere", "pd=10% schulz"))
+    compare(title, target, actual)
+COMPARISON[('sasfit', 'sphere', 'schulz')] = compare_sasfit_sphere_schulz
+
+def compare_sasfit_ellipsoid_schulz():
+    #polarradius=20, equatorialradius=10, sld=4,sld_solvent=1,volfraction=0.3,radius_polar_pd=0.1
+    #Effective radius =13.1353356684
+    #We have scaled the output from sasfit by 1e-4*volume*volfraction
+    #0.10050378152592121
+    pars = {
+        'radius_polar': 20, 'radius_polar_pd': 0.1, 'radius_polar_pd_type': 'schulz',
+        'radius_equatorial': 10, 'radius_equatorial_pd': 0., 'radius_equatorial_pd_type': 'schulz',
+        'sld': 4, 'sld_solvent': 1,
+        'volfraction': 0.3, 'radius_effective': 13.1353356684,
+    }
+
+    Q, IQ = load_sasfit(data_file('richard_test4.txt'))
+    Q, IQSD = load_sasfit(data_file('richard_test5.txt'))
+    Q, IQBD = load_sasfit(data_file('richard_test6.txt'))
+    target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=None, I=IQSD, Seff=None, Ibeta=IQBD)
+    actual = ellipsoid_pe(Q, norm="sasfit", **pars)
+    title = " ".join(("sasfit", "ellipsoid", "pd=10% schulz"))
+    compare(title, target, actual)
+COMPARISON[('sasfit', 'ellipsoid', 'schulz')] = compare_sasfit_ellipsoid_schulz
+
+
+def compare_sasfit_ellipsoid_gaussian():
+    pars = {
+        'radius_polar': 20, 'radius_polar_pd': 0, 'radius_polar_pd_type': 'gaussian',
+        'radius_equatorial': 10, 'radius_equatorial_pd': 0, 'radius_equatorial_pd_type': 'gaussian',
+        'sld': 4, 'sld_solvent': 1,
+        'volfraction': 0, 'radius_effective': None,
+    }
+
+    #Rp=20,Re=10,eta_core=4,eta_solv=1
+    Q, PQ0 = load_sasfit(data_file('sasfit_ellipsoid_IQM.txt'))
+    pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0, radius_effective=None)
+    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
+    target = Theory(Q=Q, P=PQ0)
+    compare("sasfit ellipsoid no poly", target, actual); plt.show()
+
+    #N=1,s=2,X0=20,distr 10% polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
+    Q, PQ_Rp10 = load_sasfit(data_file('sasfit_ellipsoid_IQD.txt'))
+    pars.update(volfraction=0, radius_polar_pd=0.1, radius_equatorial_pd=0.0, radius_effective=None)
+    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
+    target = Theory(Q=Q, P=PQ_Rp10)
+    compare("sasfit ellipsoid P(Q) 10% Rp", target, actual); plt.show()
+    #N=1,s=1,X0=10,distr 10% equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
+    Q, PQ_Re10 = load_sasfit(data_file('sasfit_ellipsoid_IQD2.txt'))
+    pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0.1, radius_effective=None)
+    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
+    target = Theory(Q=Q, P=PQ_Re10)
+    compare("sasfit ellipsoid P(Q) 10% Re", target, actual); plt.show()
+    #N=1,s=6,X0=20,distr 30% polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
+    Q, PQ_Rp30 = load_sasfit(data_file('sasfit_ellipsoid_IQD3.txt'))
+    pars.update(volfraction=0, radius_polar_pd=0.3, radius_equatorial_pd=0.0, radius_effective=None)
+    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
+    target = Theory(Q=Q, P=PQ_Rp30)
+    compare("sasfit ellipsoid P(Q) 30% Rp", target, actual); plt.show()
+    #N=1,s=3,X0=10,distr 30% equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
+    Q, PQ_Re30 = load_sasfit(data_file('sasfit_ellipsoid_IQD4.txt'))
+    pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0.3, radius_effective=None)
+    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
+    target = Theory(Q=Q, P=PQ_Re30)
+    compare("sasfit ellipsoid P(Q) 30% Re", target, actual); plt.show()
+    #N=1,s=12,X0=20,distr 60% polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
+    Q, PQ_Rp60 = load_sasfit(data_file('sasfit_ellipsoid_IQD5.txt'))
+    pars.update(volfraction=0, radius_polar_pd=0.6, radius_equatorial_pd=0.0, radius_effective=None)
+    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
+    target = Theory(Q=Q, P=PQ_Rp60)
+    compare("sasfit ellipsoid P(Q) 60% Rp", target, actual); plt.show()
+    #N=1,s=6,X0=10,distr 60% equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
+    Q, PQ_Re60 = load_sasfit(data_file('sasfit_ellipsoid_IQD6.txt'))
+    pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0.6, radius_effective=None)
+    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
+    target = Theory(Q=Q, P=PQ_Re60)
+    compare("sasfit ellipsoid P(Q) 60% Re", target, actual); plt.show()
+
+    #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.15
+    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq.txt'))
+    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff.txt'))
+    pars.update(volfraction=0.15, radius_polar_pd=0.1, radius_equatorial_pd=0, radius_effective=13.1354236254)
+    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
+    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
+    compare("sasfit ellipsoid P(Q) 10% Rp 15% Vf", target, actual); plt.show()
+    #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.15
+    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq2.txt'))
+    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff2.txt'))
+    pars.update(volfraction=0.15, radius_polar_pd=0.3, radius_equatorial_pd=0, radius_effective=13.0901197149)
+    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
+    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
+    compare("sasfit ellipsoid P(Q) 30% Rp 15% Vf", target, actual); plt.show()
+    #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.15
+    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq3.txt'))
+    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff3.txt'))
+    pars.update(volfraction=0.15, radius_polar_pd=0.6, radius_equatorial_pd=0, radius_effective=13.336060917)
+    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
+    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
+    compare("sasfit ellipsoid P(Q) 60% Rp 15% Vf", target, actual); plt.show()
+
+    #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.3
+    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq4.txt'))
+    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff4.txt'))
+    pars.update(volfraction=0.3, radius_polar_pd=0.1, radius_equatorial_pd=0, radius_effective=13.1354236254)
+    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
+    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
+    compare("sasfit ellipsoid P(Q) 10% Rp 30% Vf", target, actual); plt.show()
+    #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.3
+    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq5.txt'))
+    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff5.txt'))
+    pars.update(volfraction=0.3, radius_polar_pd=0.3, radius_equatorial_pd=0, radius_effective=13.0901197149)
+    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
+    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
+    compare("sasfit ellipsoid P(Q) 30% Rp 30% Vf", target, actual); plt.show()
+    #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.3
+    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq6.txt'))
+    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff6.txt'))
+    pars.update(volfraction=0.3, radius_polar_pd=0.6, radius_equatorial_pd=0, radius_effective=13.336060917)
+    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
+    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
+    compare("sasfit ellipsoid P(Q) 60% Rp 30% Vf", target, actual); plt.show()
+
+    #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.6
+    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq7.txt'))
+    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff7.txt'))
+    pars.update(volfraction=0.6, radius_polar_pd=0.1, radius_equatorial_pd=0, radius_effective=13.1354236254)
+    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
+    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
+    compare("sasfit ellipsoid P(Q) 10% Rp 60% Vf", target, actual); plt.show()
+    #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.6
+    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq8.txt'))
+    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff8.txt'))
+    pars.update(volfraction=0.6, radius_polar_pd=0.3, radius_equatorial_pd=0, radius_effective=13.0901197149)
+    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
+    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
+    compare("sasfit ellipsoid P(Q) 30% Rp 60% Vf", target, actual); plt.show()
+    #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.6
+    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq9.txt'))
+    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff9.txt'))
+    pars.update(volfraction=0.6, radius_polar_pd=0.6, radius_equatorial_pd=0, radius_effective=13.336060917)
+    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
+    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
+    compare("sasfit ellipsoid P(Q) 60% Rp 60% Vf", target, actual); plt.show()
+COMPARISON[('sasfit', 'ellipsoid', 'gaussian')] = compare_sasfit_ellipsoid_gaussian
+
+def main():
+    key = tuple(sys.argv[1:])
+    if key not in COMPARISON:
+        print("Usage: python sasfit_compare.py [sasview|sasfit|yun] [sphere|ellipsoid|vesicle|hollow_cylinder] [gaussian|schulz]")
+        print("But not for [yun sphere schulz], and vesicle & hollow_cylinder only with sasview.")
+        print("Note this compares chosen source against internal python code here, with adjustment for known scale etc issues.")
+        return
+    comparison = COMPARISON[key]
+    comparison()
+
+if __name__ == "__main__":
+    main()
Index: explore/precision.py
===================================================================
--- explore/precision.py	(revision 2a7e20e9c4ba72da5c6088e743a64f3345ee554d)
+++ explore/precision.py	(revision ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -357,10 +357,12 @@
 )
 add_function(
-    name="debye",
+    name="gauss_coil",
     mp_function=lambda x: 2*(mp.exp(-x**2) + x**2 - 1)/x**4,
     np_function=lambda x: 2*(np.expm1(-x**2) + x**2)/x**4,
     ocl_function=make_ocl("""
     const double qsq = q*q;
-    if (qsq < 1.0) { // Pade approximation
+    // For double: use O(5) Pade with 0.5 cutoff (10 mad + 1 divide)
+    // For single: use O(7) Taylor with 0.8 cutoff (7 mad)
+    if (qsq < 0.0) {
         const double x = qsq;
         if (0) { // 0.36 single
@@ -372,5 +374,5 @@
             const double B1=3./8., B2=3./56., B3=1./336.;
             return (((A3*x + A2)*x + A1)*x + 1.)/(((B3*x + B2)*x + B1)*x + 1.);
-        } else if (1) { // 1.0 for single, 0.25 for double
+        } else if (0) { // 1.0 for single, 0.25 for double
             // PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 8}]
             const double A1=1./15., A2=1./60, A3=0., A4=1./75600.;
@@ -385,5 +387,5 @@
                   /(((((B5*x + B4)*x + B3)*x + B2)*x + B1)*x + 1.);
         }
-    } else if (qsq < 1.) { // Taylor series; 0.9 for single, 0.25 for double
+    } else if (qsq < 0.8) {
         const double x = qsq;
         const double C0 = +1.;
Index: sasmodels/core.py
===================================================================
--- sasmodels/core.py	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/core.py	(revision ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -363,6 +363,8 @@
     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 radius_effective_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 ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -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 NotImplementedError()
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 ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -62,4 +62,33 @@
 }
 
+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)
+{
+    switch (mode) {
+    case 1: // equivalent sphere
+        return radius_from_volume(radius_bell, radius , length);
+    case 2: // radius
+        return radius;
+    case 3: // half length
+        return 0.5*length;
+    case 4: // half total length
+        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 ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -117,4 +117,7 @@
 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 ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -84,4 +84,33 @@
 }
 
+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)
+{
+    switch (mode) {
+    case 1: // equivalent sphere
+        return radius_from_volume(radius, radius_cap, length);
+    case 2: // radius
+        return radius;
+    case 3: // half length
+        return 0.5*length;
+    case 4: // half total length
+        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 ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -137,4 +137,7 @@
 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 ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -9,5 +9,5 @@
 
 static double
-form_volume(double core_radius, double fp_n, double thickness[])
+outer_radius(double core_radius, double fp_n, double thickness[])
 {
   double r = core_radius;
@@ -16,5 +16,22 @@
     r += thickness[i];
   }
-  return M_4PI_3 * cube(r);
+  return r;
+}
+
+static double
+form_volume(double core_radius, double fp_n, double thickness[])
+{
+  return M_4PI_3 * cube(outer_radius(core_radius, fp_n, thickness));
+}
+
+static double
+effective_radius(int mode, double core_radius, double fp_n, double thickness[])
+{
+  switch (mode) {
+  case 1: // outer radius
+    return outer_radius(core_radius, fp_n, thickness);
+  case 2: // core radius
+    return core_radius;
+  }
 }
 
Index: sasmodels/models/core_multi_shell.py
===================================================================
--- sasmodels/models/core_multi_shell.py	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/core_multi_shell.py	(revision ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -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,9 +145,4 @@
     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
-
 demo = dict(sld_core=6.4,
             radius=60,
Index: sasmodels/models/core_shell_bicelle.c
===================================================================
--- sasmodels/models/core_shell_bicelle.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/core_shell_bicelle.c	(revision ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -34,4 +34,34 @@
 
     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)
+{
+    switch (mode) {
+    case 1: // equivalent sphere
+        return radius_from_volume(radius, thick_rim, thick_face, length);
+    case 2: // outer rim radius
+        return radius + thick_rim;
+    case 3: // half outer thickness
+        return 0.5*length + thick_face;
+    case 4: // half diagonal
+        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 ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -155,4 +155,8 @@
           "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 ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -8,4 +8,39 @@
 {
     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)
+{
+    switch (mode) {
+    case 1: // equivalent sphere
+        return radius_from_volume(r_minor, x_core, thick_rim, thick_face, length);
+    case 2: // outer rim average radius
+        return 0.5*r_minor*(1.0 + x_core) + thick_rim;
+    case 3: // outer rim min radius
+        return (x_core < 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim);
+    case 4: // outer max radius
+        return (x_core > 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim);
+    case 5: // half outer thickness
+        return 0.5*length + thick_face;
+    case 6: // half diagonal
+        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 ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -147,4 +147,8 @@
           "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 d299327021d8f52877a96e17d973ae9f7ff21c73)
@@ -9,4 +9,40 @@
     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)
+{
+    switch (mode) {
+    case 1: // equivalent sphere
+        return radius_from_volume(r_minor, x_core, thick_rim, thick_face, length);
+    case 2: // outer rim average radius
+        return 0.5*r_minor*(1.0 + x_core) + thick_rim;
+    case 3: // outer rim min radius
+        return (x_core < 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim);
+    case 4: // outer max radius
+        return (x_core > 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim);
+    case 5: // half outer thickness
+        return 0.5*length + thick_face;
+    case 6: // half diagonal (this ignores the missing "corners", so may give unexpected answer if thick_face
+            //  or thick_rim is extremely large)
+        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 ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -160,4 +160,8 @@
           "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 ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -11,4 +11,38 @@
 {
     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 + 2.0*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)
+{
+    switch (mode) {
+    case 1: // equivalent sphere
+        return radius_from_volume(radius, thickness, length);
+    case 2: // outer radius
+        return radius + thickness;
+    case 3: // half outer length
+        return 0.5*length + thickness;
+    case 4: // half min outer length
+        return (radius < 0.5*length ? radius + thickness : 0.5*length + thickness);
+    case 5: // half max outer length
+        return (radius > 0.5*length ? radius + thickness : 0.5*length + thickness);
+    case 6: // half outer diagonal
+        return radius_from_diagonal(radius,thickness,length);
+    }
 }
 
Index: sasmodels/models/core_shell_cylinder.py
===================================================================
--- sasmodels/models/core_shell_cylinder.py	(revision e31b19a46729ff0f8f1a5eb9847f69524a6ed1d3)
+++ sasmodels/models/core_shell_cylinder.py	(revision ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -132,13 +132,9 @@
 source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "core_shell_cylinder.c"]
 have_Fq = True
-
-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.)
+effective_radius_type = [
+    "equivalent sphere", "outer radius", "half outer length",
+    "half min outer dimension", "half max outer dimension",
+    "half outer diagonal",
+    ]
 
 def random():
Index: sasmodels/models/core_shell_ellipsoid.c
===================================================================
--- sasmodels/models/core_shell_ellipsoid.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/core_shell_ellipsoid.c	(revision ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -36,4 +36,51 @@
     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;
+
+    switch (mode) {
+    case 1: // equivalent sphere
+        return radius_from_volume(radius_equat_core, x_core, thick_shell, x_polar_shell);
+    case 2: // average outer curvature
+        return radius_from_curvature(radius_equat_core, x_core, thick_shell, x_polar_shell);
+    case 3: // min outer radius
+        return (radius_polar_tot < radius_equat_tot ? radius_polar_tot : radius_equat_tot);
+    case 4: // max outer radius
+        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 ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -146,13 +146,8 @@
 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 random():
Index: sasmodels/models/core_shell_parallelepiped.c
===================================================================
--- sasmodels/models/core_shell_parallelepiped.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/core_shell_parallelepiped.c	(revision ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -25,4 +25,41 @@
         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 = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c);
+    return cbrt(volume/M_4PI_3);
+}
+
+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)
+{
+    switch (mode) {
+    case 1: // equivalent sphere
+        return radius_from_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c);
+    case 2: // half outer length a
+        return 0.5 * length_a + thick_rim_a;
+    case 3: // half outer length b
+        return 0.5 * length_b + thick_rim_b;
+    case 4: // half outer length c
+        return 0.5 * length_c + thick_rim_c;
+    case 5: // equivalent circular cross-section
+        return radius_from_crosssection(length_a, length_b, thick_rim_a, thick_rim_b);
+    case 6: // half outer ab diagonal
+        return 0.5*sqrt(square(length_a+ 2.0*thick_rim_a) + square(length_b+ 2.0*thick_rim_b));
+    case 7: // half outer diagonal
+        return 0.5*sqrt(square(length_a+ 2.0*thick_rim_a) + square(length_b+ 2.0*thick_rim_b) + square(length_c+ 2.0*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 ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -98,5 +98,5 @@
 is the scattering length of the solvent.
 
-.. note:: 
+.. note::
 
    the code actually implements two substitutions: $d(cos\alpha)$ is
@@ -227,18 +227,10 @@
 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)
-
-# VR defaults to 1.0
+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 random():
Index: sasmodels/models/core_shell_sphere.c
===================================================================
--- sasmodels/models/core_shell_sphere.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/core_shell_sphere.c	(revision ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -3,4 +3,15 @@
 {
     return M_4PI_3 * cube(radius + thickness);
+}
+
+static double
+effective_radius(int mode, double radius, double thickness)
+{
+    switch (mode) {
+    case 1: // outer radius
+        return radius + thickness;
+    case 2: // core radius
+        return radius;
+    }
 }
 
Index: sasmodels/models/core_shell_sphere.py
===================================================================
--- sasmodels/models/core_shell_sphere.py	(revision da1c8d14693f7a806280e49df578d46f9aa451d8)
+++ sasmodels/models/core_shell_sphere.py	(revision ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -78,15 +78,8 @@
 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 random():
@@ -103,5 +96,5 @@
 
 tests = [
-    [{'radius': 20.0, 'thickness': 10.0}, 'ER', 30.0],
+#    [{'radius': 20.0, 'thickness': 10.0}, 'ER', 30.0],
 
     # The SasView test result was 0.00169, with a background of 0.001
Index: sasmodels/models/cylinder.c
===================================================================
--- sasmodels/models/cylinder.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/cylinder.c	(revision ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -11,4 +11,35 @@
 {
     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)
+{
+    switch (mode) {
+    case 1:
+        return radius_from_volume(radius, length);
+    case 2:
+        return radius;
+    case 3:
+        return 0.5*length;
+    case 4:
+        return (radius < 0.5*length ? radius : 0.5*length);
+    case 5:
+        return (radius > 0.5*length ? radius : 0.5*length);
+    case 6:
+        return radius_from_diagonal(radius,length);
+    }
 }
 
Index: sasmodels/models/cylinder.py
===================================================================
--- sasmodels/models/cylinder.py	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/cylinder.py	(revision ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -139,11 +139,8 @@
 source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "cylinder.c"]
 have_Fq = True
-
-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.)
+effective_radius_type = [
+    "equivalent sphere", "radius",
+    "half length", "half min dimension", "half max dimension", "half diagonal",
+    ]
 
 def random():
Index: sasmodels/models/ellipsoid.c
===================================================================
--- sasmodels/models/ellipsoid.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/ellipsoid.c	(revision ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -4,4 +4,46 @@
     return M_4PI_3*radius_polar*radius_equatorial*radius_equatorial;
 }
+
+static double
+radius_from_volume(double radius_polar, double radius_equatorial)
+{
+    return cbrt(radius_polar*radius_equatorial*radius_equatorial);
+}
+
+static double
+radius_from_curvature(double radius_polar, double radius_equatorial)
+{
+    // Trivial cases
+    if (radius_polar == radius_equatorial) return radius_polar;
+    if (radius_polar * radius_equatorial == 0.)  return 0.;
+
+    // see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449
+    const double ratio = (radius_polar < radius_equatorial
+                          ? radius_polar / radius_equatorial
+                          : radius_equatorial / radius_polar);
+    const double e1 = sqrt(1.0 - ratio*ratio);
+    const double b1 = 1.0 + asin(e1) / (e1 * ratio);
+    const double bL = (1.0 + e1) / (1.0 - e1);
+    const double b2 = 1.0 + 0.5 * ratio * ratio / e1 * log(bL);
+    const double delta = 0.75 * b1 * b2;
+    const double ddd = 2.0 * (delta + 1.0) * radius_polar * radius_equatorial * radius_equatorial;
+    return 0.5 * cbrt(ddd);
+}
+
+static double
+effective_radius(int mode, double radius_polar, double radius_equatorial)
+{
+    switch (mode) {
+    case 1: // equivalent sphere
+        return radius_from_volume(radius_polar, radius_equatorial);
+    case 2: // average curvature
+        return radius_from_curvature(radius_polar, radius_equatorial);
+    case 3: // min radius
+        return (radius_polar < radius_equatorial ? radius_polar : radius_equatorial);
+    case 4: // max radius
+        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 ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -169,26 +169,7 @@
 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 ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -3,4 +3,57 @@
 {
     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 hlength)
+{
+    const double rad_min = (r_ratio > 1.0 ? radius_minor : r_ratio*radius_minor);
+    return (rad_min < hlength ? rad_min : hlength);
+}
+
+static double
+radius_from_max_dimension(double radius_minor, double r_ratio, double hlength)
+{
+    const double rad_max = (r_ratio < 1.0 ? radius_minor : r_ratio*radius_minor);
+    return (rad_max > hlength ? rad_max : hlength);
+}
+
+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)
+{
+    switch (mode) {
+    case 1: // equivalent sphere
+        return radius_from_volume(radius_minor, r_ratio, length);
+    case 2: // average radius
+        return 0.5*radius_minor*(1.0 + r_ratio);
+    case 3: // min radius
+        return (r_ratio > 1.0 ? radius_minor : r_ratio*radius_minor);
+    case 4: // max radius
+        return (r_ratio < 1.0 ? radius_minor : r_ratio*radius_minor);
+    case 5: // equivalent circular cross-section
+        return sqrt(radius_minor*radius_minor*r_ratio);
+    case 6: // half length
+        return 0.5*length;
+    case 7: // half min dimension
+        return radius_from_min_dimension(radius_minor,r_ratio,0.5*length);
+    case 8: // half max dimension
+        return radius_from_max_dimension(radius_minor,r_ratio,0.5*length);
+    case 9: // half diagonal
+        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 ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -123,20 +123,13 @@
 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,
             sld=4.0, sld_solvent=1.0, theta=10.0, phi=20, psi=30,
             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 random():
@@ -162,6 +155,6 @@
 
 tests = [
-    [{'radius_minor': 20.0, 'axis_ratio': 1.5, 'length':400.0}, 'ER', 79.89245454155024],
-    [{'radius_minor': 20.0, 'axis_ratio': 1.2, 'length':300.0}, 'VR', 1],
+#    [{'radius_minor': 20.0, 'axis_ratio': 1.5, 'length':400.0}, 'ER', 79.89245454155024],
+#    [{'radius_minor': 20.0, 'axis_ratio': 1.2, 'length':300.0}, 'VR', 1],
 
     # The SasView test result was 0.00169, with a background of 0.001
Index: sasmodels/models/fractal_core_shell.py
===================================================================
--- sasmodels/models/fractal_core_shell.py	(revision eb3eb38c009cc70439ad74fa78055d6551981520)
+++ sasmodels/models/fractal_core_shell.py	(revision 2cc8aa2a798267dbefb6309f9f93a254fa180b82)
@@ -134,6 +134,6 @@
     return radius + thickness
 
-tests = [[{'radius': 20.0, 'thickness': 10.0}, 'ER', 30.0],
-
+#tests = [[{'radius': 20.0, 'thickness': 10.0}, 'ER', 30.0],
+tests = [
 #         # At some point the SasView 3.x test result was deemed incorrect. The
           #following tests were verified against NIST IGOR macros ver 7.850.
Index: sasmodels/models/fuzzy_sphere.c
===================================================================
--- sasmodels/models/fuzzy_sphere.c	(revision ee60aa768b29c58e20b19e94578ec3a94a0889d7)
+++ sasmodels/models/fuzzy_sphere.c	(revision ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -0,0 +1,27 @@
+static double form_volume(double radius, double fuzziness)
+{
+    return M_4PI_3*cube(radius);
+}
+
+static double
+effective_radius(int mode, double radius, double fuzziness)
+{
+    switch (mode) {
+    case 1: // radius
+        return radius;
+    case 2: // radius + fuzziness
+        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 ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -78,37 +78,11 @@
               ["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
-
-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
-
-# VR defaults to 1.0
+effective_radius_type = ["radius", "radius + fuzziness"]
 
 def random():
Index: sasmodels/models/hollow_cylinder.c
===================================================================
--- sasmodels/models/hollow_cylinder.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/hollow_cylinder.c	(revision ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -15,11 +15,50 @@
 }
 
+// TODO: interface to form_volume/shell_volume not yet settled
+static double
+shell_volume(double *total, double radius, double thickness, double length)
+{
+    *total = M_PI*length*square(radius+thickness);
+    return *total - M_PI*length*radius*radius;
+}
+
 static double
 form_volume(double radius, double thickness, double length)
 {
-    double v_shell = M_PI*length*(square(radius+thickness) - radius*radius);
-    return v_shell;
+    double total;
+    return shell_volume(&total, radius, thickness, length);
 }
 
+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)
+{
+    switch (mode) {
+    case 1: // equivalent sphere
+        return radius_from_volume(radius, thickness, length);
+    case 2: // outer radius
+        return radius + thickness;
+    case 3: // half length
+        return 0.5*length;
+    case 4: // half outer min dimension
+        return (radius + thickness < 0.5*length ? radius + thickness : 0.5*length);
+    case 5: // half outer max dimension
+        return (radius + thickness > 0.5*length ? radius + thickness : 0.5*length);
+    case 6: // half outer diagonal
+        return radius_from_diagonal(radius,thickness,length);
+    }
+}
 
 static void
Index: sasmodels/models/hollow_cylinder.py
===================================================================
--- sasmodels/models/hollow_cylinder.py	(revision 455aaa1fbc231bb1caf309ca6419f0e3ffd5e19c)
+++ sasmodels/models/hollow_cylinder.py	(revision ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -100,36 +100,9 @@
 source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "hollow_cylinder.c"]
 have_Fq = True
-
-# 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 VR(radius, thickness, length):
-    """
-    :param radius:      Cylinder radius
-    :param thickness:   Cylinder wall thickness
-    :param length:      Cylinder length
-    :return:            Volf ratio for P(q)*S(q)
-    """
-    router = radius + thickness
-    vol_core = pi*radius*radius*length
-    vol_total = pi*router*router*length
-    vol_shell = vol_total - vol_core
-    return vol_total, vol_shell
+effective_radius_type = [
+    "equivalent sphere", "outer radius", "half length",
+    "half outer min dimension", "half outer max dimension",
+    "half outer diagonal",
+    ]
 
 def random():
@@ -162,5 +135,5 @@
 tests = [
     [{}, 0.00005, 1764.926],
-    [{}, 'VR', 0.55555556],
+#    [{}, 'VR', 0.55555556],
     [{}, 0.001, 1756.76],
     [{}, (qx, qy), 2.36885476192],
Index: sasmodels/models/hollow_rectangular_prism.c
===================================================================
--- sasmodels/models/hollow_rectangular_prism.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/hollow_rectangular_prism.c	(revision ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -1,4 +1,5 @@
+// TODO: interface to form_volume/shell_volume not yet settled
 static double
-form_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness)
+shell_volume(double *total, double length_a, double b2a_ratio, double c2a_ratio, double thickness)
 {
     double length_b = length_a * b2a_ratio;
@@ -8,7 +9,37 @@
     double c_core = length_c - 2.0*thickness;
     double vol_core = a_core * b_core * c_core;
-    double vol_total = length_a * length_b * length_c;
-    double vol_shell = vol_total - vol_core;
-    return vol_shell;
+    *total = length_a * length_b * length_c;
+    return *total - vol_core;
+}
+
+static double
+form_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness)
+{
+    double total;
+    return shell_volume(&total, length_a, b2a_ratio, c2a_ratio, thickness);
+}
+
+static double
+effective_radius(int mode, double length_a, double b2a_ratio, double c2a_ratio, double thickness)
+//effective_radius_type = ["equivalent sphere","half length_a", "half length_b", "half length_c",
+//                         "equivalent outer circular cross-section","half ab diagonal","half diagonal"]
+// NOTE length_a is external dimension!
+{
+    switch (mode) {
+    case 1: // equivalent sphere
+        return cbrt(0.75*cube(length_a)*b2a_ratio*c2a_ratio/M_PI);
+    case 2: // half length_a
+        return 0.5 * length_a;
+    case 3: // half length_b
+        return 0.5 * length_a*b2a_ratio;
+    case 4: // half length_c
+        return 0.5 * length_a*c2a_ratio;
+    case 5: // equivalent outer circular cross-section
+        return length_a*sqrt(b2a_ratio/M_PI);
+    case 6: // half ab diagonal
+        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio)));
+    case 7: // half diagonal
+        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 455aaa1fbc231bb1caf309ca6419f0e3ffd5e19c)
+++ sasmodels/models/hollow_rectangular_prism.py	(revision ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -132,5 +132,5 @@
                "Solvent scattering length density"],
               ["length_a", "Ang", 35, [0, inf], "volume",
-               "Shorter side of the parallelepiped"],
+               "Shortest, external, size of the parallelepiped"],
               ["b2a_ratio", "Ang", 1, [0, inf], "volume",
                "Ratio sides b/a"],
@@ -149,33 +149,9 @@
 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.)
-
-def VR(length_a, b2a_ratio, c2a_ratio, thickness):
-    """
-    Return shell volume and total volume
-    """
-    b_side = length_a * b2a_ratio
-    c_side = length_a * c2a_ratio
-    a_core = length_a - 2.0*thickness
-    b_core = b_side - 2.0*thickness
-    c_core = c_side - 2.0*thickness
-    vol_core = a_core * b_core * c_core
-    vol_total = length_a * b_side * c_side
-    vol_shell = vol_total - vol_core
-    return vol_total, vol_shell
-
+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 random():
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 ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -1,9 +1,40 @@
+// TODO: interface to form_volume/shell_volume not yet settled
 static double
-form_volume(double length_a, double b2a_ratio, double c2a_ratio)
+shell_volume(double *total, double length_a, double b2a_ratio, double c2a_ratio)
 {
     double length_b = length_a * b2a_ratio;
     double length_c = length_a * c2a_ratio;
     double vol_shell = 2.0 * (length_a*length_b + length_a*length_c + length_b*length_c);
+    *total = length_a * length_b * length_c;
     return vol_shell;
+}
+
+static double
+form_volume(double length_a, double b2a_ratio, double c2a_ratio)
+{
+    double total;
+    return shell_volume(&total, length_a, b2a_ratio, c2a_ratio);
+}
+
+
+static double
+effective_radius(int mode, double length_a, double b2a_ratio, double c2a_ratio)
+{
+    switch (mode) {
+    case 1: // equivalent sphere
+        return cbrt(0.75*cube(length_a)*b2a_ratio*c2a_ratio/M_PI);
+    case 2: // half length_a
+        return 0.5 * length_a;
+    case 3: // half length_b
+        return 0.5 * length_a*b2a_ratio;
+    case 4: // half length_c
+        return 0.5 * length_a*c2a_ratio;
+    case 5: // equivalent outer circular cross-section
+        return length_a*sqrt(b2a_ratio/M_PI);
+    case 6: // half ab diagonal
+        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio)));
+    case 7: // half diagonal
+        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 6e7d7b66fc1ec444c4e30ac47db7bbcdc554ecc3)
+++ sasmodels/models/hollow_rectangular_prism_thin_walls.py	(revision ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -109,27 +109,9 @@
 source = ["lib/gauss76.c", "hollow_rectangular_prism_thin_walls.c"]
 have_Fq = True
-
-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):
-    """
-        Return shell volume and total volume
-    """
-    b_side = length_a * b2a_ratio
-    c_side = length_a * c2a_ratio
-    vol_total = length_a * b_side * c_side
-    vol_shell = 2.0 * (length_a*b_side + length_a*c_side + b_side*c_side)
-    return vol_shell, vol_total
+effective_radius_type = [
+    "equivalent sphere", "half length_a", "half length_b", "half length_c",
+    "equivalent outer circular cross-section",
+    "half ab diagonal", "half diagonal",
+    ]
 
 
Index: sasmodels/models/mono_gauss_coil.c
===================================================================
--- sasmodels/models/mono_gauss_coil.c	(revision ee60aa768b29c58e20b19e94578ec3a94a0889d7)
+++ sasmodels/models/mono_gauss_coil.c	(revision ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -0,0 +1,60 @@
+static double form_volume(double rg)
+{
+    return 1.0;
+}
+
+static double
+effective_radius(int mode, double rg)
+{
+    switch (mode) {
+    case 1: // R_g
+        return rg;
+    case 2: // 2R_g
+        return 2.0*rg;
+    case 3: // 3R_g
+        return 3.0*rg;
+    case 4: // (5/3)^0.5*R_g
+        return sqrt(5.0/3.0)*rg;
+    }
+}
+
+static double
+gauss_coil(double qr)
+{
+    const double x = qr*qr;
+
+    // Use series expansion at low q for higher accuracy. We could use
+    // smaller polynomials if we sacrifice some digits of precision or
+    // introduce an additional series expansion around x == 1.
+    // See explore/precision.py, gauss_coil function.
+#if FLOAT_SIZE>4 // DOUBLE_PRECISION
+    // For double precision: use O(5) Pade with 0.5 cutoff (10 mad + 1 divide)
+    if (x < 0.5) {
+        // PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 8}]
+        const double A1=1./12., A2=2./99., A3=1./2640., A4=1./23760., A5=-1./1995840.;
+        const double B1=5./12., B2=5./66., B3=1./132., B4=1./2376., B5=1./95040.;
+        return (((((A5*x + A4)*x + A3)*x + A2)*x + A1)*x + 1.)
+                / (((((B5*x + B4)*x + B3)*x + B2)*x + B1)*x + 1.);
+    }
+#else
+    // For single precision: use O(7) Taylor with 0.8 cutoff (7 mad)
+    if (x < 0.8) {
+        const double C0 = +1.;
+        const double C1 = -1./3.;
+        const double C2 = +1./12.;
+        const double C3 = -1./60.;
+        const double C4 = +1./360.;
+        const double C5 = -1./2520.;
+        const double C6 = +1./20160.;
+        const double C7 = -1./181440.;
+        return ((((((C7*x + C6)*x + C5)*x + C4)*x + C3)*x + C2)*x + C1)*x + C0;
+    }
+#endif
+
+    return 2.0 * (expm1(-x) + x)/(x*x);
+}
+
+double Iq(double q, double i_zero, double rg)
+{
+    return i_zero * gauss_coil(q*rg);
+}
Index: sasmodels/models/mono_gauss_coil.py
===================================================================
--- sasmodels/models/mono_gauss_coil.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
+++ sasmodels/models/mono_gauss_coil.py	(revision ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -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,12 @@
 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"],
     ]
 # 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
+source = ["mono_gauss_coil.c"]
+have_Fq = False
+effective_radius_type = ["R_g", "2R_g", "3R_g", "(5/3)^0.5*R_g"]
 
-    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 ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -47,4 +47,11 @@
 }
 
+static double
+effective_radius(int mode, double radius, double thick_shell, double thick_solvent, double fp_n_shells)
+{
+    // case 1: outer radius
+    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 ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -146,8 +146,5 @@
 source = ["lib/sas_3j1x_x.c", "multilayer_vesicle.c"]
 have_Fq = True
-
-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
+effective_radius_type = ["outer radius"]
 
 def random():
Index: sasmodels/models/onion.c
===================================================================
--- sasmodels/models/onion.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/onion.c	(revision ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -30,5 +30,5 @@
 
 static double
-form_volume(double radius_core, double n_shells, double thickness[])
+outer_radius(double radius_core, double n_shells, double thickness[])
 {
   int n = (int)(n_shells+0.5);
@@ -37,5 +37,18 @@
     r += thickness[i];
   }
-  return M_4PI_3*cube(r);
+  return r;
+}
+
+static double
+form_volume(double radius_core, double n_shells, double thickness[])
+{
+  return M_4PI_3*cube(outer_radius(radius_core, n_shells, thickness));
+}
+
+static double
+effective_radius(int mode, double radius_core, double n_shells, double thickness[])
+{
+  // case 1: outer radius
+  return outer_radius(radius_core, n_shells, thickness);
 }
 
Index: sasmodels/models/onion.py
===================================================================
--- sasmodels/models/onion.py	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/onion.py	(revision ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -316,4 +316,5 @@
 single = False
 have_Fq = True
+effective_radius_type = ["outer radius"]
 
 profile_axes = ['Radius (A)', 'SLD (1e-6/A^2)']
@@ -365,9 +366,4 @@
 
     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
 
 demo = {
Index: sasmodels/models/parallelepiped.c
===================================================================
--- sasmodels/models/parallelepiped.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/parallelepiped.c	(revision ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -5,4 +5,24 @@
 }
 
+static double
+effective_radius(int mode, double length_a, double length_b, double length_c)
+{
+    switch (mode) {
+    case 1: // equivalent sphere
+        return cbrt(0.75*length_a*length_b*length_c/M_PI);
+    case 2: // half length_a
+        return 0.5 * length_a;
+    case 3: // half length_b
+        return 0.5 * length_b;
+    case 4: // half length_c
+        return 0.5 * length_c;
+    case 5: // equivalent circular cross-section
+        return sqrt(length_a*length_b/M_PI);
+    case 6: // half ab diagonal
+        return 0.5*sqrt(length_a*length_a + length_b*length_b);
+    case 7: // half diagonal
+        return 0.5*sqrt(length_a*length_a + length_b*length_b + length_c*length_c);
+    }
+}
 
 static void
Index: sasmodels/models/parallelepiped.py
===================================================================
--- sasmodels/models/parallelepiped.py	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/parallelepiped.py	(revision ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -140,5 +140,5 @@
    ensuring that the inequality $A < B < C$ is not violated,  The calculation
    will not report an error, but the results may be not correct.
-   
+
 .. _parallelepiped-orientation:
 
@@ -231,23 +231,8 @@
 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.)
-
-# VR defaults to 1.0
-
+effective_radius_type = [
+    "equivalent sphere", "half length_a", "half length_b", "half length_c",
+    "equivalent circular cross-section", "half ab diagonal", "half diagonal",
+    ]
 
 def random():
Index: sasmodels/models/pearl_necklace.py
===================================================================
--- sasmodels/models/pearl_necklace.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
+++ sasmodels/models/pearl_necklace.py	(revision 2cc8aa2a798267dbefb6309f9f93a254fa180b82)
@@ -140,4 +140,7 @@
             thick_string_pd=0.2, thick_string_pd_n=5,
            )
-
-tests = [[{}, 0.001, 17380.245], [{}, 'ER', 115.39502]]
+# ER function is not being used here, not that it is likely very sensible to 
+# include an S(Q) with this model, the default in sasview 5.0 would be to the 
+# "unconstrained" radius_effective.
+#tests = [[{}, 0.001, 17380.245], [{}, 'ER', 115.39502]]
+tests = [[{}, 0.001, 17380.245]]
Index: sasmodels/models/pringle.c
===================================================================
--- sasmodels/models/pringle.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
+++ sasmodels/models/pringle.c	(revision ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -104,4 +104,15 @@
 }
 
+static double
+effective_radius(int mode, double radius, double thickness, double alpha, double beta)
+{
+    switch (mode) {
+    case 1: // equivalent sphere
+        return cbrt(0.75*radius*radius*thickness);
+    case 2: // radius
+        return radius;
+    }
+}
+
 double Iq(
     double q,
Index: sasmodels/models/pringle.py
===================================================================
--- sasmodels/models/pringle.py	(revision ef07e950965796e2a3a8156fcade98ce5e403001)
+++ sasmodels/models/pringle.py	(revision ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -72,14 +72,7 @@
 
 
-source = ["lib/polevl.c", "lib/sas_J0.c", "lib/sas_J1.c", \
+source = ["lib/polevl.c", "lib/sas_J0.c", "lib/sas_J1.c",
           "lib/sas_JN.c", "lib/gauss76.c", "pringle.c"]
-
-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.)
+effective_radius_type = ["equivalent sphere", "radius"]
 
 def random():
Index: sasmodels/models/raspberry.c
===================================================================
--- sasmodels/models/raspberry.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/raspberry.c	(revision ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -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,15 @@
     double volume=1.0;
     return volume;
+}
+
+static double
+effective_radius(int mode, double radius_lg, double radius_sm, double penetration)
+{
+    switch (mode) {
+    case 1: // radius_large
+        return radius_lg;
+    case 2: // radius_outer
+        return radius_lg + 2.0*radius_sm - penetration;
+    }
 }
 
Index: sasmodels/models/raspberry.py
===================================================================
--- sasmodels/models/raspberry.py	(revision ef07e950965796e2a3a8156fcade98ce5e403001)
+++ sasmodels/models/raspberry.py	(revision ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -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 ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -3,4 +3,25 @@
 {
     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)
+{
+    switch (mode) {
+    case 1: // equivalent sphere
+        return cbrt(0.75*cube(length_a)*b2a_ratio*c2a_ratio/M_PI);
+    case 2: // half length_a
+        return 0.5 * length_a;
+    case 3: // half length_b
+        return 0.5 * length_a*b2a_ratio;
+    case 4: // half length_c
+        return 0.5 * length_a*c2a_ratio;
+    case 5: // equivalent circular cross-section
+        return length_a*sqrt(b2a_ratio/M_PI);
+    case 6: // half ab diagonal
+        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio)));
+    case 7: // half diagonal
+        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 ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -136,17 +136,8 @@
 source = ["lib/gauss76.c", "rectangular_prism.c"]
 have_Fq = True
-
-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.)
+effective_radius_type = [
+    "equivalent sphere", "half length_a", "half length_b", "half length_c",
+    "equivalent circular cross-section", "half ab diagonal", "half diagonal",
+    ]
 
 def random():
Index: sasmodels/models/sphere.c
===================================================================
--- sasmodels/models/sphere.c	(revision ee60aa768b29c58e20b19e94578ec3a94a0889d7)
+++ sasmodels/models/sphere.c	(revision ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -0,0 +1,20 @@
+static double form_volume(double radius)
+{
+    return M_4PI_3*cube(radius);
+}
+
+static double
+effective_radius(int mode, double radius)
+{
+    // case 1: 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 ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -67,30 +67,7 @@
              ]
 
-source = ["lib/sas_3j1x_x.c"]
+source = ["lib/sas_3j1x_x.c","sphere.c"]
 have_Fq = True
-
-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
-
-# VR defaults to 1.0
+effective_radius_type = ["radius"]
 
 def random():
@@ -106,5 +83,5 @@
       "radius": 120., "radius_pd": 0.2, "radius_pd_n":45},
      0.2, 0.228843],
-    [{"radius": 120., "radius_pd": 0.2, "radius_pd_n":45}, "ER", 120.],
-    [{"radius": 120., "radius_pd": 0.2, "radius_pd_n":45}, "VR", 1.],
+#    [{"radius": 120., "radius_pd": 0.2, "radius_pd_n":45}, "ER", 120.],
+#    [{"radius": 120., "radius_pd": 0.2, "radius_pd_n":45}, "VR", 1.],
 ]
Index: sasmodels/models/spherical_sld.c
===================================================================
--- sasmodels/models/spherical_sld.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/spherical_sld.c	(revision ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -1,6 +1,4 @@
-static double form_volume(
-    double fp_n_shells,
-    double thickness[],
-    double interface[])
+static double
+outer_radius(double fp_n_shells, double thickness[], double interface[])
 {
     int n_shells= (int)(fp_n_shells + 0.5);
@@ -9,5 +7,20 @@
         r += thickness[i] + interface[i];
     }
-    return M_4PI_3*cube(r);
+    return r;
+}
+
+static double form_volume(
+    double fp_n_shells,
+    double thickness[],
+    double interface[])
+{
+    return M_4PI_3*cube(outer_radius(fp_n_shells, thickness, interface));
+}
+
+static double
+effective_radius(int mode, double fp_n_shells, double thickness[], double interface[])
+{
+    // case 1: outer radius
+    return outer_radius(fp_n_shells, thickness, interface);
 }
 
Index: sasmodels/models/spherical_sld.py
===================================================================
--- sasmodels/models/spherical_sld.py	(revision 5601947b596b5a6e24ccffd7b45469b40ae22f9d)
+++ sasmodels/models/spherical_sld.py	(revision ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -22,11 +22,11 @@
 
     0: erf($\nu z$)
-    
+
     1: Rpow($z^\nu$)
-    
+
     2: Lpow($z^\nu$)
-    
+
     3: Rexp($-\nu z$)
-    
+
     4: Lexp($-\nu z$)
 
@@ -217,5 +217,5 @@
               ["interface[n_shells]",  "Ang",        50.0,   [0, inf],       "volume", "thickness of the interface"],
               ["shape[n_shells]",      "",           0,      [SHAPES],       "", "interface shape"],
-              ["nu[n_shells]",         "",           2.5,    [0, inf],       "", "interface shape exponent"],
+              ["nu[n_shells]",         "",           2.5,    [1, inf],       "", "interface shape exponent"],
               ["n_steps",              "",           35,     [0, inf],       "", "number of steps in each interface (must be an odd integer)"],
              ]
@@ -224,4 +224,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)']
@@ -269,12 +270,4 @@
 
 
-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
-
-
 demo = {
     "n_shells": 5,
Index: sasmodels/models/triaxial_ellipsoid.c
===================================================================
--- sasmodels/models/triaxial_ellipsoid.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/triaxial_ellipsoid.c	(revision ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -7,4 +7,36 @@
 }
 
+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)
+{
+    switch (mode) {
+    case 1: // equivalent sphere
+        return radius_from_volume(radius_equat_minor,radius_equat_major, radius_polar);
+    case 2: // min radius
+        return radius_from_min_dimension(radius_equat_minor,radius_equat_major, radius_polar);
+    case 3: // max radius
+        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 ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -158,19 +158,5 @@
 source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "triaxial_ellipsoid.c"]
 have_Fq = True
-
-def ER(radius_equat_minor, radius_equat_major, radius_polar):
-    """
-    Returns the effective radius used in the S*P calculation
-    """
-    from .ellipsoid import ER as ellipsoid_ER
-
-    # now that radii can be in any size order, radii need sorting a,b,c
-    # where a~b and c is either much smaller or much larger
-    radii = np.vstack((radius_equat_major, radius_equat_minor, radius_polar))
-    radii = np.sort(radii, axis=0)
-    selector = (radii[1] - radii[0]) > (radii[2] - radii[1])
-    polar = np.where(selector, radii[0], radii[2])
-    equatorial = np.sqrt(np.where(~selector, radii[0]*radii[1], radii[1]*radii[2]))
-    return ellipsoid_ER(polar, equatorial)
+effective_radius_type = ["equivalent sphere", "min radius", "max radius"]
 
 def random():
Index: sasmodels/models/vesicle.c
===================================================================
--- sasmodels/models/vesicle.c	(revision 71b751d6a3a4ef4164f3b322841ddf3875c025c5)
+++ sasmodels/models/vesicle.c	(revision ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -1,9 +1,25 @@
+// TODO: interface to form_volume/shell_volume not yet settled
+static double
+shell_volume(double *total, double radius, double thickness)
+{
+    //note that for the vesicle model, the volume is ONLY the shell volume
+    *total = M_4PI_3 * cube(radius+thickness);
+    return *total - M_4PI_3*cube(radius);
+}
+
 static double
 form_volume(double radius, double thickness)
 {
     //note that for the vesicle model, the volume is ONLY the shell volume
-    return M_4PI_3*(cube(radius+thickness) - cube(radius));
+    double total;
+    return shell_volume(&total, radius, thickness);
 }
 
+static double
+effective_radius(int mode, double radius, double thickness)
+{
+    // case 1: outer radius
+    return radius + thickness;
+}
 
 static void
Index: sasmodels/models/vesicle.py
===================================================================
--- sasmodels/models/vesicle.py	(revision b4776059692723f870409527b8054b2a8a5e6797)
+++ sasmodels/models/vesicle.py	(revision ee60aa768b29c58e20b19e94578ec3a94a0889d7)
@@ -100,28 +100,5 @@
 source = ["lib/sas_3j1x_x.c", "vesicle.c"]
 have_Fq = True
-
-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):
-    '''
-    returns the volumes of the shell and of the whole sphere including the
-    core plus shell - is used to normalize when including polydispersity.
-
-    :param radius: core radius
-    :param thickness: shell thickness
-    :return whole: volume of core and shell
-    :return whole-core: volume of the shell
-    '''
-
-    whole = 4./3. * pi * (radius + thickness)**3
-    core = 4./3. * pi * radius**3
-    return whole, whole - core
+effective_radius_type = ["outer radius"]
 
 def random():
@@ -151,5 +128,5 @@
          [{}, 0.100600200401, 1.77063682331],
          [{}, 0.5, 0.00355351388906],
-         [{}, 'ER', 130.],
-         [{}, 'VR', 0.54483386436],
+#         [{}, 'ER', 130.],
+#         [{}, 'VR', 0.54483386436],
         ]
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 c952e591d03e500c509f1a9ba5a1b76e32985aa5)
@@ -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
