Index: sasmodels/compare.py
===================================================================
--- sasmodels/compare.py	(revision 31df0c91cf57174acd8e9cbe5f9caadffdc880ba)
+++ sasmodels/compare.py	(revision 8f04da467b56e4e22291910958adebc417b57c19)
@@ -763,4 +763,6 @@
     for k in range(opts['sets']):
         opts['pars'] = parse_pars(opts)
+        if opts['pars'] is None:
+            return
         result = run_models(opts, verbose=True)
         if opts['plot']:
@@ -1395,4 +1397,6 @@
         elif opts['explore']:
             opts['pars'] = parse_pars(opts)
+            if opts['pars'] is None:
+                return
             explore(opts)
         else:
Index: sasmodels/models/adsorbed_layer.py
===================================================================
--- sasmodels/models/adsorbed_layer.py	(revision 0bdddc2288f8fc5bc5e416af84e456e811c3b733)
+++ sasmodels/models/adsorbed_layer.py	(revision 8f04da467b56e4e22291910958adebc417b57c19)
@@ -99,5 +99,5 @@
     #    scale = 6 pi/100 (contrast/density*absorbed_amount)^2 * Vf/radius
     # the remaining parameters can be randomly generated from zero to
-    # twice the default value.
+    # twice the default value as done by default in compare.py
     import numpy as np
     pars = dict(
Index: sasmodels/models/bcc_paracrystal.py
===================================================================
--- sasmodels/models/bcc_paracrystal.py	(revision 1511c37c7ea796eff9bd41e33beabb6efce39d91)
+++ sasmodels/models/bcc_paracrystal.py	(revision 8f04da467b56e4e22291910958adebc417b57c19)
@@ -154,8 +154,7 @@
 # april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct!
 # add 2d test later
-q =4.*pi/220.
+q = 4.*pi/220.
 tests = [
-    [{ },
-     [0.001, q, 0.215268], [1.46601394721, 2.85851284174, 0.00866710287078]],
+    [{}, [0.001, q, 0.215268], [1.46601394721, 2.85851284174, 0.00866710287078]],
     [{'theta': 20.0, 'phi': 30, 'psi': 40.0}, (-0.017, 0.035), 2082.20264399],
     [{'theta': 20.0, 'phi': 30, 'psi': 40.0}, (-0.081, 0.011), 0.436323144781],
Index: sasmodels/models/be_polyelectrolyte.py
===================================================================
--- sasmodels/models/be_polyelectrolyte.py	(revision a151caa16f645cbb8ff93c008adabc6d6a7f5136)
+++ sasmodels/models/be_polyelectrolyte.py	(revision 8f04da467b56e4e22291910958adebc417b57c19)
@@ -141,4 +141,5 @@
 def random():
     import numpy as np
+    # TODO: review random be_polyelectrolyte model generation
     pars = dict(
         scale=10000, #background=0,
Index: sasmodels/models/binary_hard_sphere.py
===================================================================
--- sasmodels/models/binary_hard_sphere.py	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
+++ sasmodels/models/binary_hard_sphere.py	(revision 8f04da467b56e4e22291910958adebc417b57c19)
@@ -110,4 +110,24 @@
 source = ["lib/sas_3j1x_x.c", "binary_hard_sphere.c"]
 
+def random():
+    import numpy as np
+    # TODO: binary_hard_sphere fails at low qr
+    radius_lg = 10**np.random.uniform(2, 4.7)
+    radius_sm = 10**np.random.uniform(2, 4.7)
+    volfraction_lg = 10**np.random.uniform(-3, -0.3)
+    volfraction_sm = 10**np.random.uniform(-3, -0.3)
+    # TODO: Get slightly different results if large and small are swapped
+    # modify the model so it doesn't care which is which
+    if radius_lg < radius_sm:
+        radius_lg, radius_sm = radius_sm, radius_lg
+        volfraction_lg, volfraction_sm = volfraction_sm, volfraction_lg
+    pars = dict(
+        radius_lg=radius_lg,
+        radius_sm=radius_sm,
+        volfraction_lg=volfraction_lg,
+        volfraction_sm=volfraction_sm,
+    )
+    return pars
+
 # parameters for demo and documentation
 demo = dict(sld_lg=3.5, sld_sm=0.5, sld_solvent=6.36,
Index: sasmodels/models/core_shell_bicelle_elliptical.py
===================================================================
--- sasmodels/models/core_shell_bicelle_elliptical.py	(revision 9802ab3945a154b871264715dcb463d4ec644afa)
+++ sasmodels/models/core_shell_bicelle_elliptical.py	(revision 8f04da467b56e4e22291910958adebc417b57c19)
@@ -6,5 +6,5 @@
 core-shell scattering length density profile. Thus this is a variation
 of the core-shell bicelle model, but with an elliptical cylinder for the core.
-Outer shells on the rims and flat ends may be of different thicknesses and 
+Outer shells on the rims and flat ends may be of different thicknesses and
 scattering length densities. The form factor is normalized by the total particle volume.
 
@@ -17,6 +17,6 @@
 .. figure:: img/core_shell_bicelle_parameters.png
 
-   Cross section of model used here. Users will have 
-   to decide how to distribute "heads" and "tails" between the rim, face 
+   Cross section of model used here. Users will have
+   to decide how to distribute "heads" and "tails" between the rim, face
    and core regions in order to estimate appropriate starting parameters.
 
@@ -27,6 +27,6 @@
 .. math::
 
-    \rho(r) = 
-      \begin{cases} 
+    \rho(r) =
+      \begin{cases}
       &\rho_c \text{ for } 0 \lt r \lt R; -L \lt z\lt L \\[1.5ex]
       &\rho_f \text{ for } 0 \lt r \lt R; -(L+2t) \lt z\lt -L;
@@ -48,11 +48,11 @@
 .. math::
 
-        \begin{align}    
-    F(Q,\alpha,\psi) = &\bigg[ 
+        \begin{align}
+    F(Q,\alpha,\psi) = &\bigg[
     (\rho_c - \rho_f) V_c \frac{2J_1(QR'sin \alpha)}{QR'sin\alpha}\frac{sin(QLcos\alpha/2)}{Q(L/2)cos\alpha} \\
     &+(\rho_f - \rho_r) V_{c+f} \frac{2J_1(QR'sin\alpha)}{QR'sin\alpha}\frac{sin(Q(L/2+t_f)cos\alpha)}{Q(L/2+t_f)cos\alpha} \\
     &+(\rho_r - \rho_s) V_t \frac{2J_1(Q(R'+t_r)sin\alpha)}{Q(R'+t_r)sin\alpha}\frac{sin(Q(L/2+t_f)cos\alpha)}{Q(L/2+t_f)cos\alpha}
     \bigg]
-    \end{align} 
+    \end{align}
 
 where
@@ -61,14 +61,14 @@
 
     R'=\frac{R}{\sqrt{2}}\sqrt{(1+X_{core}^{2}) + (1-X_{core}^{2})cos(\psi)}
-    
-    
-and $V_t = \pi.(R+t_r)(Xcore.R+t_r)^2.(L+2.t_f)$ is the total volume of the bicelle, 
-$V_c = \pi.Xcore.R^2.L$ the volume of the core, $V_{c+f} = \pi.Xcore.R^2.(L+2.t_f)$ 
+
+
+and $V_t = \pi.(R+t_r)(Xcore.R+t_r)^2.(L+2.t_f)$ is the total volume of the bicelle,
+$V_c = \pi.Xcore.R^2.L$ the volume of the core, $V_{c+f} = \pi.Xcore.R^2.(L+2.t_f)$
 the volume of the core plus the volume of the faces, $R$ is the radius
-of the core, $Xcore$ is the axial ratio of the core, $L$ the length of the core, 
-$t_f$ the thickness of the face, $t_r$ the thickness of the rim and $J_1$ the usual 
-first order bessel function. The core has radii $R$ and $Xcore.R$ so is circular, 
-as for the core_shell_bicelle model, for $Xcore$ =1. Note that you may need to 
-limit the range of $Xcore$, especially if using the Monte-Carlo algorithm, as 
+of the core, $Xcore$ is the axial ratio of the core, $L$ the length of the core,
+$t_f$ the thickness of the face, $t_r$ the thickness of the rim and $J_1$ the usual
+first order bessel function. The core has radii $R$ and $Xcore.R$ so is circular,
+as for the core_shell_bicelle model, for $Xcore$ =1. Note that you may need to
+limit the range of $Xcore$, especially if using the Monte-Carlo algorithm, as
 setting radius to $R/Xcore$ and axial ratio to $1/Xcore$ gives an equivalent solution!
 
@@ -76,5 +76,5 @@
 bicelles is then given by integrating over all possible $\alpha$ and $\psi$.
 
-For oriented bicelles the *theta*, *phi* and *psi* orientation parameters will appear when fitting 2D data, 
+For oriented bicelles the *theta*, *phi* and *psi* orientation parameters will appear when fitting 2D data,
 see the :ref:`elliptical-cylinder` model for further information.
 
@@ -82,5 +82,5 @@
 .. figure:: img/elliptical_cylinder_angle_definition.png
 
-    Definition of the angles for the oriented core_shell_bicelle_elliptical particles.   
+    Definition of the angles for the oriented core_shell_bicelle_elliptical particles.
 
 
@@ -105,6 +105,6 @@
 description = """
     core_shell_bicelle_elliptical
-    Elliptical cylinder core, optional shell on the two flat faces, and shell of 
-    uniform thickness on its rim (extending around the end faces).    
+    Elliptical cylinder core, optional shell on the two flat faces, and shell of
+    uniform thickness on its rim (extending around the end faces).
     Please see full documentation for equations and further details.
     Involves a double numerical integral around the ellipsoid diameter
@@ -136,17 +136,30 @@
           "core_shell_bicelle_elliptical.c"]
 
-demo = dict(scale=1, background=0,
-            radius=30.0,
-            x_core=3.0,
-            thick_rim=8.0,
-            thick_face=14.0,
-            length=50.0,
-            sld_core=4.0,
-            sld_face=7.0,
-            sld_rim=1.0,
-            sld_solvent=6.0,
-            theta=90,
-            phi=0,
-            psi=0)
+def random():
+    import numpy as np
+    outer_major = 10**np.random.uniform(1, 4.7)
+    outer_minor = 10**np.random.uniform(1, 4.7)
+    # Use a distribution with a preference for thin shell or thin core,
+    # limited by the minimum radius. Avoid core,shell radii < 1
+    min_radius = min(outer_major, outer_minor)
+    thick_rim = np.random.beta(0.5, 0.5)*(min_radius-2) + 1
+    radius_major = outer_major - thick_rim
+    radius_minor = outer_minor - thick_rim
+    radius = radius_major
+    x_core = radius_minor/radius_major
+    outer_length = 10**np.random.uniform(1, 4.7)
+    # Caps should be a small percentage of the total length, but at least one
+    # angstrom long.  Since outer length >= 10, the following will suffice
+    thick_face = 10**np.random.uniform(-np.log10(outer_length), -1)*outer_length
+    length = outer_length - thick_face
+    pars = dict(
+        radius=radius,
+        x_core=x_core,
+        thick_rim=thick_rim,
+        thick_face=thick_face,
+        length=length
+    )
+    return pars
+
 
 q = 0.1
Index: sasmodels/models/core_shell_cylinder.py
===================================================================
--- sasmodels/models/core_shell_cylinder.py	(revision 9b79f2966c4083a48b4894b6bf0587bef7a5be09)
+++ sasmodels/models/core_shell_cylinder.py	(revision 8f04da467b56e4e22291910958adebc417b57c19)
@@ -142,4 +142,19 @@
     return whole, whole - core
 
+def random():
+    import numpy as np
+    outer_radius = 10**np.random.uniform(1, 4.7)
+    # Use a distribution with a preference for thin shell or thin core
+    # Avoid core,shell radii < 1
+    radius = np.random.beta(0.5, 0.5)*(outer_radius-2) + 1
+    thickness = outer_radius - core
+    length = np.random.uniform(1, 4.7)
+    pars = dict(
+        radius=radius,
+        thickness=thickness,
+        length=length,
+    )
+    return pars
+
 demo = dict(scale=1, background=0,
             sld_core=6, sld_shell=8, sld_solvent=1,
Index: sasmodels/models/core_shell_ellipsoid.py
===================================================================
--- sasmodels/models/core_shell_ellipsoid.py	(revision 31df0c91cf57174acd8e9cbe5f9caadffdc880ba)
+++ sasmodels/models/core_shell_ellipsoid.py	(revision 8f04da467b56e4e22291910958adebc417b57c19)
@@ -154,10 +154,12 @@
     import numpy as np
     V = 10**np.random.uniform(5, 12)
-    radius_polar = 10**np.random.uniform(1.3, 4)
-    radius_equatorial = np.sqrt(V/radius_polar) # ignore 4/3 pi
-    thickness_polar = np.random.uniform(0.01, 1)*radius_polar
-    thickness_equatorial = np.random.uniform(0.01, 1)*radius_equatorial
-    radius_polar -= thickness_polar
-    radius_equatorial -= thickness_equatorial
+    outer_polar = 10**np.random.uniform(1.3, 4)
+    outer_equatorial = np.sqrt(V/radius_polar) # ignore 4/3 pi
+    # Use a distribution with a preference for thin shell or thin core
+    # Avoid core,shell radii < 1
+    thickness_polar = np.random.beta(0.5, 0.5)*(outer__polar-2) + 1
+    thickness_equatorial = np.random.beta(0.5, 0.5)*(outer_equatorial-2) + 1
+    radius_polar = outer_polar - thickness_polar
+    radius_equatorial = outer_equatorial - thickness_equatorial
     x_core = radius_polar/radius_equatorial
     x_polar_shell = thickness_polar/thickness_equatorial
Index: sasmodels/models/core_shell_parallelepiped.py
===================================================================
--- sasmodels/models/core_shell_parallelepiped.py	(revision 9b79f2966c4083a48b4894b6bf0587bef7a5be09)
+++ sasmodels/models/core_shell_parallelepiped.py	(revision 8f04da467b56e4e22291910958adebc417b57c19)
@@ -4,5 +4,5 @@
 
 Calculates the form factor for a rectangular solid with a core-shell structure.
-The thickness and the scattering length density of the shell or 
+The thickness and the scattering length density of the shell or
 "rim" can be different on each (pair) of faces. However at this time
 the 1D calculation does **NOT** actually calculate a c face rim despite the presence of
@@ -42,5 +42,5 @@
 
 **meaning that there are "gaps" at the corners of the solid.**  Again note that
-$t_C = 0$ currently. 
+$t_C = 0$ currently.
 
 The intensity calculated follows the :ref:`parallelepiped` model, with the
@@ -82,5 +82,5 @@
 based on the the averaged effective radius $(=\sqrt{(A+2t_A)(B+2t_B)/\pi})$
 and length $(C+2t_C)$ values, after appropriately
-sorting the three dimensions to give an oblate or prolate particle, to give an 
+sorting the three dimensions to give an oblate or prolate particle, to give an
 effective radius, for $S(Q)$ when $P(Q) * S(Q)$ is applied.
 
@@ -126,5 +126,5 @@
 title = "Rectangular solid with a core-shell structure."
 description = """
-     P(q)= 
+     P(q)=
 """
 category = "shape:parallelepiped"
@@ -179,4 +179,19 @@
 # VR defaults to 1.0
 
+def random():
+    import numpy as np
+    outer = 10**np.random.uniform(1, 4.7, size=3)
+    thick = np.random.beta(0.5, 0.5, size=3)*(outer-2) + 1
+    length = outer - thick
+    pars = dict(
+        length_a=length[0],
+        length_b=length[1],
+        length_c=length[2],
+        thick_rim_a=thick[0],
+        thick_rim_b=thick[1],
+        thick_rim_c=thick[2],
+    )
+    return pars
+
 # parameters for demo
 demo = dict(scale=1, background=0.0,
Index: sasmodels/models/core_shell_sphere.py
===================================================================
--- sasmodels/models/core_shell_sphere.py	(revision 31df0c91cf57174acd8e9cbe5f9caadffdc880ba)
+++ sasmodels/models/core_shell_sphere.py	(revision 8f04da467b56e4e22291910958adebc417b57c19)
@@ -101,7 +101,9 @@
 def random():
     import numpy as np
-    total_radius = 10**np.random.uniform(1.3, 4.3)
-    radius = np.random.uniform(0, 1)*total_radius
-    thickness = total_radius - radius
+    outer_radius = 10**np.random.uniform(1.3, 4.3)
+    # Use a distribution with a preference for thin shell or thin core
+    # Avoid core,shell radii < 1
+    radius = np.random.beta(0.5, 0.5)*(outer_radius-2) + 1
+    thickness = outer_radius - core
     pars = dict(
         radius=radius,
@@ -112,11 +114,10 @@
 tests = [
     [{'radius': 20.0, 'thickness': 10.0}, 'ER', 30.0],
-     # TODO: VR test suppressed until we sort out new product model
-     # and determine what to do with volume ratio.
-     #[{'radius': 20.0, 'thickness': 10.0}, 'VR', 0.703703704],
+    # TODO: VR test suppressed until we sort out new product model
+    # and determine what to do with volume ratio.
+    #[{'radius': 20.0, 'thickness': 10.0}, 'VR', 0.703703704],
 
-     # The SasView test result was 0.00169, with a background of 0.001
-     [{'radius': 60.0, 'thickness': 10.0, 'sld_core': 1.0, 'sld_shell':2.0,
-       'sld_solvent':3.0, 'background':0.0},
-      0.4, 0.000698838],
+    # The SasView test result was 0.00169, with a background of 0.001
+    [{'radius': 60.0, 'thickness': 10.0, 'sld_core': 1.0, 'sld_shell': 2.0,
+      'sld_solvent': 3.0, 'background': 0.0}, 0.4, 0.000698838],
 ]
Index: sasmodels/models/fcc_paracrystal.py
===================================================================
--- sasmodels/models/fcc_paracrystal.py	(revision 1511c37c7ea796eff9bd41e33beabb6efce39d91)
+++ sasmodels/models/fcc_paracrystal.py	(revision 8f04da467b56e4e22291910958adebc417b57c19)
@@ -125,5 +125,5 @@
     d_factor = 10**np.random.uniform(-2, -0.7)  # sigma_d in 0.01-0.7
     dnn_fraction = np.random.beta(a=10, b=1)
-    dnn = radius*4/np.sqrt(3)/dnn_fraction
+    dnn = radius*4/np.sqrt(2)/dnn_fraction
     pars = dict(
         #sld=1, sld_solvent=0, scale=1, background=1e-32,
@@ -135,9 +135,8 @@
 
 # april 10 2017, rkh add unit tests, NOT compared with any other calc method, assume correct!
-q =4.*pi/220.
+q = 4.*pi/220.
 tests = [
-    [{ },
-     [0.001, q, 0.215268], [0.275164706668, 5.7776842567, 0.00958167119232]],
-     [{}, (-0.047, -0.007), 238.103096286 ],
-     [{}, (0.053, 0.063), 0.863609587796 ],
+    [{}, [0.001, q, 0.215268], [0.275164706668, 5.7776842567, 0.00958167119232]],
+    [{}, (-0.047, -0.007), 238.103096286],
+    [{}, (0.053, 0.063), 0.863609587796],
 ]
Index: sasmodels/models/fractal_core_shell.py
===================================================================
--- sasmodels/models/fractal_core_shell.py	(revision 1511c37c7ea796eff9bd41e33beabb6efce39d91)
+++ sasmodels/models/fractal_core_shell.py	(revision 8f04da467b56e4e22291910958adebc417b57c19)
@@ -100,15 +100,12 @@
 def random():
     import numpy as np
-    total_radius = 10**np.random.uniform(0.7, 4)
-    core_portion = np.random.uniform(0, 1)
-    radius = total_radius * core_portion
-    thickness = total_radius - radius
-    #radius = 5
-    cor_length = 10**np.random.uniform(0.7, 2)*radius
-    #cor_length = 20*radius
+    outer_radius = 10**np.random.uniform(0.7, 4)
+    # Use a distribution with a preference for thin shell or thin core
+    # Avoid core,shell radii < 1
+    thickness = np.random.beta(0.5, 0.5)*(outer_radius-2) + 1
+    radius = outer_radius - thickness
+    cor_length = 10**np.random.uniform(0.7, 2)*outer_radius
     volfraction = 10**np.random.uniform(-3, -1)
-    #volfraction = 0.05
     fractal_dim = 2*np.random.beta(3, 4) + 1
-    #fractal_dim = 2
     pars = dict(
         #background=0, sld_block=1, sld_solvent=0,
@@ -119,5 +116,4 @@
     )
     return pars
-
 
 demo = dict(scale=0.05,
Index: sasmodels/models/hardsphere.py
===================================================================
--- sasmodels/models/hardsphere.py	(revision 40a87fa5096610e24c3ead10299f609e99f093f2)
+++ sasmodels/models/hardsphere.py	(revision 8f04da467b56e4e22291910958adebc417b57c19)
@@ -75,5 +75,5 @@
 Iq = r"""
       double D,A,B,G,X,X2,X4,S,C,FF,HARDSPH;
-      // these are c compiler instructions, can also put normal code inside the "if else" structure 
+      // these are c compiler instructions, can also put normal code inside the "if else" structure
       #if FLOAT_SIZE > 4
       // double precision    orig had 0.2, don't call the variable cutoff as PAK already has one called that! Must use UPPERCASE name please.
@@ -82,5 +82,5 @@
       #else
       // 0.1 bad, 0.2 OK, 0.3 good, 0.4 better, 0.8 no good
-      #define CUTOFFHS 0.4  
+      #define CUTOFFHS 0.4
       #endif
 
@@ -110,7 +110,7 @@
       if(X < CUTOFFHS) {
       // RKH Feb 2016, use Taylor series expansion for small X
-      // else no obvious way to rearrange the equations to avoid needing a very high number of significant figures. 
-      // Series expansion found using Mathematica software. Numerical test in .xls showed terms to X^2 are sufficient 
-      // for 5 or 6 significant figures, but I put the X^4 one in anyway 
+      // else no obvious way to rearrange the equations to avoid needing a very high number of significant figures.
+      // Series expansion found using Mathematica software. Numerical test in .xls showed terms to X^2 are sufficient
+      // for 5 or 6 significant figures, but I put the X^4 one in anyway
             //FF = 8*A +6*B + 4*G - (0.8*A +2.0*B/3.0 +0.5*G)*X2 +(A/35. +B/40. +G/50.)*X4;
             // refactoring the polynomial makes it very slightly faster (0.5 not 0.6 msec)
@@ -153,4 +153,13 @@
    """
 
+def random():
+    import numpy as np
+    pars = dict(
+        scale=1, background=0,
+        radius_effective=10**np.random.uniform(1, 4),
+        volfraction=10**np.random.uniform(-2, 0),  # high volume fraction
+    )
+    return pars
+
 # ER defaults to 0.0
 # VR defaults to 1.0
Index: sasmodels/models/hayter_msa.py
===================================================================
--- sasmodels/models/hayter_msa.py	(revision 5702f5291bd20fe3c0e9feff9bbfc2b4f621f306)
+++ sasmodels/models/hayter_msa.py	(revision 8f04da467b56e4e22291910958adebc417b57c19)
@@ -64,5 +64,5 @@
         Routine takes absolute value of charge, use HardSphere if charge
         goes to zero.
-        In sasview the effective radius and volume fraction may be calculated 
+        In sasview the effective radius and volume fraction may be calculated
         from the parameters used in P(Q).
 """
@@ -89,4 +89,22 @@
 # ER defaults to 0.0
 # VR defaults to 1.0
+
+def random():
+    import numpy as np
+    # TODO: too many failures for random hayter_msa parameters
+    pars = dict(
+        scale=1, background=0,
+        radius_effective=10**np.random.uniform(1, 4.7),
+        volfraction=10**np.random.uniform(-2, 0),  # high volume fraction
+        charge=min(int(10**np.random.uniform(0, 1.3)+0.5), 200),
+        temperature=10**np.random.uniform(0, np.log10(450)), # max T = 450
+        #concentration_salt=10**np.random.uniform(-3, 1),
+        dialectconst=10**np.random.uniform(0, 6),
+        #charge=10,
+        #temperature=318.16,
+        concentration_salt=0.0,
+        #dielectconst=71.08,
+    )
+    return pars
 
 # default parameter set,  use  compare.sh -midQ -linear
Index: sasmodels/models/hollow_cylinder.py
===================================================================
--- sasmodels/models/hollow_cylinder.py	(revision 31df0c91cf57174acd8e9cbe5f9caadffdc880ba)
+++ sasmodels/models/hollow_cylinder.py	(revision 8f04da467b56e4e22291910958adebc417b57c19)
@@ -123,14 +123,25 @@
 def random():
     import numpy as np
-    length = 10**np.random.uniform(2, 6)
-    radius = 10**np.random.uniform(1, 3)
-    kuhn_length = 10**np.random.uniform(-2, -0.7)*length  # at least 10 segments
+    length = 10**np.random.uniform(1, 4.7)
+    outer_radius = 10**np.random.uniform(1, 4.7)
+    # Use a distribution with a preference for thin shell or thin core
+    # Avoid core,shell radii < 1
+    thickness = np.random.beta(0.5, 0.5)*(outer_radius-2) + 1
+    radius = outer_radius - thickness
     pars = dict(
         length=length,
         radius=radius,
-        kuhn_length=kuhn_length,
+        thickness=thickness,
     )
     return pars
 
+# parameters for demo
+demo = dict(scale=1.0, background=0.0, length=400.0, radius=20.0,
+            thickness=10, sld=6.3, sld_solvent=1, theta=90, phi=0,
+            thickness_pd=0.2, thickness_pd_n=9,
+            length_pd=.2, length_pd_n=10,
+            radius_pd=.2, radius_pd_n=9,
+            theta_pd=10, theta_pd_n=5,
+           )
 q = 0.1
 # april 6 2017, rkh added a 2d unit test, assume correct!
Index: sasmodels/models/hollow_rectangular_prism.py
===================================================================
--- sasmodels/models/hollow_rectangular_prism.py	(revision 31df0c91cf57174acd8e9cbe5f9caadffdc880ba)
+++ sasmodels/models/hollow_rectangular_prism.py	(revision 8f04da467b56e4e22291910958adebc417b57c19)
@@ -149,5 +149,10 @@
     import numpy as np
     a, b, c = 10**np.random.uniform(1, 4.7, size=3)
-    thickness = np.random.uniform(0.01, 0.49) * min(a, b, c)
+    # Thickness is limited to 1/2 the smallest dimension
+    # Use a distribution with a preference for thin shell or thin core
+    # Avoid core,shell radii < 1
+    min_dim = 0.5*min(a, b, c)
+    thickness = np.random.beta(0.5, 0.5)*(min_dim-2) + 1
+    #print(a, b, c, thickness, thickness/min_dim)
     pars = dict(
         length_a=a,
Index: sasmodels/models/lamellar_stack_paracrystal.py
===================================================================
--- sasmodels/models/lamellar_stack_paracrystal.py	(revision 1511c37c7ea796eff9bd41e33beabb6efce39d91)
+++ sasmodels/models/lamellar_stack_paracrystal.py	(revision 8f04da467b56e4e22291910958adebc417b57c19)
@@ -140,7 +140,5 @@
     # Let polydispersity peak around 15%; 95% < 0.4; max=100%
     sigma_d = np.random.beta(1.5, 7)
-    scale = thickness * 10**np.random.uniform(-7, -4)
     pars = dict(
-        #scale=1,
         thickness=thickness,
         Nlayers=Nlayers,
Index: sasmodels/models/linear_pearls.py
===================================================================
--- sasmodels/models/linear_pearls.py	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
+++ sasmodels/models/linear_pearls.py	(revision 8f04da467b56e4e22291910958adebc417b57c19)
@@ -65,10 +65,15 @@
 source = ["lib/sas_3j1x_x.c", "linear_pearls.c"]
 
-demo = dict(scale=1.0, background=0.0,
-            radius=80.0,
-            edge_sep=350.0,
-            num_pearls=3,
-            sld=1.0,
-            sld_solvent=6.3)
+def random():
+    import numpy as np
+    radius = 10**np.random.uniform(1, 3) # 1 - 1000
+    edge_sep = 10**np.random.uniform(0, 3)  # 1 - 1000
+    num_pearls = np.round(10**np.random.uniform(0.3, 3)) # 2 - 1000
+    pars = dict(
+        radius=radius,
+        edge_sep=edge_sep,
+        num_pearls=num_pearls,
+    )
+    return pars
 
 """
Index: sasmodels/models/multilayer_vesicle.py
===================================================================
--- sasmodels/models/multilayer_vesicle.py	(revision 5d23de23656b990e141f01721eb85cd0afaa4dc2)
+++ sasmodels/models/multilayer_vesicle.py	(revision 8f04da467b56e4e22291910958adebc417b57c19)
@@ -149,12 +149,25 @@
     return radius + n_shells * (thick_shell + thick_solvent) - thick_solvent
 
-demo = dict(scale=1, background=0,
-            volfraction=0.05,
-            radius=60.0,
-            thick_shell=10.0,
-            thick_solvent=10.0,
-            sld_solvent=6.4,
-            sld=0.4,
-            n_shells=2.0)
+def random():
+    import numpy as np
+    volfraction = 10**np.random.uniform(-3, -0.5)  # scale from 0.1% to 30%
+    radius = 10**np.random.uniform(0, 2.5) # core less than 300 A
+    total_thick = 10**np.random.uniform(2, 4) # up to 10000 A of shells
+    # random number of shells, with shell+solvent thickness > 10 A
+    n_shells = int(10**np.random.uniform(0, np.log10(total_thick)-1)+0.5)
+    # split total shell thickness with preference for shell over solvent;
+    # make sure that shell thickness is at least 1 A
+    one_thick = total_thick/n_shells
+    thick_solvent = 10**np.random.uniform(-2, 0)*(one_thick - 1)
+    thick_shell = one_thick - thick_solvent
+    pars = dict(
+        scale=1,
+        volfraction=volfraction,
+        radius=radius,
+        thick_shell=thick_shell,
+        thick_solvent=thick_solvent,
+        n_shells=n_shells,
+    )
+    return pars
 
 tests = [
Index: sasmodels/models/parallelepiped.py
===================================================================
--- sasmodels/models/parallelepiped.py	(revision 31df0c91cf57174acd8e9cbe5f9caadffdc880ba)
+++ sasmodels/models/parallelepiped.py	(revision 8f04da467b56e4e22291910958adebc417b57c19)
@@ -245,9 +245,9 @@
 def random():
     import numpy as np
-    a, b, c = 10**np.random.uniform(1, 4.7, size=3)
+    length = 10**np.random.uniform(1, 4.7, size=3)
     pars = dict(
-        length_a=a,
-        length_b=b,
-        length_c=c,
+        length_a=length[0],
+        length_b=length[1],
+        length_c=length[2],
     )
     return pars
Index: sasmodels/models/pearl_necklace.py
===================================================================
--- sasmodels/models/pearl_necklace.py	(revision 1bd1ea22a4ccb6ec87819647e470bbc002b2a1ac)
+++ sasmodels/models/pearl_necklace.py	(revision 8f04da467b56e4e22291910958adebc417b57c19)
@@ -117,4 +117,18 @@
     return rad_out
 
+def random():
+    import numpy as np
+    radius = 10**np.random.uniform(1, 3) # 1 - 1000
+    thick_string = 10**np.random.uniform(0, np.log10(radius)-1) # 1 - radius/10
+    edge_sep = 10**np.random.uniform(0, 3)  # 1 - 1000
+    num_pearls = np.round(10**np.random.uniform(0.3, 3)) # 2 - 1000
+    pars = dict(
+        radius=radius,
+        edge_sep=edge_sep,
+        thick_string=thick_string,
+        num_pearls=num_pearls,
+    )
+    return pars
+
 # parameters for demo
 demo = dict(scale=1, background=0, radius=80.0, edge_sep=350.0,
Index: sasmodels/models/pringle.py
===================================================================
--- sasmodels/models/pringle.py	(revision 30fbe2eda4d6f97acc840bd5e605de3829caf3c7)
+++ sasmodels/models/pringle.py	(revision 8f04da467b56e4e22291910958adebc417b57c19)
@@ -85,12 +85,16 @@
     return 0.5 * (ddd) ** (1. / 3.)
 
-demo = dict(background=0.0,
-            scale=1.0,
-            radius=60.0,
-            thickness=10.0,
-            alpha=0.001,
-            beta=0.02,
-            sld=1.0,
-            sld_solvent=6.35)
+def random():
+    import numpy as np
+    alpha, beta = 10**np.random.uniform(-1, 1, size=2)
+    radius = 10**np.random.uniform(1, 3)
+    thickness = 10**np.random.uniform(0.7, 2)
+    pars = dict(
+        radius=radius,
+        thickness=thickness,
+        alpha=alpha,
+        beta=beta,
+    )
+    return pars
 
 tests = [
Index: sasmodels/models/raspberry.py
===================================================================
--- sasmodels/models/raspberry.py	(revision 8e68ea06a7246c4bd95ce82c755719ced6fa5a60)
+++ sasmodels/models/raspberry.py	(revision 8f04da467b56e4e22291910958adebc417b57c19)
@@ -154,4 +154,24 @@
 source = ["lib/sas_3j1x_x.c", "raspberry.c"]
 
+def random():
+    import numpy as np
+    # Limit volume fraction to 20% each
+    volfraction_lg = 10**np.random.uniform(-3, -0.3)
+    volfraction_sm = 10**np.random.uniform(-3, -0.3)
+    # Prefer most particles attached (peak near 60%), but not all or none
+    surface_fraction = np.random.beta(6, 4)
+    radius_lg = 10**np.random.uniform(1.7, 4.7)  # 500 - 50000 A
+    radius_sm = 10**np.random.uniform(-3, -0.3)*radius_lg  # 0.1% - 20%
+    penetration = np.random.beta(1, 10) # up to 20% pen. for 90% of examples
+    pars = dict(
+        volfraction_lg=volfraction_lg,
+        volfraction_sm=volfraction_sm,
+        surface_fraction=surface_fraction,
+        radius_lg=radius_lg,
+        radius_sm=radius_sm,
+        penetration=penetration,
+    )
+    return pars
+
 # parameters for demo
 demo = dict(scale=1, background=0.001,
Index: sasmodels/models/sc_paracrystal.py
===================================================================
--- sasmodels/models/sc_paracrystal.py	(revision 1511c37c7ea796eff9bd41e33beabb6efce39d91)
+++ sasmodels/models/sc_paracrystal.py	(revision 8f04da467b56e4e22291910958adebc417b57c19)
@@ -144,5 +144,5 @@
     d_factor = 10**np.random.uniform(-2, -0.7)  # sigma_d in 0.01-0.7
     dnn_fraction = np.random.beta(a=10, b=1)
-    dnn = radius*4/np.sqrt(3)/dnn_fraction
+    dnn = radius*4/np.sqrt(4)/dnn_fraction
     pars = dict(
         #sld=1, sld_solvent=0, scale=1, background=1e-32,
Index: sasmodels/models/squarewell.py
===================================================================
--- sasmodels/models/squarewell.py	(revision 3330bb44faf647c84e605c5a75bed3f385a25a52)
+++ sasmodels/models/squarewell.py	(revision 8f04da467b56e4e22291910958adebc417b57c19)
@@ -91,10 +91,10 @@
 	double S,C,SL,CL;
 	x= q;
-	
+
 	req = radius_effective;
 	phis = volfraction;
 	edibkb = welldepth;
 	lambda = wellwidth;
-	
+
 	sigma = req*2.;
 	eta = phis;
@@ -110,7 +110,7 @@
 	beta = -(eta/3.0) * ( 18. + 20.*eta - 12.*eta2 + eta4 )/etam14;
 	gamm = 0.5*eta*( sk + eta3*(eta-4.) )/etam14;
-	
+
 	//  calculate the structure factor
-	
+
 	sk = x*sigma;
 	sk2 = sk*sk;
@@ -125,5 +125,5 @@
 	ck =  -24.0*eta*( t1 + t2 + t3 + t4 )/sk3/sk3;
 	struc  = 1./(1.-ck);
-	
+
 	return(struc);
 """
@@ -132,12 +132,22 @@
 # VR defaults to 1.0
 
+def random():
+    import numpy as np
+    pars = dict(
+        scale=1, background=0,
+        radius_effective=10**np.random.uniform(1, 4.7),
+        volfraction=np.random.uniform(0.00001, 0.08),
+        welldepth=np.random.uniform(0, 1.5),
+        wellwidth=np.random.uniform(1, 1.2),
+    )
+    return pars
+
 demo = dict(radius_effective=50, volfraction=0.04, welldepth=1.5,
             wellwidth=1.2, radius_effective_pd=0, radius_effective_pd_n=0)
 #
 tests = [
-    [{'scale': 1.0, 'background' : 0.0, 'radius_effective' : 50.0,
-      'volfraction' : 0.04,'welldepth' : 1.5, 'wellwidth' : 1.2,
-      'radius_effective_pd' : 0},
-     [0.001], [0.97665742]],
+    [{'scale': 1.0, 'background': 0.0, 'radius_effective': 50.0,
+      'volfraction': 0.04,'welldepth': 1.5, 'wellwidth': 1.2,
+      'radius_effective_pd': 0}, [0.001], [0.97665742]],
     ]
 # ADDED by: converting from sasview RKH  ON: 16Mar2016
Index: sasmodels/models/stacked_disks.py
===================================================================
--- sasmodels/models/stacked_disks.py	(revision 9d50a1e22d29f2f0e48aa22713dc40b14277b429)
+++ sasmodels/models/stacked_disks.py	(revision 8f04da467b56e4e22291910958adebc417b57c19)
@@ -137,4 +137,23 @@
 
 source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "stacked_disks.c"]
+
+def random():
+    import numpy as np
+    radius = 10**np.random.uniform(1, 4.7)
+    total_stack = 10**np.random.uniform(1, 4.7)
+    n_stacking = int(10**np.random.uniform(0, np.log10(total_stack)-1) + 0.5)
+    d = total_stack/n_stacking
+    thick_core = np.random.uniform(0, d-2)  # at least 1 A for each layer
+    thick_layer = (d - thick_core)/2
+    # Let polydispersity peak around 15%; 95% < 0.4; max=100%
+    sigma_d = d * np.random.beta(1.5, 7)
+    pars = dict(
+        thick_core=thick_core,
+        thick_layer=thick_layer,
+        radius=radius,
+        n_stacking=n_stacking,
+        sigma_d=sigma_d,
+    )
+    return pars
 
 demo = dict(background=0.001,
Index: sasmodels/models/stickyhardsphere.py
===================================================================
--- sasmodels/models/stickyhardsphere.py	(revision 40a87fa5096610e24c3ead10299f609e99f093f2)
+++ sasmodels/models/stickyhardsphere.py	(revision 8f04da467b56e4e22291910958adebc417b57c19)
@@ -98,4 +98,15 @@
     ]
 
+def random():
+    import numpy as np
+    pars = dict(
+        scale=1, background=0,
+        radius_effective=10**np.random.uniform(1, 4.7),
+        volfraction=np.random.uniform(0.00001, 0.74),
+        perturb=10**np.random.uniform(-2, -1),
+        stickiness=np.random.uniform(0, 1),
+    )
+    return pars
+
 # No volume normalization despite having a volume parameter
 # This should perhaps be volume normalized?
Index: sasmodels/product.py
===================================================================
--- sasmodels/product.py	(revision f88e2480698b5dcd90f16c16ecd9d445273a1b9a)
+++ sasmodels/product.py	(revision 8f04da467b56e4e22291910958adebc417b57c19)
@@ -93,4 +93,6 @@
     # Remember the component info blocks so we can build the model
     model_info.composition = ('product', [p_info, s_info])
+    # TODO: delegate random to p_info, s_info
+    #model_info.random = lambda: {}
     model_info.demo = demo
 
