Index: compare-new.py
===================================================================
--- compare-new.py	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
+++ compare-new.py	(revision f4cf580e329add99780558498dbc8a8ab076ad24)
@@ -3,4 +3,5 @@
 
 import sys
+import math
 
 import numpy as np
@@ -15,5 +16,8 @@
     Load a sasview model given the model name.
     """
-    modelname = modelname
+    # convert model parameters from sasmodel form to sasview form
+    #print "old",sorted(pars.items())
+    modelname, pars = revert_model(modelname, pars)
+    #print "new",sorted(pars.items())
     sans = __import__('sans.models.'+modelname)
     ModelClass = getattr(getattr(sans.models,modelname,None),modelname,None)
@@ -41,5 +45,5 @@
     return kernel
 
-def load_dll(modelname, dtype='single'):
+def load_ctypes(modelname, dtype='single'):
     sasmodels = __import__('sasmodels.models.'+modelname)
     module = getattr(sasmodels.models, modelname, None)
@@ -49,15 +53,25 @@
 def randomize(p, v):
     """
-    Randomizing pars using sasview names.
-    Guessing at angles and slds.
-    """
-    if any(p.endswith(s) for s in ('_pd_n','_pd_nsigmas','_pd_type')):
+    Randomizing parameter.
+
+    Guess the parameter type from name.
+    """
+    if any(p.endswith(s) for s in ('_pd_n','_pd_nsigma','_pd_type')):
         return v
-    if any(s in p for s in ('theta','phi','psi')):
-        return np.random.randint(0,90)
-    elif any(s in p for s in ('sld','Sld','SLD')):
-        return np.random.rand()*1e-5
+    elif any(s in p for s in ('theta','phi','psi')):
+        # orientation in [-180,180], orientation pd in [0,45]
+        if p.endswith('_pd'):
+            return 45*np.random.rand()
+        else:
+            return 360*np.random.rand() - 180
+    elif 'sld' in p:
+        # sld in in [-0.5,10]
+        return 10.5*np.random.rand() - 0.5
+    elif p.endswith('_pd'):
+        # length pd in [0,1]
+        return np.random.rand()
     else:
-        return np.random.rand()
+        # length, scale, background in [0,200]
+        return 200*np.random.rand()
 
 def parlist(pars):
@@ -74,83 +88,89 @@
         if p.endswith("_pd"): pars[p] = 0
 
-def compare(gpuname, gpupars, Ncpu, Ngpu, opts):
-    #from sasmodels.core import load_data
-    #data = load_data('December/DEC07098.DAT')
+def compare(name, pars, Ncpu, Ngpu, opts, set_pars):
+
+    # Sort out data
     qmax = 1.0 if '-highq' in opts else (0.2 if '-midq' in opts else 0.05)
     if "-1d" in opts:
         from sasmodels.bumps_model import empty_data1D
-        qmax = np.log10(qmax)
+        qmax = math.log10(qmax)
         data = empty_data1D(np.logspace(qmax-3, qmax, 128))
+        index = slice(None, None)
     else:
         from sasmodels.bumps_model import empty_data2D, set_beam_stop
         data = empty_data2D(np.linspace(-qmax, qmax, 128))
         set_beam_stop(data, 0.004)
+        index = ~data.mask
     is2D = hasattr(data, 'qx_data')
+
+    # modelling accuracy is determined by dtype and cutoff
     dtype = 'double' if '-double' in opts else 'single'
     cutoff_opts = [s for s in opts if s.startswith('-cutoff')]
     cutoff = float(cutoff_opts[0].split('=')[1]) if cutoff_opts else 1e-5
 
-
-    if '-random' in opts:
-        gpupars = dict((p,randomize(p,v)) for p,v in gpupars.items())
-    if '-pars' in opts: print "pars",parlist(gpupars)
-    if '-mono' in opts: suppress_pd(gpupars)
-
-    cpuname, cpupars = revert_model(gpuname, gpupars)
-
-    try:
-        gpumodel = load_opencl(gpuname, dtype=dtype)
-    except Exception,exc:
-        print exc
-        print "... trying again with single precision"
-        gpumodel = load_opencl(gpuname, dtype='single')
-    model = BumpsModel(data, gpumodel, cutoff=cutoff, **gpupars)
+    # randomize parameters
+    random_opts = [s for s in opts if s.startswith('-random')]
+    if random_opts:
+        if '=' in random_opts[0]:
+            seed = int(random_opts[0].split('=')[1])
+        else:
+            seed = int(np.random.rand()*10000)
+        np.random.seed(seed)
+        print "Randomize using -random=%i"%seed
+        # Note: the sort guarantees order of calls to random number generator
+        pars = dict((p,randomize(p,v)) for p,v in sorted(pars.items()))
+        # The capped cylinder model has a constraint on its parameters
+        if name == 'capped_cylinder' and pars['cap_radius'] < pars['radius']:
+            pars['radius'],pars['cap_radius'] = pars['cap_radius'],pars['radius']
+    pars.update(set_pars)
+
+    # parameter selection
+    if '-mono' in opts:
+        suppress_pd(pars)
+    if '-pars' in opts:
+        print "pars",parlist(pars)
+
+    # OpenCl calculation
     if Ngpu > 0:
+        try:
+            model = load_opencl(name, dtype=dtype)
+        except Exception,exc:
+            print exc
+            print "... trying again with single precision"
+            model = load_opencl(name, dtype='single')
+        problem = BumpsModel(data, model, cutoff=cutoff, **pars)
         toc = tic()
-        for i in range(Ngpu):
+        for _ in range(Ngpu):
             #pars['scale'] = np.random.rand()
-            model.update()
-            gpu = model.theory()
+            problem.update()
+            gpu = problem.theory()
         gpu_time = toc()*1000./Ngpu
-        print "ocl t=%.1f ms, intensity=%f"%(gpu_time, sum(gpu[~np.isnan(gpu)]))
+        print "opencl t=%.1f ms, intensity=%.0f"%(gpu_time, sum(gpu[index]))
         #print max(gpu), min(gpu)
 
-    comp = None
-    if Ncpu > 0 and "-dll" in opts:
-        dllmodel = load_dll(gpuname, dtype=dtype)
-        model = BumpsModel(data, dllmodel, cutoff=cutoff, **gpupars)
+    # ctypes/sasview calculation
+    if Ncpu > 0 and "-ctypes" in opts:
+        model = load_ctypes(name, dtype=dtype)
+        problem = BumpsModel(data, model, cutoff=cutoff, **pars)
         toc = tic()
-        for i in range(Ncpu):
-            model.update()
-            cpu = model.theory()
+        for _ in range(Ncpu):
+            problem.update()
+            cpu = problem.theory()
         cpu_time = toc()*1000./Ncpu
-        comp = "dll"
-
-    elif 0: # Hack to check new vs old for GpuCylinder
-        from Models.code_cylinder_f import GpuCylinder as oldgpu
-        from sasmodel import SasModel
-        oldmodel = SasModel(data, oldgpu, dtype=dtype, **cpupars)
+        comp = "ctypes"
+        print "ctypes t=%.1f ms, intensity=%.0f"%(cpu_time, sum(cpu[index]))
+    elif Ncpu > 0:
+        model = sasview_model(name, **pars)
         toc = tic()
-        for i in range(Ngpu):
-            oldmodel.update()
-            cpu = oldmodel.theory()
-        cpu_time = toc()*1000./Ngpu
-        comp = "old"
-
-    elif Ncpu > 0:
-        cpumodel = sasview_model(cpuname, **cpupars)
-        toc = tic()
-        if is2D:
-            for i in range(Ncpu):
-                cpu = cpumodel.evalDistribution([data.qx_data, data.qy_data])
-        else:
-            for i in range(Ncpu):
-                cpu = cpumodel.evalDistribution(data.x)
+        for _ in range(Ncpu):
+            if is2D:
+                cpu = model.evalDistribution([data.qx_data, data.qy_data])
+            else:
+                cpu = model.evalDistribution(data.x)
         cpu_time = toc()*1000./Ncpu
-        comp = 'sasview'
-        #print "sasview t=%.1f ms, intensity=%.0f"%(cpu_time, sum(cpu[model.index]))
-
-    if comp:
-        print "%s t=%.1f ms, intensity=%f"%(comp, cpu_time, sum(cpu[model.index]))
+        comp = "sasview"
+        print "sasview t=%.1f ms, intensity=%.0f"%(cpu_time, sum(cpu[index]))
+
+    # Compare, but only if computing both forms
     if Ngpu > 0 and Ncpu > 0:
         #print "speedup %.2g"%(cpu_time/gpu_time)
@@ -158,11 +178,11 @@
         #cpu *= max(gpu/cpu)
         resid, relerr = np.zeros_like(gpu), np.zeros_like(gpu)
-        resid[model.index] = (gpu - cpu)[model.index]
-        relerr[model.index] = resid[model.index]/cpu[model.index]
-        print "max(|ocl-%s|)"%comp, max(abs(resid[model.index]))
-        print "max(|(ocl-%s)/ocl|)"%comp, max(abs(relerr[model.index]))
-
+        resid[index] = (gpu - cpu)[index]
+        relerr[index] = resid[index]/cpu[index]
+        print "max(|ocl-%s|)"%comp, max(abs(resid[index]))
+        print "max(|(ocl-%s)/ocl|)"%comp, max(abs(relerr[index]))
+
+    # Plot if requested
     if '-noplot' in opts: return
-
     import matplotlib.pyplot as plt
     if Ncpu > 0:
@@ -173,10 +193,12 @@
         if Ncpu > 0: plt.subplot(132)
         plot_data(data, gpu, scale='log')
-        plt.title("ocl t=%.1f ms"%gpu_time)
+        plt.title("opencl t=%.1f ms"%gpu_time)
     if Ncpu > 0 and Ngpu > 0:
         plt.subplot(133)
         err = resid if '-abs' in opts else relerr
+        errstr = "abs err" if '-abs' in opts else "rel err"
+        #err,errstr = gpu/cpu,"ratio"
         plot_data(data, err, scale='linear')
-        plt.title("max rel err = %.3g"%max(abs(err[model.index])))
+        plt.title("max %s = %.3g"%(errstr, max(abs(err[index]))))
         if is2D: plt.colorbar()
     plt.show()
@@ -200,7 +222,7 @@
     -lowq/-midq/-highq use q values up to 0.05, 0.2 or 1.0
     -2d/-1d uses 1d or 2d random data
-    -preset/-random randomizes the parameters
+    -preset/-random[=seed] preset or random parameters
     -poly/-mono force monodisperse/polydisperse
-    -sasview/-dll whether cpu is tested using sasview or dll
+    -sasview/-ctypes whether cpu is tested using sasview or ctypes
     -cutoff=1e-5/value cutoff for including a point in polydispersity
     -nopars/-pars prints the parameter set or not
@@ -214,4 +236,16 @@
     %s
 """
+
+VALID_OPTIONS = [
+    'plot','noplot',
+    'single','double',
+    'lowq','midq','highq',
+    '2d','1d',
+    'preset','random',
+    'poly','mono',
+    'sasview','ctypes',
+    'nopars','pars',
+    'rel','abs',
+    ]
 
 def main():
@@ -227,4 +261,13 @@
         sys.exit(1)
 
+    valid_opts = set(VALID_OPTIONS)
+    invalid = [o[1:] for o in opts
+               if o[1:] not in valid_opts
+                  and not o.startswith('-cutoff=')
+                  and not o.startswith('-random=')]
+    if invalid:
+        print "Invalid options: %s"%(", ".join(invalid))
+        sys.exit(1)
+
     name, pars = MODELS[args[0]]()
     Nopencl = int(args[1]) if len(args) > 1 else 5
@@ -238,15 +281,15 @@
 
     # Fill in parameters given on the command line
+    set_pars = {}
     for arg in args[3:]:
         k,v = arg.split('=')
         if k not in pars:
             # extract base name without distribution
-            # style may be either a.d or a_pd_d
-            s = set((p.split('_pd')[0]).split('.')[0] for p in pars)
+            s = set(p.split('_pd')[0] for p in pars)
             print "%r invalid; parameters are: %s"%(k,", ".join(sorted(s)))
             sys.exit(1)
-        pars[k] = float(v) if not v.endswith('type') else v
-
-    compare(name, pars, Nsasview, Nopencl, opts)
+        set_pars[k] = float(v) if not v.endswith('type') else v
+
+    compare(name, pars, Nsasview, Nopencl, opts, set_pars)
 
 # ===========================================================================
@@ -265,11 +308,11 @@
     pars = dict(
         scale=1, background=0,
-        sld=6e-6, solvent_sld=1e-6,
+        sld=6, solvent_sld=1,
         #radius=5, length=20,
         radius=260, length=290,
-        theta=30, phi=15,
+        theta=30, phi=0,
         radius_pd=.2, radius_pd_n=1,
         length_pd=.2,length_pd_n=1,
-        theta_pd=15, theta_pd_n=25,
+        theta_pd=15, theta_pd_n=45,
         phi_pd=15, phi_pd_n=1,
         )
@@ -280,11 +323,11 @@
     pars = dict(
         scale=1, background=0,
-        sld=6e-6, solvent_sld=1e-6,
-        radius=260, cap_radius=80000, length=290,
+        sld=6, solvent_sld=1,
+        radius=260, cap_radius=290, length=290,
         theta=30, phi=15,
         radius_pd=.2, radius_pd_n=1,
         cap_radius_pd=.2, cap_radius_pd_n=1,
         length_pd=.2, length_pd_n=1,
-        theta_pd=15, theta_pd_n=25,
+        theta_pd=15, theta_pd_n=45,
         phi_pd=15, phi_pd_n=1,
         )
@@ -296,12 +339,12 @@
     pars = dict(
         scale=1, background=0,
-        core_sld=6e-6, shell_sld=8e-6, solvent_sld=1e-6,
-        radius=325, thickness=25, length=34.2709,
-        theta=90, phi=0,
-        radius_pd=0.1, radius_pd_n=10,
-        length_pd=0.1, length_pd_n=5,
-        thickness_pd=0.1, thickness_pd_n=5,
-        theta_pd=15.8, theta_pd_n=5,
-        phi_pd=0.0008748, phi_pd_n=5,
+        core_sld=6, shell_sld=8, solvent_sld=1,
+        radius=45, thickness=25, length=340,
+        theta=30, phi=15,
+        radius_pd=.2, radius_pd_n=1,
+        length_pd=.2, length_pd_n=1,
+        thickness_pd=.2, thickness_pd_n=1,
+        theta_pd=15, theta_pd_n=45,
+        phi_pd=15, phi_pd_n=1,
         )
     return 'core_shell_cylinder', pars
@@ -312,11 +355,11 @@
     pars = dict(
         scale=1, background=0,
-        sld=6e-6, solvent_sld=1e-6,
+        sld=6, solvent_sld=1,
         rpolar=50, requatorial=30,
-        theta=0, phi=0,
-        rpolar_pd=0.3, rpolar_pd_n=10,
-        requatorial_pd=0, requatorial_pd_n=10,
-        theta_pd=0, theta_pd_n=45,
-        phi_pd=0, phi_pd_n=45,
+        theta=30, phi=15,
+        rpolar_pd=.2, rpolar_pd_n=1,
+        requatorial_pd=.2, requatorial_pd_n=1,
+        theta_pd=15, theta_pd_n=45,
+        phi_pd=15, phi_pd_n=1,
         )
     return 'ellipsoid', pars
@@ -327,13 +370,13 @@
     pars = dict(
         scale=1, background=0,
-        sld=6e-6, solvent_sld=1e-6,
-        theta=0, phi=0, psi=0,
+        sld=6, solvent_sld=1,
+        theta=30, phi=15, psi=5,
         req_minor=25, req_major=36, rpolar=50,
-        theta_pd=0, theta_pd_n=5,
-        phi_pd=0, phi_pd_n=5,
-        psi_pd=0, psi_pd_n=5,
-        req_minor_pd=0, req_minor_pd_n=5,
-        req_major_pd=0, req_major_pd_n=5,
-        rpolar_pd=.3, rpolar_pd_n=25,
+        req_minor_pd=0, req_minor_pd_n=1,
+        req_major_pd=0, req_major_pd_n=1,
+        rpolar_pd=.2, rpolar_pd_n=1,
+        theta_pd=15, theta_pd_n=45,
+        phi_pd=15, phi_pd_n=1,
+        psi_pd=15, psi_pd_n=1,
         )
     return 'triaxial_ellipsoid', pars
@@ -343,7 +386,7 @@
     pars = dict(
         scale=1, background=0,
-        sld=6e-6, solvent_sld=1e-6,
+        sld=6, solvent_sld=1,
         radius=120,
-        radius_pd=.3, radius_pd_n=5,
+        radius_pd=.2, radius_pd_n=45,
         )
     return 'sphere', pars
@@ -353,7 +396,7 @@
     pars = dict(
         scale=1, background=0,
-        sld=6e-6, solvent_sld=1e-6,
+        sld=6, solvent_sld=1,
         thickness=40,
-        thickness_pd= 0.3, thickness_pd_n=40,
+        thickness_pd= 0.2, thickness_pd_n=40,
         )
     return 'lamellar', pars
Index: sasmodels/convert.py
===================================================================
--- sasmodels/convert.py	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
+++ sasmodels/convert.py	(revision f4cf580e329add99780558498dbc8a8ab076ad24)
@@ -37,5 +37,5 @@
         name='sphere',
         sldSph='sld', sldSolv='solvent_sld',
-        #radius='radius',  # DO NOT LIST IDENTICAL PARAMETERS!
+        radius='radius',  # listing identical parameters is optional
         ),
     }
@@ -59,4 +59,5 @@
     newpars = pars.copy()
     for old,new in mapping.items():
+        if old == new: continue
         # Bumps style parameter names
         for variant in ("", "_pd", "_pd_n", "_pd_nsigma", "_pd_type"):
Index: sasmodels/gen.py
===================================================================
--- sasmodels/gen.py	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
+++ sasmodels/gen.py	(revision f4cf580e329add99780558498dbc8a8ab076ad24)
@@ -352,5 +352,5 @@
   const real I = %(fn)s(%(qcall)s, %(pcall)s);
   if (I>=REAL(0.0)) { // scattering cannot be negative
-    ret += weight*I;
+    ret += weight*I%(sasview_spherical)s;
     norm += weight;
     %(volume_norm)s
@@ -364,5 +364,10 @@
 SPHERICAL_CORRECTION="""\
 // Correction factor for spherical integration p(theta) I(q) sin(theta) dtheta
-real spherical_correction = (Ntheta>1 ? fabs(cos(M_PI_180*phi)) : REAL(1.0));\
+real spherical_correction = (Ntheta>1 ? fabs(sin(M_PI_180*theta)) : REAL(1.0));\
+"""
+# Use this to reproduce sasview behaviour
+SASVIEW_SPHERICAL_CORRECTION="""\
+// Correction factor for spherical integration p(theta) I(q) sin(theta) dtheta
+real spherical_correction = (Ntheta>1 ? fabs(cos(M_PI_180*theta))*M_PI_2 : REAL(1.0));\
 """
 
@@ -472,10 +477,16 @@
     fq_pars = [p[0] for p in info['parameters'][len(COMMON_PARAMETERS):]
                if p[0] in set(fixed_pars+pd_pars)]
-    if False and "phi" in pd_pars:
+    if False and "theta" in pd_pars:
         spherical_correction = [indent(SPHERICAL_CORRECTION, depth)]
         weights = [p+"_w" for p in pd_pars]+['spherical_correction']
+        sasview_spherical = ""
+    elif "theta" in pd_pars:
+        spherical_correction = [indent(SASVIEW_SPHERICAL_CORRECTION,depth)]
+        weights = [p+"_w" for p in pd_pars]
+        sasview_spherical = "*spherical_correction"
     else:
         spherical_correction = []
         weights = [p+"_w" for p in pd_pars]
+        sasview_spherical = ""
     subst = {
         'weight_product': "*".join(weights),
@@ -484,4 +495,5 @@
         'qcall': q_pars['qcall'],
         'pcall': ", ".join(fq_pars), # skip scale and background
+        'sasview_spherical': sasview_spherical,
         }
     loop_body = [indent(LOOP_BODY%subst, depth)]
Index: sasmodels/gpu.py
===================================================================
--- sasmodels/gpu.py	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
+++ sasmodels/gpu.py	(revision f4cf580e329add99780558498dbc8a8ab076ad24)
@@ -22,6 +22,4 @@
 devices, where it can be combined with other structure factors and form
 factors and have instrumental resolution effects applied.
-
-
 """
 import numpy as np
@@ -59,5 +57,5 @@
     source, info = gen.make(kernel_module)
     ## for debugging, save source to a .cl file, edit it, and reload as model
-    open(info['name']+'.cl','w').write(source)
+    #open(info['name']+'.cl','w').write(source)
     #source = open(info['name']+'.cl','r').read()
     return GpuModel(source, info, dtype)
Index: sasmodels/models/capped_cylinder.c
===================================================================
--- sasmodels/models/capped_cylinder.c	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
+++ sasmodels/models/capped_cylinder.c	(revision f4cf580e329add99780558498dbc8a8ab076ad24)
@@ -22,7 +22,14 @@
     const real upper = REAL(1.0);
     const real lower = h/cap_radius; // integral lower bound
-    const real m = q*cos_alpha*cap_radius; // cos argument slope
-    const real b = q*cos_alpha*(REAL(0.5)*length-h); // cos argument intercept
-    const real qrst = q*sin_alpha*cap_radius; // Q*R*sin(theta)
+    // cos term in integral is:
+    //    cos (q (R t - h + L/2) cos(alpha))
+    // so turn it into:
+    //    cos (m t + b)
+    // where:
+    //    m = q R cos(alpha)
+    //    b = q(L/2-h) cos(alpha)
+    const real m = q*cap_radius*cos_alpha; // cos argument slope
+    const real b = q*(REAL(0.5)*length-h)*cos_alpha; // cos argument intercept
+    const real qrst = q*cap_radius*sin_alpha; // Q*R*sin(theta)
     real total = REAL(0.0);
     for (int i=0; i<76 ;i++) {
@@ -31,6 +38,6 @@
         const real t = REAL(0.5)*(Gauss76Z[i]*(upper-lower)+upper+lower);
         const real radical = REAL(1.0) - t*t;
-        const real caparg = qrst*sqrt(radical); // cap bessel function arg
-        const real be = (caparg == REAL(0.0) ? REAL(0.5) : J1(caparg)/caparg);
+        const real arg = qrst*sqrt(radical); // cap bessel function arg
+        const real be = (arg == REAL(0.0) ? REAL(0.5) : J1(arg)/arg);
         const real Fq = cos(m*t + b) * radical * be;
         total += Gauss76Wt[i] * Fq;
@@ -84,5 +91,5 @@
         // faster since we don't multiply and divide sin(alpha).
         const real besarg = q*radius*sn;
-        const real siarg = REAL(0.5)*q*length*cn;
+        const real siarg = q*REAL(0.5)*length*cn;
         // lim_{x->0} J1(x)/x = 1/2,   lim_{x->0} sin(x)/x = 1
         const real bj = (besarg == REAL(0.0) ? REAL(0.5) : J1(besarg)/besarg);
@@ -131,9 +138,6 @@
     const real cap_Fq = _cap_kernel(q, h, cap_radius, length, sn, cn);
 
-    // The following is CylKernel() / sin(alpha), but we are doing it in place
-    // to avoid sin(alpha)/sin(alpha) for alpha = 0.  It is also a teensy bit
-    // faster since we don't multiply and divide sin(alpha).
     const real besarg = q*radius*sn;
-    const real siarg = REAL(0.5)*q*length*cn;
+    const real siarg = q*REAL(0.5)*length*cn;
     // lim_{x->0} J1(x)/x = 1/2,   lim_{x->0} sin(x)/x = 1
     const real bj = (besarg == REAL(0.0) ? REAL(0.5) : J1(besarg)/besarg);
Index: sasmodels/models/capped_cylinder.py
===================================================================
--- sasmodels/models/capped_cylinder.py	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
+++ sasmodels/models/capped_cylinder.py	(revision f4cf580e329add99780558498dbc8a8ab076ad24)
@@ -1,5 +1,5 @@
 r"""
 Calculates the scattering from a cylinder with spherical section end-caps.
-This model simply becomes the `convex_lens`_ when the length of the cylinder
+This model simply becomes the a convex lens when the length of the cylinder
 $L=0$, that is, a sphereocylinder with end caps that have a radius larger
 than that of the cylinder and the center of the end cap radius lies within
@@ -128,4 +128,13 @@
     [ "radius", "Ang",  20, [0, inf], "volume",
       "Cylinder radius" ],
+    # TODO: use an expression for cap radius with fixed bounds.
+    # The current form requires cap radius R bigger than cylinder radius r.
+    # Could instead use R/r in [1,inf], r/R in [0,1], or the angle between
+    # cylinder and cap in [0,90].  The problem is similar for the barbell
+    # model.  Propose r/R in [0,1] in both cases, with the model specifying
+    # cylinder radius in the capped cylinder model and sphere radius in the
+    # barbell model.  This leads to the natural value of zero for no cap
+    # in the capped cylinder, and zero for no bar in the barbell model.  In
+    # both models, one would be a pill.
     [ "cap_radius", "Ang", 20, [0, inf], "volume",
       "Cap radius" ],
Index: sasmodels/models/sphere.py
===================================================================
--- sasmodels/models/sphere.py	(revision 5d4777d1f3042abc2b98b66609abf23d3bca1828)
+++ sasmodels/models/sphere.py	(revision f4cf580e329add99780558498dbc8a8ab076ad24)
@@ -30,10 +30,10 @@
 
 Iq = """
+    const real qr = q*radius;
     real sn, cn;
-    const real qr = q*radius;
     SINCOS(qr, sn, cn);
-    const real bes = (qr==REAL(0.0) ? REAL(1.0) : (REAL(3.0)*(sn-qr*cn))/(qr*qr*qr));
-    const real f = bes * (sld - solvent_sld) * form_volume(radius);
-    return REAL(1.0e-4)*f*f;
+    const real bes = qr==REAL(0.0) ? REAL(1.0) : REAL(3.0)*(sn-qr*cn)/(qr*qr*qr);
+    const real fq = bes * (sld - solvent_sld) * form_volume(radius);
+    return REAL(1.0e-4)*fq*fq;
     """
 
